2.1 CIWS 센서 구성 및 측정치 특성

CIWS와 표적 간 상대기하는 그림 1에 도시된 바와 같다. 그림에서 관성좌표계 $(X_{I}, Y_{I}, Z_{I})$는 원점이 CIWS인 NED 좌표계이다. $\vec{R}^{I}\equiv [\begin{matrix}x&y&z\end{matrix}]^{T}$는 관성좌표계 표적 위치벡터, $r$은 상대거리, $\psi$는 방위각, $\theta$는 고각을 의미한다. 표적 시선방향을 $X_{L}$축으로 하는 오른손좌표계를 시선좌표계라 정의하면, 관성좌표계에서 시선좌표계로의 좌표변환행렬은 $C_{I}^{L}=R_{y}(\theta)R_{z}(\psi)$로 표현된다.

그림 1. 상대기하

Fig. 1. Relative geometry

일반적으로 CIWS에는 상호보완적 특성을 지닌 레이더와 전자광학센서가 장착되어 있다. 능동센서인 레이더는 탐지거리가 길고 기상변화에 둔감하지만 다중경로의 영향을 많이 받는 반면, 전자광학센서는 거리정보 획득이 불가능한 피동센서이지만 다중경로의 영향을 거의 받지 않는다.

대부분의 CIWS 추적레이더는 MPRF(medium pulse repetition frequency) 신호처리 방식을 채택하고 있으므로, 방위각, 고각 뿐만 아니라 상대거리 및 상대거리 변화율 측정치를 제공하는 것으로 가정할 수 있다.

여기서 $r$은 상대거리, $\dot{r}$은 상대거리 변화율, $\psi$는 방위각, $\theta$는 고각 참값을 의미한다. 만일 다중경로 오차가 없다면, 측정잡음 $\delta r$, $\delta\dot{r}$, $\delta\psi$, $\delta\theta$는 영평균 정규잡음으로 가정할 수 있다.

위의 식에서 $\rho$는 상관계수를 의미한다.

불행하게도 해상환경에서 운용되는 CIWS 레이더의 표적 측정치는 해면반사의 영향을 받아 통계적 특성이 크게 달라진다 ^[18]. 해면반사로 인해 간접경로로 수신되는 신호는 레이더-표적 간의 직접경로로 수신되는 신호와 상쇄/보강 간섭을 일으켜, 그림 2에 도시된 바와 같이 SNR(Signal-to-Noise Ratio)의 증감을 초래한다. 따라서, 다중경로 효과는 표적 탐지 성능에도 직접적인 영향을 줄 뿐만 아니라, 레이더 수신신호의 위상을 왜곡시켜 결과적으로 모노펄스 연산 과정에서 비정상적인 고각 측정 오차를 야기한다. 그림 3은 다중경로에 의한 방위각 및 고각 측정오차의 평균(편향오차)과 표준편차의 변화 양상을 보여준다. 그림 3에서 좌측 y축은 방위각, 우측 y축은 고각에 대한 측정오차 수준을 의미한다. 고각 측정오차 표준차의 최대값은 0.34 deg로 방위각(0.013 deg) 대비 약 25배 이상 큰 값을 보인다. 이는 참고문헌 ^[20]과 대동소이한 결과로, 다중경로가 주로 고각 오차에 영향을 미친다는 것을 시사한다.

그림 2. 다중경로 유무에 따른 신호 대 잡음비 비교

Fig. 2. The impact of multi-path on signal-to-noise ratio

그림 3. 다중경로 상황에서의 각도 측정오차 평균 및 표준편차

Fig. 3. Error statistics of angle measurement errors under multipath environments

Remark 1. CIWS 레이더와 같이 빔폭이 상당히 작은(pencil beam) 경우, 수평면에서는 해수면 반사가 거의 대칭적으로 발생하기 때문에 다중경로 효과가 방위각 측정치에 미치는 영향은 미미하다 ^[19-21]. 또한 간접경로에 의한 거리 차는 레이더 거리 분해능에 비해 무시할 만하므로, 일단 표적이 탐지되면 상대거리 측정치도 정상적으로 산출된다. 즉, 다중경로 상황에서 고각을 제외한 나머지 레이더 측정치들은 식 (1)과 동일한 통계적 특성을 갖는 것으로 근사해도 무방하다.

레이더가 다중경로 효과로 인해 비정상적으로 큰 오차를 포함한 고각 측정치를 사용하여 추적필터를 구동하면, 표적 상태추정치가 더 이상 영평균 특성을 유지하지 못하고 크게 편향될 수밖에 없다. 따라서 CIWS는 해상환경 변화 시에도 안정적으로 대함유도탄 표적을 추적하기 위해 다중경로에 상대적으로 강인한 전자광학센서를 보조수단으로 활용한다.

전자광학센서는 표적 방위각, 고각의 측정 정확도가 높아 정보융합 시 다중경로로 인한 레이더 추적성능 열화를 효과적으로 억제할 수 있다. 그러나 영상처리의 높은 연산량으로 인해 레이더 대비 측정주기가 길다. 즉, 통신채널을 통해 융합센터로 전송된 전자광학센서 측정치는 현재 시점 $t_{k}$가 아니라, 훨씬 이전 시점 $t_{d}$에 산출되었을 가능성이 농후하다. 또한, 통상 레이더와 전자광학센서는 정보 제공주기가 상이하므로, 융합필터 입장에서는 전자광학센서로부터 시간지연(latency)을 포함한 비동기 측정치(OOSM)가 획득되는 것으로 간주할 수 있다.

여기서 $(\widetilde{\beta}_{d}, \widetilde{\epsilon}_{d})$는 현재 시점 $t_{k}$에서 전자광학센서로부터 전송받은 방위각 및 고각 측정치이다. 편의상, 측정잡음 $(\delta\beta_{d}, \delta\epsilon_{d})$는 영평균 정규분포를 따르는 것으로 가정한다.

식 (2)에서 전자광학센서가 측정치를 실제로 산출한 시점 $t_{d}$는 현재 시점 $t_{k}$를 기준으로 $T_{d}$ 만큼 지연된 시점이다. 즉, $t_{d}=t_{k}-T_{d}$. 통상 레이더 측정주기 $T_{r}= t_{k}- t_{k-1}$에 비해, 전자광학센서 측정치의 시간지연이 매우 크므로 $(T_{d}\gg T_{r})$, 양의 정수 $l\gg 1$에 대해 다음 조건이 만족된다.

2.2 비동기 이종센서 융합 문제정의

CIWS는 이종센서인 레이더와 전자광학센서 정보를 융합하기 위해 그림 4의 트랙-측정치 하이브리드 정보융합 구조를 사용한다. 주 센서인 레이더는 자체 표적추적 알고리듬을 내장하고 있어 레이더 표적트랙, 다시 말해 표적 상태 추정치 $\mathbf{z}_{k}$ 및 추정오차 공분산 $\Xi_{k}$을 융합필터에 제공한다. 이와 달리, 전자광학센서는 식 (2)에 정의된 방위각, 고각 측정치 $\mathbf{y}_{d}$를 융합센터로 직접 송신한다. 융합필터는 통상 주기가 빠른 레이더와 동기화되어 구현되며, 레이더 자체트랙과 전자광학센서 측정치를 이용해 개선된 표적트랙 $(\hat{\mathbf{x}}_{k}, P_{k})$를 생성한다.

그림 4. CIWS의 정보융합 구조

Fig. 4. Data fusion architecture of the CIWS

전술한 바와 같이 CIWS 정보융합 필터 설계 목표는 레이더 자체트랙과 전자광학센서 측정치를 이용하여 다음 운동방정식을 만족하는 표적 상태변수 $\mathbf{x}_{k}$를 추정하는 것이다.

여기서 $\mathbf{x}_{k}$는 시간 $t_{k}$에서의 3차원 표적 위치, 속도, 가속도, $\mathbf{w}_{k, k-1}$는 표적 운동모델의 불확실성을 반영하기 위한 공정잡음으로 영평균 정규분포를 따른다. 시스템 행렬 $\Phi_{k, k-1}$와 공정잡음 분산 $Q_{k, k-1}$은 표적 운동을 어떻게 모델링하느냐에 따라 달라지는데, 본 논문에서는 편의상 참고문헌 ^[22]의 등가속 운동모델을 사용하였다. 참고로, 샘플링주기가 충분히 짧은 경우 등가속도 운동모델로도 충분히 추적필터 설계를 위한 표적 운동을 모사할 수 있다.

한편, 융합필터 설계를 위한 측정방정식은 $t_{k}$ 시점에서의 전자광학센서 측정치 (2)의 수신 여부에 따라 달라진다.

위의 식에서 레이더는 자체트랙이 융합센터로 제공하므로 행렬 $H_{k}=I$이고, 식 (2)로부터 측정함수 $h(\bullet)$는 $t_{d}$시점에서의 표적 상대위치 $(x_{d}, y_{d}, z_{d})$의 비선형 함수로 기술된다.

이제, 비동기 이종센서 정보융합 필터 설계 문제를 명확히 정의하기 위해 융합필터 내에서 이루어지는 정보처리 과정을 살펴보자. 그림 5에서 레이더는 매 시점 자체 표적트랙 $\mathbf{z}_{k}$(흑색 세모)를 제공하며, 전자광학센서 측정치 $\mathbf{y}_{d}$(적색 세모)는 매우 드물게 융합필터에 도달한다. 실제로 $t_{k}$ 시점에 획득된 전자광학센서의 측정치는 $t_{d}$시점에 획득된 것이므로, 융합필터는 그 사이, 즉 $t_{k-l}\sim t_{k}$ 시점에 주어진 측정치로 융합트랙을 생성할 수 밖에 없으며, 시간지연 측정치 $\mathbf{y}_{d}$가 도달하면 그제서야 이를 처리하게 된다. 따라서 이종센서 융합은 사실상 $t_{k-l}$ 시점에서의 확률밀도함수를 시구간 $t_{k-l}\sim t_{k}$에서 획득된 레이더 정보 $Z^{k}=\left\{\mathbf{z}_{k-l+1}, \mathbf{z}_{k-l+2}, \cdots , \mathbf{z}_{k}\right\}$와 $t_{d}$시점에서 생성되었으나 뒤늦게 $t_{k}$시점에서야 주어진 전자광학센서 측정치 $y_{d}$를 이용해 현재시점 $t_{k}$에서의 확률밀도함수 $p(\mathbf{x}_{k}|Z^{k}, \mathbf{y}_{d})$를 산출하는 OOSM 문제의 일종이라 할 수 있다.

그림 5. 융합필터 시간 다이어그램

Fig. 5. Fusion filter time diagram

Remark 2. CIWS의 이종센서 정보융합 문제는 가용 센서 정보의 불편향성을 전제하는 일반적인 비동기 센서 정보융합, 즉 OOSM 문제와는 본질적으로 성격이 다르다. 다중경로 효과가 두드러지는 구간에서는 CIWS 레이더 고각 측정치에 비정상적인 편향오차가 포함되므로, 이를 이용해 산출된 레이더의 자체 표적트랙은 더 이상 불편향성(unbiasedness)을 유지하지 못한다. 따라서 전자광학센서에서 불편향 고각 측정치 (2)가 간헐적으로 제공되더라도 역전파(retrodiction) 개념을 활용하는 기존 OOSM 필터로는 수학적으로 최적 성능을 담보할 수 없다.

3.1. 역전파 기반 융합기법의 한계

이제 그림 5의 사후확률밀도 함수 갱신 과정을 수식으로 기술해보자. $k$시점까지 융합센터에 제공된 측정치는 시간 순으로 $\left\{\mathbf{z}_{1}, \cdots , \mathbf{z}_{k-l}, \mathbf{y}_{d}, \mathbf{z}_{k-l+1}, \cdots , \mathbf{z}_{k}\right\}$와 같이 재정렬할 수 있으며, 상태변수 $\mathbf{x}$에 대한 추정치는 시간 정렬된 측정치로 확률공간을 순차 갱신함으로써 산출 가능하다.

먼저 레이더 측정치 집합 $Z^{k-l}\equiv \left\{z_{1}, \cdots , z_{k-l}\right\}$을 통해 갱신된 상태추정치와 추정오차 공분산 행렬을 다음과 같이 정의하자.

시간지연 측정치 $\mathbf{y}_{d}$는 상태변수 $\mathbf{x}_{d}$의 함수로 기술되므로, $\mathbf{y}_{d}$를 이용하여 확률공간을 갱신하기 위해서는 $t_{k-l}$시점 상태추정치를 기반으로 $t_{d}$시점의 상태를 예측하는 과정이 필요하다. 이러한 시스템 전파 및 확률공간의 갱신 과정은 채프먼-콜모고로프(Chapman-Kolmogorov) 방정식과 베이즈 정리(Bayes’ theorem)에 따라 아래와 같이 기술된다.

여기서 $c_{d}$는 정규화 상수이다.

마찬가지로, 이후에 레이더 측정치 $\mathbf{z}_{k-l+1}$가 주어지면 확률공간 (10)은 다음과 같이 갱신된다.

$k$시점까지 레이더 측정치가 매 시점 주어지면, 위 과정을 반복하여 확률공간을 갱신할 수 있으며 그 결과는 아래와 같다.

참고로, 만일 최종적으로 산출된 확률밀도함수 (14)가 정규분포 혹은 단봉(unimodal) 특성을 갖는다면, $t_{k}$시점 상태추정치 $\hat{\mathbf{x}}_{k|k}\equiv E\left\{p(\mathbf{x}_{k}| Z^{k}, \mathbf{y}_{d})\right\}$ 및 추정오차공분산 $P_{k|k}\equiv var\left\{p(\mathbf{x}_{k}| Z^{k}, \mathbf{y}_{d})\right\}$을 손쉽게 얻을 수 있다.

전술한 확률공간 갱신 결과는 식 (4)의 시스템 모델을 이용해 역전파를 활용한 OOSM 필터(표 1)로도 구현 가능하다 ^[17]. OOSM 필터는 식 (7), (8)로 기술되는 확률공간을 측정치 획득 시점 순으로 순차 갱신하는 것이 아니라, 일단 식 (9), (10)을 건너뛴 후 $t_{k}$시점에 이르러서야 역전파를 거쳐 $t_{d}$ 시점 정보 $\mathbf{y}_{d}$를 이용해 확률공간을 최종 갱신하는 방식을 취한다.

전술한 바와 같이 OOSM 필터는 측정치 획득시점에 따라 확률공간을 순차 갱신해 가는 표준 칼만필터와 구현 방법에 차이가 있을 뿐 사실상 동일한 역할을 수행한다. 즉, OOSM 필터는 표준 칼만필터 문제와 동일한 가정이 충족될 때만 비동기 이종센서 융합문제의 수학적 최적 해가 된다. 잘 알려져 있듯이, 표준 칼만필터는 세 가지 가정, 즉 1) 주어진 시스템 모델에 파라미터 불확실성이 존재하지 않고, 2) 공정잡음과 측정잡음 간 상관성이 존재하지 않으며, 3) 가용 측정치가 영평균 정규분포를 따르는 경우 상태추정치의 불편향성을 만족하는 최소분산 추정기로 동작한다 ^[23].

이러한 관점에서 볼 때, 기존 OOSM 필터는 레이더가 다중경로의 영향을 빈번하게 받는 CIWS에 응용하기에는 성능측면에서 그 한계가 명확하다. 이는 식 (13)에서 확인 가능하다. 식 (13)에서 $p(\mathbf{x}_{i}|\mathbf{x}_{i-1})$는 $t_{i-1}$에서 $t_{i}$시점으로 상태를 예측하는 과정을 나타내고, $p(\mathbf{z}_{i}| \mathbf{x}_{i})$는 측정치 $\mathbf{z}_{i}$를 이용하여 $t_{i}$시점 상태예측치를 보정하는 것을 의미한다. $p(\mathbf{x}_{k-l +1}| Z^{k-l +1}, \mathbf{y}_{d})$는 식 (11)을 통해 계산되는 값으로, 전자광학센서 측정치 $\mathbf{y}_{d}$로 보정된 확률공간 $p(\mathbf{x}_{d}| Z^{k-l}, \mathbf{y}_{d})$을 레이더 측정치 $\mathbf{z}_{k-l+1}$로 한 번 더 갱신한 결과이다. 즉, 식 (13)으로부터 전자광학센서 측정치 $\mathbf{y}_{d}$를 이용해 확률공간을 갱신하더라도, 측정치 산출시점 $t_{d}$이후 여러 시점에 걸쳐 획득된 레이더 정보 $\left\{\mathbf{z}_{k-l+1}, \cdots , \mathbf{z}_{k}\right\}$가 다중경로의 영향을 받아 편향되어 있다면 만족할 만한 표적추적 성능을 얻기란 사실상 불가능하다.

Remark 3. 표 1의 OOSM 필터는 보조 센서 측정치가 실제로 생성된 시점 $t_{d}$으로 역전파를 수행하기 위해 현재 시점 $t_{k}$에서 $l$-step 지연된 시점 $t_{k-l}(\le t_{d})$에서의 상태추정치 $\hat{\mathbf{x}}_{k-l|k-l}$와 공분산행렬 $P_{k-l|k-l}$를 필요로 한다. 따라서 성능문제를 차치하더라도 이 기법은 보조센서의 시간지연 $T_{d}$가 길어질수록 더 많은 메모리를 요구하는 문제를 지니고 있다.

3.2. 시간불일치 보상 비동기 이종센서 융합필터 설계

3.2.1. 시간불일치 보상 필터 설계

만일 보조센서로부터 매 시점 현재 시점 정보가 지속적으로 제공된다면 다중경로의 영향을 받은 레이더 정보에 의해 추정성능 열화를 효과적으로 억제할 수 있다. 이에 착안하여, 본 논문에서는 그림 6과 같이 시간 불일치 보상 필터를 설계하여 전자광학센서의 시간지연 측정치 $\mathbf{y}_{d}=(\widetilde{\beta_{d}}, \widetilde{\epsilon_{d}})$로부터 현재 시점의 방위각 및 고각 추정치 $\hat{\psi}_{k}$와 $\hat{\theta}_{k}$를 매 시점 산출한 후, 이를 융합하는 방법을 OOSM 필터의 대안으로 제시한다.

그림 6. 제안하는 CIWS 정보융합 구조

Fig. 6. Proposed data fusion architecture of the CIWS

시간불일치 보상필터 설계를 위해, 식 (2)로부터 시간지연 된 각도 $\phi_{d}$$=(\beta$ 혹은 $\epsilon)$와 현재 시점 각도 $\phi =(\psi$ 혹은 $\theta)$ 간의 상관관계를 살펴보자. 시간지연 시스템의 입출력 관계는 1차 Pade 근사를 적용하여 전달함수의 형태로 모델링이 가능하다.

(15)

$G_{d}(s)=\frac{\phi_{d}(s)}{\phi(s)}= e^{-T_{d}s}\approx\frac{-0.5T_{d}s+1}{0.5T_{d}s+1}$

식 (15)에 역 라플라스 변환을 적용하면, 다음 결과를 얻는다.

(16)

$\begin{bmatrix}\dot{\beta}\\\dot{\epsilon}\end{bmatrix}= -\begin{bmatrix}\dot{\psi}\\\dot{\theta}\end{bmatrix}+\frac{2}{T_{d}}\left(\begin{bmatrix}\psi \\\theta\end{bmatrix}-\begin{bmatrix}\beta \\\epsilon\end{bmatrix}\right)$

한편, 코리올리 방정식으로부터 시선좌표계 표적 가속도 벡터는 다음과 같이 계산된다.

(17)

$\vec{A}^{L}=\frac{d^{2}\vec{R}^{L}}{dt^{2}} +\frac{d\vec{\omega}_{IL}^{L}}{dt}\times\vec{R}^{L} +\vec{\omega}_{IL}^{L}\times\left(\frac{d\vec{R}^{L}}{dt}+\vec{\omega}_{IL}^{L}\times\vec{R}^{L}\right)$

여기서 $\vec{R}^{L}\equiv C_{I}^{L}\vec{R}^{I}=\begin{bmatrix}r\\0\\0\end{bmatrix}$, $\vec{A}^{L}\equiv \begin{bmatrix}a_{x}\\ a_{y}\\ a_{z}\end{bmatrix}= C_{I}^{L}\frac{d^{2}\vec{R}}{dt}$, $\vec{\omega}_{IL}^{L}\equiv \begin{bmatrix}\omega_{x}\\\omega_{y}\\\omega_{z}\end{bmatrix}=\begin{bmatrix}-\psi\sin\theta \\\dot{\theta}\\\psi\cos\theta\end{bmatrix}$.

위 식의 두 번째 및 세 번째 행을 정리하면 다음과 같다.

(18)

$\begin{bmatrix}\dot{\omega}_{y}\\\dot{\omega}_{z}\end{bmatrix}=\begin{bmatrix}\frac{-2\dot{r}}{r}\omega_{y}-\frac{1}{r}a_{z}+\omega_{x}\omega_{y}\\\frac{-2\dot{r}}{r}\omega_{z}+\frac{1}{r}a_{y}-\omega_{x}\omega_{y}\end{bmatrix}$

대함유도탄 표적은 종말유도 단계에서 함정으로 접근하면서 수평면 상에서는 비례항법유도법칙을 적용하며, 수직면 상에서는 해면밀착 후 급상승/하강(pop-up & dive) 등의 회피기동을 수행하는 경우가 대부분이다 ^[24]. 따라서 CIWS의 교전 범위 내에서는 $|\psi |\ll 1, |\theta |\ll 1, |\dot{\psi}|\ll 1$이 만족된다. 따라서 식 (17)에서 $|\omega_{x}|\ll 1$ 이므로, 식 (18)을 근사하여 표적 방위각 고각에 대한 동특성 방정식을 유도할 수 있다 ^[25].

(19)

$\begin{bmatrix}\ddot{\psi}\\\ddot{\theta}\end{bmatrix}=\begin{bmatrix}-\frac{2\dot{r}}{r}\dot{\psi}+\dot{\theta}\tan\theta \\ -\frac{2\dot{r}}{r}\dot{\theta}\end{bmatrix}+\begin{bmatrix}a_{y}\\ -a_{z}\end{bmatrix}$

여기서 $(a_{y}, a_{z})$는 시선벡터에 수직방향인 표적 가속도이다.

대함유도탄 표적이 해면밀착 비행하는 경우 $|\dot{\theta}|\ll 1$이므로, 식 (16)과 식 (19)로부터 다음 운동모델을 정의할 수 있다.

(20)

$\dot{\boldsymbol{\eta}} =A \boldsymbol{\eta} + B \mathbf{e} , \mathbf{e}\sim N(0, \mathcal{E})$

여기서 $\boldsymbol{\eta}\equiv \begin{bmatrix}\boldsymbol{\eta}_{\psi}\\ \boldsymbol{\eta}_{\theta}\end{bmatrix}$, $A\equiv A_{\psi}\oplus A_{\theta}$, $B\equiv B_{\psi}\oplus B_{\theta}$,

$\boldsymbol{\eta}_{\psi}=\begin{bmatrix}\psi \\\dot{\psi}\\ a_{y}\\\beta\end{bmatrix}$, $A_{\psi}=\begin{bmatrix}0 & 1& 0 & 0 \\ 0 & -\frac{2\dot{r}}{r}&\frac{1}{r}& 0 \\ 0 & 0& 0 & 0 \\\frac{2}{T_{d}}& -1 & 0 & -\frac{2}{T_{d}}\end{bmatrix}$, $B_{\psi}=\begin{bmatrix}0 &0 \\ 0 & 0 \\1 & 0 \\0 & 1\end{bmatrix}$,

$\boldsymbol{\eta}_{\theta}=\begin{bmatrix}\theta \\\dot{\theta}\\ a_{z}\\\epsilon\end{bmatrix}$, $A_{\theta}=\begin{bmatrix}0 & 1& 0 & 0 \\ 0 & -\frac{2\dot{r}}{r}& -\frac{1}{r}& 0 \\ 0 & 0& 0 & 0 \\\frac{2}{T_{d}}& -1 & 0 & -\frac{2}{T_{d}}\end{bmatrix}$, $B_{\theta}=\begin{bmatrix}0 &0 \\ 0 & 0 \\1 & 0 \\0 & 1\end{bmatrix}$.

위의 식에서 연산자 $\oplus$는 행렬 직교합(direct sum)을 의미한다. 시스템 행렬 $A$에 사용된 시간지연 $T_{d}$는 전자광학센서에서 제공되는 시간 소인(time stamp)을 통해 계산 가능하며, 상대거리 $r$ 및 거리변화율 $\dot{r}$은 레이더의 거리-도플러 프로세싱 과정에서 상당한 정확도로 알 수 있는 값들이다.

일반적으로 샘플링 시간 $T_{s}$가 매우 작은 값을 가지므로, 샘플링 주기 동안 $r$, $\dot{r}$, $T_{d}$는 상수(piecewise constant)로 가정해도 무방하다. 이 경우, 이산화된 선형 시스템 모델을 다음과 같이 기술할 수 있다.

(21)

$\begin{cases} \boldsymbol{\eta}_{k+1} &=F_{k}\boldsymbol{\eta}_{k}+ \mathbf{e}_{k}, & \mathbf{e}_{k}\sim N(0, \mathcal{E}_{k})\\ \mathbf{y}_{d}&= C \boldsymbol{\eta}_{k}+ \boldsymbol{\nu}_{d}, & \boldsymbol{\nu}_{d}\sim N(0, R_{d})\end{cases}$

여기서 $F_{k}\equiv \exp\left(\int_{k T_{s}}^{(k+1)T_{s}}A(\tau)d\tau\right)= e^{AT_{s}}$, $\mathcal{E}_{k}\equiv \int_{0}^{T_{s}}e^{A\tau}B\mathcal{E}B^{T}(e^{A\tau})^{T}d\tau$,

$\mathbf{y}_{d}=\begin{bmatrix}\widetilde{\beta}_{d}\\\widetilde{\epsilon}_{d}\end{bmatrix}$, $C\equiv C_{\psi}\oplus C_{\theta}$, , $C_{\psi}= C_{\theta}=\begin{bmatrix}0 &0 & 0 & 1\end{bmatrix}$, $R_{d}\equiv diag(\sigma_{\beta}^{2}, \sigma_{\epsilon}^{2})$.

참고로 시스템 모델 (21)에 대한 가관측 행렬은 full-rank 이다. 따라서 시간지연 측정치 $\mathbf{y}_{d}$를 이용해 현 시점 고각 및 방위각 정보를 추정할 수 있다.

이상의 결과로부터 시간지연 측정치 $(\widetilde{\beta}, \widetilde{\epsilon})$으로부터 현재 시점 방위각 및 고각 $(\psi , \theta)$을 산출하는 시간지연 보상 필터는 선형 시스템 모델 (21)에 대한 선형 칼만필터가 된다.

측정치 갱신식

(22a)

$K_{k}=\gamma_{k|k-1}C^{T}(C\gamma_{k|k-1}C^{T}+ R_{d})^{-1}$

(22b)

$\gamma_{k|k}=(I- K_{k}C)\gamma_{k|k-1}$

(22c)

$\hat{\boldsymbol{\eta}}_{k|k}=\hat{\boldsymbol{\eta}}_{k|k-1}+ K_{k}(\mathbf{y}_{d}- C\hat{\boldsymbol{\eta}}_{k|k-1})$

여기서 $\hat{\boldsymbol{\eta}}_{k|k-1}$와 $\gamma_{k|k-1}$는 각각 상태변수 $\boldsymbol{\eta}_{k}$에 대한 사전추정치와 추정오차공분산을 의미한다.

시스템 전파식

(23a)

$\gamma_{k+1|k}= F_{k}\gamma_{k|k}F_{k}^{T}+ \mathcal{E}_{k}$

(23b)

$\hat{\boldsymbol{\eta}}_{k+1|k}= F_{k}\hat{\boldsymbol{\eta}}_{k|k}$

여기서 $\hat{\boldsymbol{\eta}}_{k|k}$와 $\gamma_{k|k}$는 각각 상태변수 $\boldsymbol{\eta}_{k}$에 대한 사후추정치와 추정오차공분산을 의미한다. $t_{k}$시점 방위각 및 고각 추정치는 $\hat{\psi}_{k}\equiv \hat{\boldsymbol{\eta}}_{k|k}(1)$, $\hat{\theta}_{k}\equiv \hat{\boldsymbol{\eta}}_{k|k}(5)$, 추정오차공분산은 $\hat{\sigma}_{\psi}^{2}\equiv \gamma_{k|k}(1, 1)$, $\hat{\sigma}_{\theta}^{2}\equiv \gamma_{k|k}(5, 5)$ 이다. 참고로 시스템모델 (21)에서 방위각 $\psi$과 고각 $\theta$ 운동방정식은 서로 독립적이므로, 수직면과 수평면으로 필터를 분리하여 설계해도 무방하다.

참고로 앞서와 동일한 방법으로 2차 Pade 근사를 적용하여 시간 불일치 보상 필터를 설계할 수 있다. 그러나 Pade 근사의 차수가 높아질수록 필터 차수 역시 증가하므로 표적 동특성 모델의 대역폭, 시간 지연 근사 정확도, 그리고 구현 복잡도를 종합적으로 고려하여 차수를 결정하는 것이 바람직하다. 그림 7은 시간 지연 시스템과 1차 및 2차 Pade 근사 전달함수의 위상 보드 다이어그램을 나타낸 것이다. 일반적으로 대함 유도탄 표적의 오토파일럿 대역폭은 1 Hz 미만이며, 이 주파수 영역에서는 1차와 2차 Pade 근사가 모두 시간 지연을 충분히 정확하게 모사한다. 따라서 1차 Pade 근사만으로도 2차 Pade 근사와 동일한 시간 지연 보상 효과를 기대할 수 있다.

그림 7. 시간지연 시스템의 위상 보드 선도

Fig. 7. Bode phase plot of time-delay system.