1. 서 론
최근 재진입 비행체 기술 발전에 따라, 이에 대응하기 위한 탄도탄 방어체계 고도화 연구가 심도 있게 진행되고 있다. MaRV(maneuvering
reentry vehicle)는 탄도비행 궤적을 그리는 전통적인 BRV(ballistic reentry vehicle)와 달리 재진입 단계에서 변칙적인
상승(pull-up) 기동을 통해 탐지·추적 및 요격을 회피할 수 있어 대공 방어체계에 큰 위협이 되고 있다[1,2]. MaRV에 효과적으로 대응하기 위해서는 재진입 단계에서의 궤적을 장시간 예측한 후, 이를 기반으로 요격탄을 운용해야한다. 하지만, MaRV의 기동
시점과 공력학적 특성을 정확히 알 수 없으므로 MaRV의 상승기동 예측은 아직까지 해결이 매우 어려운 도전적 과제로 인식되고 있다[2].
지금까지 진행되어온 재진입 비행체의 궤적 예측 연구는 모델 기반(model-based) 기법과 데이터 기반(data-driven) 기법으로 분류할
수 있다[3]. 모델 기반 기법은 재진입 비행체의 비선형 운동방정식과 추적필터에서 제공되는 현시점 항적을 이용해 미래 궤적을 산출하는 소위 IVP(initial
value problem)를 푸는 방식을 취하고 있다. 여기서 핵심적인 기술적 난제는 중력과 항력뿐만 아니라 횡기동 수행 시 발생하는 양력을 포괄하는
정교한 운동모델의 구축에 있다[4-7]. 재진입 비행체에 작용하는 공력학적 특성은 속력 및 받음각의 비선형 함수로 모델링되는데, 비행 고도가 낮아짐에 따라 궤적 예측에 영향을 끼치는 공력학적
특성이 크게 변화한다. 이를 고려하기 위해 비행조건별 공력 데이터베이스를 구축하여 궤적 예측에 활용하는 방안이 제시되었다[6-8]. 하지만 이 방법은 지나치게 많은 공력 데이터베이스를 사전지식으로 요구하므로 재진입 비행체의 공력학적 형상 불확실성에 매우 취약하다는 구조적 단점을
지니고 있다.
그 대안으로 물리적 모델에 의존하지 않고 실 데이터로부터 직접 재진입 비행체의 비행패턴과 동특성을 학습하는 소위 데이터 기반 궤적 예측 기법이 연구되었다.
가장 대표적인 사례 중 하나인 GP(Gaussian process) 기법은 커널 함수를 통해 비행체의 비선형 동특성을 비모수적으로 학습하는 방식을
취함으로써 운동모델 불확실성이 장시간 궤적 예측 성능에 끼치는 영향을 완화하고자 했다[3,9]. 최근에는 복잡한 기동패턴을 보이는 HGV(hypersonic gliding vehicle)의 궤적 예측을 위해 RNN(recurrent neural
network)이나 LSTM(long short-term memory)을 활용한 궤적 예측 방법도 제시된 바 있다[10-16]. 하지만, 이들 데이터 기반 궤적 예측 기법은 학습을 위한 재진입 비행체의 궤적 데이터를 충분히 확보하기 어려운 실제 상황에서는 만족할 만한 궤적
예측 성능을 제공하기 어렵다. 더욱이 학습데이터에 의존적인 궤적 예측성능, 과적합으로 인한 궤적 예측의 안정성 저하 문제 등 실제 응용을 가로막는
여러 기술적 이슈를 내포하고 있다[10,11].
이러한 한계를 극복하기 위해 본 논문에서는 MaRV의 상승기동 특성을 활용한 궤적 예측 기법을 새롭게 제안한다. 제안한 알고리듬은 그림 1과 같이 상승기동 탐지 및 궤적 예측의 두 부분으로 구성된다. MaRV는 기준 고도 이하로 진입하는 $Q_{i}$ 지점에서부터 미지의 상승기동을 수행하므로,
MaRV의 기동을 효과적으로 탐지하기 위해 탄도비행 및 상승기동 구간에서의 운동모델에 해당하는 상태추정기들을 설계하여 필터 뱅크를 구성한다. 잘 알려진
바와 같이 운동모델의 불확실성은 필터 추정치의 편향오차를 야기하므로, 탄도비행 운동모델을 채택한 부필터의 잔차에 대한 카이제곱 검정을 통해 손쉽게
MaRV의 상승기동 시점을 결정할 수 있다. 일단 기동이 탐지되면 MaRV가 사거리 연장을 위해 공력학적 효율성이 확보되는 양항비를 유지한다는 점에
착안하여[17], $Q_{i}\sim Q_{o}$ 상승기동 구간에서의 궤적 특성을 분석한다. 상승기동에 의해 MaRV가 대기권계면을 재이탈하는 경우 대기 밀도가
0이 되므로, 이를 경계조건으로 삼아 MaRV 궤적의 닫힌 해(closed-form solution)를 유도한다. 이러한 수학적 해석결과를 활용할
경우, 제한된 관측정보와 기동탐지 시점에서의 양항비 추정치만으로도 운동모델의 영향을 크게 받는 대기권 내에서의 상승기동 궤적을 효과적으로 예측할 수
있다. 즉, 제안한 기법은 운동모델 불확실성이 존재하는 상황에서도 안정적인 궤적 예측 결과를 제공할 수 있다는 점에서 주어진 사전지식에 전적으로 의존하는
기존 모델기반 궤적 예측 기법과 명확히 대비된다. 더욱이, 제안 기법은 MaRV의 상승기동 궤적 예측과정에서 수치적분의 비중이 낮아 기존 기법 대비
연산량 부담이 적다는 장점이 있다. 전형적인 MaRV 궤적 시나리오에 대한 시뮬레이션을 통해 제안 기법이 기존 기법에 비해 우수한 상승기동 궤적 예측
성능을 보임을 확인한다.
그림 1. 재진입 비행체의 비행궤적
Fig. 1. Typical trajectory of a MaRV
2. 재진입 비행체의 상승기동 탐지
본 절에서는 MaRV의 비행구간 별 운동방정식을 수립하고 이를 바탕으로 MaRV의 상승기동 탐지 알고리듬을 설계한다.
2.1 상승기동 여부에 따른 운동방정식
MaRV의 운동을 기술하기 위해 도입된 좌표계들의 정의는 다음과 같다.
• 지구좌표계 (E-frame)
원점이 지구 중심에 위치하고 $X_{E}$축은 본초자오선과 적도면의 교점 방향, $Z_{E}$축은 지구 자전축과 일치하는 오른손 좌표계 이다. 지구좌표계
상에서 MaRV의 위치벡터와 속도벡터는 각각 $ R^{E}=[xyz]^{T}$와 $ V^{E}=[\dot{x}\dot{y}\dot{z}]^{T}$이다.
• 항법좌표계 (N-frame)
MaRV의 위도($L$)와 경도($l$)을 원점으로 하는 국지수평 좌표계로, $X_{N}$축은 정북 방향, $Y_{N}$축은 정동 방향, $Z_{N}$축은
지구 중심 방향으로 정의된다. 지구좌표계와 항법좌표계간의 좌표변환 행렬은 식 (1)로 계산된다.
• 속도좌표계 (V-frame)
원점은 MaRV의 무게중심, $X_{V}$축은 속도벡터 방향, $Y_{V}$축은 항법좌표계의 $X_{N}- Y_{N}$ 평면에 투영된 속도벡터의 우측
방향, $Z_{V}$축은 $X_{V}- Y_{V}$ 평면에 수직인 방향으로 정의되는 오른손 좌표계이다. 속도좌표계와 항법좌표계의 관계는 그림 2에 도시한 바와 같다.
그림 2. 속도 좌표계 정의
Fig. 2. Definition of V-frame
항법좌표계 상에서 표현된 MaRV 속도벡터를 $ V^{N}=[\dot{x}^{N}\dot{y}^{N}\dot{z}^{N}]^{T}$라 하면 속도좌표계와
항법좌표계 간의 좌표변환 행렬은 다음과 같이 정의된다.
여기서 $v=\left . ∥ V^{N}\right .∥$와 $v_{g}=\sqrt{\left(\dot{x}^{N}\right)^{2}+\left(\dot{y}^{N}\right)^{2}}$는
각각 MaRV의 속력과 지면 속력을 의미한다.
• 동체좌표계 (B-frame)
MaRV 동체에 고정된 좌표계로 원점은 MaRV 무게중심에 위치하며 $X_{B}$축은 기수방향, $Y_{B}$축은 우측 날개방향과 일치하는 오른손
좌표계이다. 참고로 그림 3에 도시한 바와 같이 속도좌표계와 동체좌표계 간의 관계는 수직/수평 받음각 $\alpha$, $\beta$ 및 총 받음각 $\alpha_{t}\approx\sqrt{\alpha^{2}+\beta^{2}}$로
설명된다.
그림 3. 동체 좌표계 정의
Fig. 3. Definition of body coordinate systems
재진입 단계에서 MaRV에 작용하는 외력은 공력과 중력이다. 만일, 지구 자전에 의한 코리올리 효과가 크지 않다면 지구좌표계는 MaRV의 운동을 지구좌표계
상에서 식 (3)과 같이 쓸 수 있다[2].
여기서 $q(h,\: v)=0.5\rho v^{2}$는 동압(dynamic pressure), $(\alpha_{D},\: \alpha_{T},\:
\alpha_{C})$는 공력에 의한 기동 파라미터, $\rho$는 공기밀도, $GM$은 표준중력상수(standard gravitational parameter)
이다. 기동 파라미터는 무차원 공력계수와 다음 관계를 갖는다[2].
식 (4)에서 $m$, $S_{ref}$는 MaRV의 질량과 단면적, $C_{D}$와 $C_{L}$은 MaRV의 형상에 의해 결정되는 항력/양력 계수를 나타낸다.
불행하게도 실제 상황에서 MaRV의 제원 및 공력계수를 정확히 알 수 없으므로, 본 논문에서는 가용정보를 이용한 기동 파라미터 추정 기법을 설계하고,
그 결과를 활용하여 궤적 예측을 수행한다.
대부분의 MaRV는 십자 대칭형상(cruciform)을 가지므로, 양력과 항력의 상관 관계를 다음 식으로 모델링할 수 있다[18].
여기서 $C_{D0}$는 양력이 없을 경우의 기저 항력(base drag), $K$는 유도항력계수(induced drag parameter), $n$은
MaRV의 형상에 따라 달라지는 양항곡선의 지수 항이다[19].
Remark 1. 식 (5)에서 확인할 수 있듯이, MaRV가 재진입 후 상승기동 궤적을 그리기 위해 양력을 발생시키면, 불가피하게 유도항력에 의한 감속이 야기된다. 이 경우
탄도계수(ballistic coefficient), 즉 항력에 해당하는 기동파라미터 $\alpha_{D}$만으로는 MaRV의 비행특성을 제대로 기술할
수 없으며, 3축 기동파라미터를 모두 고려해야 한다[2].
식 (4)로부터 재진입시 MaRV에 작용하는 항력 및 양력은 다음과 같이 계산된다.
이때 공기밀도 $\rho$은 고도 $h$의 함수로 근사할 수 있다.
여기서 $\rho_{0}$는 해수면에서의 공기밀도, $k_{\beta}$는 감쇠계수이다.
식 (5)를 식 (6)에 대입하여 정리하면 양항비 $\eta$를 무차원 공력계수 $C_{D_{0}}$와 $C_{L}$로 쓸 수 있다.
Remark 2. 식 (8)에서 MaRV의 비행거리를 최대화하기 위해서는 상승기동을 수행해야 하지만, 이로 인해 지나치게 큰 유도항력이 발생하는 경우 속력 감소로 인해 기대한
만큼의 비행거리 증대 효과를 거둘 수 없을 것임을 유추해 볼 수 있다. 따라서, 상승기동 구간에서는 MaRV가 공력학적 효율을 따져 적절한 양항비를
유지하는 것이 바람직하다.
2.2 상승기동 탐지 알고리듬 설계
MaRV의 운동은 횡기동의 수행여부에 따라 크게 달라지므로 상승기동 예측을 위해 먼저 비행 구간별 운동모델이 반영된 궤적추정기들로 필터 뱅크를 구성한
후, 각 부필터의 잔차 특성을 활용해 상승기동 여부를 탐지할 필요가 있다.
• 탄도비행 운동모델 (BMM: ballistic motion model)
양력이 작용하지 않는 경우($\alpha_{T}=\alpha_{C}= 0$), MaRV의 비행궤적은 탄도비행 운동모델로 모사된다. MaRV의 고도와
속력에 따른 기동 파라미터 $\alpha_{D}$의 변화는 급격하지 않으므로 이를 $\dot{\alpha}_{D}\approx 0$로 근사해도 무방하다.
이 경우, 식 (3)은 다음 7차 비선형 미분방정식으로 간략화 된다.
여기서
샘플링주기 $T$로 연속시간 비선형 미분방정식 (9)를 이산화하면 추적필터 설계를 위한 탄도비행 운동모델을 얻는다.
위의 식에서 사용된 행렬의 정의는 다음과 같다.
식 (10)에서 운동모델의 불확실성을 고려하기 위한 백색 공정잡음 $ w_{b}\sim N(0,\: Q_{b})$은 분산이 $Q_{b}$인 영평균 정규분포를
따르는 것으로 가정한다.
• 상승기동 운동모델 (SMM: skip motion model)
3축 기동 파라미터를 모두 고려하는 경우, 식 (3)으로부터 MaRV의 운동모델을 9차 비선형 미분방정식으로 쓸 수 있다.
여기서
앞서와 마찬가지로 식 (11)을 샘플링주기 $T$로 이산화하면, 상승기동 운동모델을 유도할 수 있다.
위의 식에 사용된 행렬 및 변수의 정의는 다음과 같다.
• 측정모델
궤적 예측을 위해 레이더에서 제공된 표적 항적을 활용할 수 있다면, MaRV의 비행구간 별 측정방정식을 다음과 같이 모델링할 수 있다.
여기서 측정행렬 $H_{b}$와 $H_{s}$는 다음과 같이 정의된다.
또한 식 (13)에서 측정잡음 $ v\sim N(0,\: R)$은 분산이 $R$인 영평균 정규분포를 따르는 것으로 간주한다. 통상 레이더에서 표적 항적과 더불어 분산이
제공되므로, 이를 이용해 $R$을 설정할 수 있다.
• 비행구간 별 MaRV 상태추정기
앞서 유도된 운동모델 (10)과 (12), 그리고 측정모델 (13)을 결합하면 이산시간 비선형 상태공간 방정식을 얻을 수 있다. 따라서 표 1에 정리된 UKF(Unscented Kalman Filter)를 적용하면 MaRV의 상태변수를 추정할 수 있다[19]. 지금부터는 탄도비행 및 상승기동 구간 운동모델을 이용해 설계된 상태추정기를 각각 BMM 부필터와 SMM 부필터로 명명한다.
표 1 UKF 순환식
Table 1 Unscented Kalman filter recursion
|
system model
|
$\begin{cases}
x_{k+1}&= f(x_{k})+ w_{k}\\
z_{k}&= h(x_{k})+ v_{k}
\end{cases},\: w_{k}\sim N(0,\: Q_{k}),\: v_{k}\sim N(0,\: R_{k})$
|
|
filter
recursion
|
• system propagation
$\begin{align*}
\hat{x}_{k \vert k-1}
& =\sum_{i=0}^{2L}w_{i}^{m}\chi_{i,\: k}^{x}\\
P_{k \vert k-1}& =\sum_{i=0}^{2L}w_{i}^{c}(\chi_{i,\: k}^{x}-\hat{x}_{k
\vert k-1})(\chi_{i,\: k}^{x}-\hat{x}_{k \vert k-1})^{T}+ Q_{k}
\end{align*}$
- $i=0$
$\begin{align*}
\chi_{0,\: k}^{x}
& = f(\hat{x}_{k-1 \vert k-1}),\: \\
w_{0}^{m}& =\kappa /(L+\kappa),\: w_{0}^{c}= w_{0}^{m}+ 1 -\alpha^{2}+\beta
\end{align*}$
- $i neq 0$
$ \begin{align*}
\chi_{i,\: k}^{x}
& = f(\hat{x}_{k-1 \vert k-1})+\left(\sqrt{(L+\kappa)P_{k-1 \vert k-1}}\right)_{i},\:
\\
\chi_{i,\: k}^{x}& = f(\hat{x}_{k-1 \vert k-1})-\left(\sqrt{(L+\kappa)P_{k-1
\vert k-1}}\right)_{i-L},\: \\
w_{i}^{m}& = w_{i}^{c}= 1/[2(L+\kappa)]
\end{align*}$
• measurement update
$\begin{align*}
K_{k}
& = P_{k \vert k-1}H_{k}^{T}(H_{k}P_{k \vert k-1}H_{k}^{T}+ R_{k})^{-1}\\
\hat{x}_{k \vert k}& =\hat{x}_{k \vert k-1}+ K_{k}(z_{k}-H_{k}\hat{x}_{k
\vert k-1})\\
P_{k \vert k}& =(I - K_{k}H_{k})P_{k \vert k-1}
\end{align*}$
|
대부분의 MaRV는 상승기동을 반복적으로 수행하지 않으므로 그림 4와 같이 운동모델 간 천이 가능성이 배제된 필터 뱅크를 구성할 수 있다. 이 구조 하에서 BMM 부필터는 상승기동을 탐지하는 역할을, SMM 부필터는
상승기동 탐지 후 궤적 예측 알고리듬의 입력변수인 양항비를 추정하는 역할을 수행한다.
그림 4. 기동탐지 및 상태추정을 위한 필터 뱅크
Fig. 4. Filter bank for maneuver detection and state estimation
• 상승기동 탐지
BMM 부필터의 잔차와 잔차 공분산은 다음과 같이 계산된다.
식 (14) 및 식 (15)을 이용해 정규화된 잔차 제곱합(NIS: normalized innovation squared)을 정의할 수 있다.
잘 알려져 있듯이 이론적으로 BMM 부필터 설계에 사용된 모델에 불확실성이 없는 경우, 다시 말해 MaRV가 상승기동을 수행하지 않으면 잔차 $\nu_{k}$는
평균이 0이고 분산이 $S_{k}$인 백색 정규잡음이다. 이 경우, 식 (15)에서 정의한 정규화된 잔차 제곱합 $ϵ_{k}$는 자유도가 dim($z_{k}$)$=3$인 카이제곱 분포를 따른다. 반면, MaRV가 상승기동을 수행하는
경우에는 운동모델 불확실성으로 인해 BMM 부필터의 잔차가 편향되어(biased) $\epsilon_{k}$의 평균적 특성이 크게 달라진다.
이러한 사실에 기초하여 MaRV이 상승기동 한다는 대립가설(alternative hypothesis)을 $H_{1}$, 귀무가설(null hypothesis)을
$H_{0}$라 하면, 다음과 같이 BMM 부필터 잔차에 대한 카이제곱 검정을 통해 MaRV의 상승기동 여부를 손쉽게 판별할 수 있다.
여기서 $\tau$는 상승기동 탐지를 위한 임계값을 의미한다. 카이제곱 분포의 자유도가 3이므로, 허용 가능한 오경보율(false alarm rate)이
0.1%일 때 임계값은 $\tau =16.27$로 설정된다.
3. 재진입 비행체의 상승기동 궤적 예측
본 절에서는 상승기동 궤적 해석 결과에 기반한 궤적 예측 알고리듬을 설계한다. 상승기동 구간에서의 MaRV 궤적 예측 프로세스는 그림 5에 도시한 바와 같다. 제안한 알고리듬은 2절의 필터뱅크로부터 매 시점 산출된 MaRV 위치/속도 및 기동 파라미터 추정치를 입력받아 대기권 재이탈
시점 $Q_{o}$에서의 MaRV의 위치/속도 벡터를 계산하는 Step 1, 이 결과를 초기조건으로 하여 수치적분을 통해 상승기동 궤적을 전파하는
Step 2로 구성되어 상승기동 전 구간에 대한 궤적 예측 결과를 도출한다.
그림 5. 상승궤적 예측 프로세스
Fig. 5. Process of skip trajectory prediction
• Step 1. 대기권 이탈 위치/속도 산출 ($Q_{i}\sim Q_{o}$ 구간)
상승기동에 의해 MaRV가 대기권을 재이탈하는 위치와 속도는 상승기동 궤적의 닫힌 해를 통해 산출된다. 수학적 분석에 사용되는 기본 가정은 다음과
같다.
A1. 대기권 재진입 후 MaRV는 공력학적 효율을 감안하여 비행평면 상에서 특정 양항비로 상승기동 한다.
A2. 대기밀도는 고도에 따라 지수적으로 감소한다 (식 (7)).
A3. 상승기동 구간에서 지구는 반지름이 $R_{e}$인 구형 지구로 근사할 수 있으며, 지구자전 효과도 무시할 만 하다.
$Q_{i}\sim Q_{o}$ 구간에 대한 MaRV의 궤적은 그림 6에 도시한 바와 같다. 케플러 제 1 법칙과 제 2 법칙을 적용하면 극좌표계 MaRV 위치 $(r,\:\theta)$를 얻을 수 있다[20].
그림 6. $Q_{i}\sim Q_{o}$구간 비행궤적
Fig. 6. Skip trajectory from $Q_{i}$ to $Q_{o}$
이때, 단위질량당 각 운동량(angular momentum) $\lambda$와 상승기동 궤적의 이심률 $e$은 다음과 같이 정의된다.
위 식에서 $\gamma$는 비행경로각(flight path angle), 아랫첨자 $i$는 $Q_{i}$ 지점에서의 물리량을 의미한다. 식 (18)을 식 (17)에 대입하여 정리하면 $Q_{i}\sim Q_{o}$ 각거리(angular distance) $\Phi_{s}$를 계산할 수 있다.
위의 식에서 $v_{s}=\sqrt{GM/R_{e}}$는 지표면에서의 위성속도(satellite speed) 이다.
상승기동 시 MaRV에 작용하는 공력학적 힘은 중력에 비해 상대적으로 매우 크므로, 그림 6으로부터 수직면 운동방정식을 다음과 같이 근사적으로 쓸 수 있다.
여기서 $r_{c}$는 MaRV 비행궤적의 곡률반경이다. MaRV의 고도를 $h$라고 하면 비행경로 $s$의 증분 $ds$는 다음 식으로 근사된다.
식 (6) 식 (7), 식 (22)를 식 (21)에 대입하여 정리하면 다음식을 얻는다.
위 식의 양변을 재진입 지점 $Q_{i}$에서 대기권 이탈 지점 $Q_{o}$까지 적분하면 식 (24)을 얻는다.
대기권계면, 즉 $Q_{i}$ 및 $Q_{o}$에서는 공기가 매우 희박하므로 위 식의 좌변이 0이된다. 즉,
이상의 결과로부터 MaRV가 상승기동을 수행하는 경우 $Q_{o}$과 $Q_{i}$에서의 비행경로각은 크기는 동일하고 부호가 반대가 됨을 알 수 있다.
한편, 식 (20) ~ (22)를 통해 MaRV의 양항비 $\eta$를 계산하면 다음과 같다.
연쇄법칙(chain rule)을 적용하면 MaRV의 속력 변화율 $\dot{v}$을 다시 쓸 수 있다.
식 (26)을 식 (27)에 대입하면 다음 식을 얻는다.
따라서, $Q_{i}\sim Q_{o}$ 구간에 대해 식 (28)의 양변을 적분하면 $Q_{o}$에서의 속력 $v_{o}= | V_{o}^{E}|$을 대기권 진입속력 $v_{i}= | V_{i}^{E}|$와 양항비
$\eta$의 함수로 쓸 수 있다.
이제 앞서 유도된 각거리 $\Phi_{s}$와 속력 $v_{o}= | V_{o}^{E}|$을 이용해 대기권 재이탈 지점 $Q_{o}$에서의 위치벡터
$ R_{o}^{E}$와 속도벡터 $ V_{o}^{E}$를 계산해 보자. 가정 A1에 의해 $ R_{o}^{E}$와 $ V_{o}^{E}$는 대기권
재진입 지점 $Q_{i}$에서의 위치벡터 $ R_{i}^{E}$ 및 속도벡터$ V_{i}^{E}$를 포함하는 평면 상에 놓여야 한다. 가정 A3에
의해 $| R_{o}^{E}| = | R_{i}^{E}|=r_{i}$이므로, $ R_{o}^{E}$와 $ V_{o}^{E}$를 다음과 같이 정의된다.
여기서 $(c_{1},\: c_{2})$, $(d_{1},\: d_{2})$는 $\left | c_{1}R_{i}^{E}+ c_{2}V_{i}^{E}\right
| =$$\left | d_{1}R_{i}^{E}+ d_{2}V_{i}^{E}\right | = 1$을 만족하는 임의의 상수이다.
그림 6으로부터 $ R_{i}^{E}$와 $ R_{o}^{E}$ 사이의 기하학적 관계를 다음과 같이 쓸 수 있다.
식 (30)을 식 (31)에 대입하면 위치벡터 $ R_{o}^{E}$에 사용된 계수를 다음과 같이 결정할 수 있다.
속도벡터 $ V_{o}^{E}$에 사용된 미지계수 $(d_{1},\: d_{2})$ 역시 앞서와 동일한 방법으로 계산할 수 있다.
• Step 2. 대기권 이탈 후 궤적 전파 ($Q_{o}\sim Q_{p}$ 구간)
대기권 이탈 후에는 공력학적 힘이 거의 작용하지 않으므로, MaRV의 운동은 탄도비행 모델로 근사할 수 있다. 식 (3)에서 동압을 0으로 설정하면 대기권 이탈 후 MaRV의 비행궤적 산출에 필요한 운동방정식을 쓸 수 있다.
위 미분방정식의 초기치를 식 (31)으로 설정한 후 수치적분하면 대기권 재이탈 후의 궤적을 예측할 수 있다. 상승기동 궤적의 정점은 수직면 비행경로각이 0이 될 때이므로, 해당 조건이
충족할 때까지 궤적 전파를 수행하여 상승기동 구간 예측 결과를 최종 출력한다.
Remark 3. 앞서 설명한 바와 같이 제안한 궤적 예측 기법은 추정된 양항비를 활용하여 $Q_{i}\sim Q_{o}$ 구간 궤적을 효과적으로
재구성할 수 있다. 비록 대기권 이탈 후, 즉 $Q_{o}\sim Q_{p}$ 구간에서의 궤적 예측을 위해 운동모델 (9)를 활용한 수치적분이 불가피하지만, 해당 구간에서 공기밀도가 무시할만한 수준이므로 모델 불확실성이 예측 결과에 미치는 영향은 미미한 수준에 그칠 것임을
예상해 볼 수 있다. 즉, 제안한 기법은 MaRV의 공력 모델 불확실성에 대한 궤적 예측 알고리듬의 견실성 확보에 큰 도움이 된다.
4. 성능분석
제안한 궤적 예측 알고리듬의 유효성을 평가하기 위해, 그림 7~9에 도시된 바와 같이 비행사거리가 각각 300km, 400km 인 두가지 전형적인 MaRV 궤적에 대하여 모의실험을 수행하였다. 그림 7~8에서 검은색 실선은 MaRV의 발사부터 탄착까지의 전체 궤적을, 녹색은 상승궤적 구간, 적색 점은 상승궤적 정점을 도시한 것이다. 그림 8의 세로축의 값은 상승기동 개시 고도 $h_{i}$에 대한 상대 고도이다. 그림 9의 세로 점선은 상승궤적 구간 시작 및 종료 시점을 나타낸다.
상승기동 구간 궤적 예측 성능을 비교하기 위해 제안한 기법과 함께, 필터에서 산출된 기동 파라미터가 계속 유지된다는 가정 하에 운동방정식 (11)을 수치적분하는 기존 방법을 함께 시뮬레이션 하였다. 궤적 예측성능은 100회 반복시행을 통해 평가되었으며, 모의실험 환경은 표 2에 정리된 바와 같다. 궤적 예측 알고리듬의 안정적인 동작을 위해, MaRV의 기동 탐지를 위한 BMM 부필터의 공정잡음 $Q_{b}$와 MaRV
상태 및 기동파라미터 추정을 위한 SMM 부필터의 공정잡음 $Q_{s}$을 충분히 작게 설정하였다. 레이더는 자체 시스템에서 E-frame 추적정보
및 탄종 식별정보를 0.2초 간격으로 제공하며, 제공되는 추적정보의 위치 오차는 15m 수준으로 설정하였다.
그림 7. 지구좌표계에서의 3차원 궤적
Fig. 7. 3-dimensional target trajectory in E-frame
그림 8. 수직면 궤적
Fig. 8. Vertical trajectory
그림 9. 속력 프로파일
Fig. 9. Speed profile
표 2 모의실험 조건
Table 2 Simulation condition
|
item
|
specification
|
|
noise
statistics
|
$Q_{b}=diag(0.01,\: 0.01,\: 0.01)$
$Q_{s}=diag(4.47,\: 7.28,\: 7.28,\: 0.0,\: 0.0,\: 0.0)\times 10^{-6}$
$R = diag(15.0,\: 15.0,\: 15.0)[m]$
|
|
sampling
period
|
$T=0.2[\sec]$
|
4.1 상승기동 탐지 및 양항비 추정성능
그림 10은 2절에서 설계된 필터 뱅크의 NIS(normalized innovation squared)를 도시한 결과이다. 기동 탐지 알고리듬은 각 부필터의
NIS가 임계치 $\tau =16.27$를 초과하는 경우 상승기동이 이루어진 것으로 판단한다. 그림에서 확인할 수 있듯이 BMM 부필터는 기동 개시
약 3초 이후에 NIS가 임계치를 초과하며, 이 시점에 검출 로직에 따라 기동이 탐지된다. SMM 부필터는 기동 개시 직후에는 BMM 부필터와 마찬가지로
NIS가 증가하여 임계치를 초과하지만, 기동 개시 이후 약 10초 후부터는 필터가 수렴하면서 NIS가 감소하는 경향을 보인다. 이 시점부터 SMM
부필터의 위치/속도 오차와 더불어 양항비 추정오차 또한 수렴하는 것을 그림 11을 통해 확인할 수 있다. 위치 추정치의 수렴속도는 약 15초, 속도 추정치의 수렴속도는 약 13초, 양항비 추정치의 수렴속도는 약 9초로 MaRV의
운동에 직접적으로 관여하는 양항비 추정치가 상대적으로 빠르게 수렴함을 알 수 있다. 더욱이 상태추정 필터가 수렴 상태에 도달한 이후에는 양항비 추정결과가
상당한 정확도 수준을 보이므로, 이를 상승궤적 예측에 활용하는데에 무리가 없다.
그림 10. 상승기동 전후 필터 뱅크 정규화 잔차
Fig. 10. Normalized residual of filter bank
그림 11. 표적 상태추정치 RMSE
Fig. 11. RMSE of target state estimates
4.2 상승기동 궤적 예측 성능
기동 탐지 이후 SMM 부필터가 수렴하게 되면, 궤적 예측을 수행할 수 있다. SMM 부필터의 출력을 사용하여 궤적예측기의 초기치를 설정하며, 마찬가지로
SMM 부필터에서 기동필터의 기동 탐지 필터 수렴 후 양항비 추정치를 산출하여 상승기동시의 재이탈 위치 및 속도를 산출하고 궤적 전파를 통해 상승궤적
정점을 예측할 수 있다.
궤적 예측 성능분석의 편의를 도모하기 위해, MaRV의 상승기동 궤적의 정점을 주요 특징점으로 설정하고, 이를 기준으로 수직면에서의 CEP(circular
error probable)를 분석한다. CEP는 특징점 예측오차의 50%를 포함하는 원의 반지름을 나타내는 것이다[21]. 따라서, CEP가 작다는 것은 궤적 예측 알고리듬의 정확도가 우수하다는 것을 의미한다.
• 상승기동 궤적 예측 오차
시뮬레이션 결과 각 알고리듬의 MaRV 상승기동 궤적 예측 성능은 그림 12~13와 같다. 그림 12는 기동 탐지 후 상태추정 필터가 수렴한 시점, 즉 상승기동 후 약 10초가 경과한 시점(그림 11 참조)에 상승궤적 정점을 예측한 결과를 나타낸다. 시뮬레이션 결과로부터 기존 방법 대비, 제안 기법의 CEP의 반경과 편향오차 모두 대폭 감소 함을
알 수 있다. 제안한 기법의 특징점 예측오차 평균과 CEP 반경을 더한 값은 약 0.45~0.47km 수준인데 반해, 기동파라미터 추정치를 이용해
운동모델을 단순 수치 적분하는 기존 방법은 2.97~3.06km 수준이다. 이는 약 4~7배 가량의 성능 향상 효과에 해당한다.
그림 13은 매 시점 획득된 양항비 추정치를 이용해 상승기동 궤적을 예측하고, 특징점을 기준으로 한 RMSE를 도시한 것이다. 상태추정 필터가 완전히 수렴하여
신뢰성 있는 양항비 추정결과가 획득되더라도 기존 기법의 예측오차는 최대 약 2.20km까지 증가하다가 상승기동 후 약 50초가 경과해야 만족할 만한
예측결과를 보이는 반면, 제안한 기법은 예측오차의 최대값이 0.50km 이내로 유지될 뿐만 아니라 예측오차의 수렴속도도 매우 빠름을 알 수 있다.
• 모델 불확실성의 영향
운동모델의 부정확성은 궤적 전파 과정에 직접적인 영향을 끼치므로 이를 통해 궤적 예측방법의 모델 의존성을 확인할 수 있다. 이를 위해, 의도적인 파라미터
불확실성을 다음과 같이 부가한 후 시뮬레이션을 수행하였다.
- 경우 1 : 공력 불확실성이 없는 경우
- 경우 2 : 공력 불확실성 –10%
- 경우 3 : 공력 불확실성 +10%
모델 불확실성이 특징점 예측 성능에 미치는 영향은 표 3에 정리한 바와 같다. 표에서 확인할 수 있듯이 기존 방법에 의해 예측된 특징점 CEP는 모델 불확실성이 없는 경우 1에서 약 3.0km 수준이었다가
모델 불확실성이 추가되면 최대 약 4.8km까지 크게 증가한다. 기존 기법의 궤적 예측 성능 저하는 기동탐지 시점에서 운동방정식 (11)을 수치적분 하는 과정에서 운동모델 불확실성이 누적된 결과로 볼 수 있다. 이와 달리, 제안한 기법은 공력 불확실성이 존재하더라도 예측 오차의 변화가
약 0.6km를 넘지 않아 기존 기법에 비해 견실한 궤적 예측성능을 보여준다. 이는 Remark 3에서 언급한 바와 같이 제안한 기법은 대기권 안쪽에서
상승궤적의 닫힌 해를 활용하여 궤적 예측을 수행하므로 공력 불확실성에 상대적으로 둔감한 특성을 지니기 때문이다.
이상의 결과를 종합할 때 제안 기법이 MaRV의 미지 기동과 모델 불확실성이 존재하는 상황에서도 장시간 상승궤적 예측에 매우 효과적인 방법이라는 결론을
내릴 수 있다.
그림 12. 상승기동 궤적 예측 알고리듬의 성능비교
Fig. 12. Performance comparison of the skip trajectory prediction algorithms
그림 13. 시간에 따른 상승궤적 예측 CEP
Fig. 13. Time history prediction CEP of skip trajectory
표 3 모델 불확실성에 따른 CEP 변화
Table 3 CEP variations along model uncertainties
|
Scenario
|
previous method [km]
|
proposed method [km]
|
|
case 1
|
case 2
|
case 3
|
case 1
|
case 2
|
case 3
|
|
#1
|
2.966
|
4.050
|
1.461
|
0.471
|
0.867
|
0.240
|
|
#2
|
3.056
|
4.834
|
1.926
|
0.448
|
1.051
|
0.766
|