해상용 다기능레이더의 대기굴절 오차보정에 관한 연구

Ⅰ. 서 론

레이더 RF(Radio Frequency) 신호는 대기중에 방사될 경우, 고도에 따라 밀도가 다른 대기 특성에 의해 굴절되어 진행한다[1]. 이는 흔히 대기 굴절현상이라고 지칭되며, 굴절현상에 의해 굽어진 레이더 RF 신호의 진행방향은 거리와 고각 측정치에 편향된 오차를 야기한다. 대기굴절 현상에 의한 오차는 근거리 표적 탐지 목적의 레이더에서는 무시 가능한 수준이지만, 중장거리 표적 탐지를 목표로 하는 레이더는 대기굴절 오차를 고려해야만 한다[2].

레이더에서 대기 굴절오차를 보상하는 알고리즘에는 RT(Ray Tracing) 방식이 널리 쓰이고 있다. RT의 굴절 관계식은 PE(Parabolic Equation)에서 근사화를 통해 얻을 수 있으며, 이 수식은 Snell의 법칙으로도 알려져 있다[3]. Snell의 법칙은 굴절률이 다른 두 매개체의 진행방향에 대해 설명하고 있으므로, 대기 굴절현상은 고도에 따른 대기 굴절율을 Snell의 법칙에 적용하여 전파의 진행방향을 수학적으로 표현할 수 있다. 이와 같이 RT 방식을 사용하여 대기굴절 오차 보상 알고리즘들은 [4~9]에서 찾아볼 수 있다.

문헌 [4~7]의 ray를 반복적으로 계산하여 표적위치를 찾는 방식은 정확도가 높지만 계산 복잡도가 매우 높다. 계산 복잡도 때문에 미리 오차를 다항식으로 근사한 방식이 [8]에서 소개되었으나 고각 범위가 10도 이하로 한정적이며, 고도에 따른 실제 기상데이터의 적용이 불가능하다. 문헌 [9]에서는 고도에 따른 굴절지수와 표적 위치를 운동방정식으로 모델링하고 칼만필터를 활용해 오차를 추정하여, 기상데이터가 반영되며 실시간 계산이 될 수 있도록 하였다. 또한 문헌 [10]에서는 기상정보가 없을 때, 기체를 띄워 고도에 따른 굴절율을 실시간으로 모델링 하여 표적 위치를 추정하는 방식을 소개하였다. 하지만 [9][10]에서는 모두 고도값을 알고 있어야 한다는 전제조건이 있으므로 해상용 레이더에는 직접적으로 적용이 어렵다.

정리하여, 앞서 살펴보았던 기존의 RT 기반의 오차보상 알고리즘들은 계산 복잡도, 한정적인 빔 운용 범위, 표적의 고도정보를 필요로 하는 전제조건 때문에 실제 해상용 다기능 레이더가 단독으로 처리할 수 있는 알고리즘으로 선택되기 어렵다. 해상용 다기능레이더는 제공되는 기상정보로 해상에서 함정 단독으로 수백에서 수천의 표적에 대해 한정된 프레임 시간 내에 대기굴절 오차를 보상해야 한다. 따라서 본 논문에서는 해상용 다기능레이더에 적용가능한, 정확하고 단순한 방식의 오차 테이블 방식과 피팅 방식의 보상 방안들을 제안한다. 제공되는 기상정보를 통해 RT방식을 반복기법을 적용하여 운용 이전 혹은 운용 사이사이 미리 오차 테이블을 정확히 만든 후, 실시간 운용시에는 테이블을 활용하거나 근사화 수식을 활용하여 보상하는 방식이다.

본문 2장에서는 기존의 표적중심이 아닌 레이더 입장에서의 오차 계산 방식을 소개하고, 3장에서는 오차 보정 방안에 대해 소개한다. 4장에서는 실제 해상 환경 데이터를 통해 오차 보상 결과와 수행 시간에 대해 확인하고, 5장에서 결론을 통해 마무리한다.

Ⅱ. 레이더 측정치 기준 오차 테이블 계산

대기굴절 현상에 영향을 받는 레이더 측정치는 거리와 고각이기 때문에, 오차는 거리와 고각 측정치를 입력받아 계산되어야 한다. 아래는 측정거리R_radar와 측정고각 θ_radar에 대한 정의를 나타낸다.

위 수식에서 R_LOS는 레이더와 표적간의 실제 거리인 가시거리를 나타내고 ΔR은 거리오차를 나타낸다. 변수 θ_LOS는 레이더와 표적간의 실제 상대고각을 나타내고 Δθ는 고각 오차를 나타낸다.

거리오차와 고각 오차는 [7]에서 소개한 RT로 표현한 ray 궤적과 지구중심간의 기하관계를 사인법칙을 활용하여 계산할 수 있다. 그림 1은 [7]에서의 기하관계를 표현한 그림이다.

Fig. 1.

Geometry of altitude layers, ray, and earth

그림 1에서는 고도를 K개로 계층화하여 각 계층별 굴절지수 n_i, Ray의 기하거리 R_i, 고도에 따라 굴절된 거리 R_b, 굴절각 θ_i, 지구 중심과의 거리 l_i, 지구 중심각 ϕ_i으로 표현하였다. 위의 변수들은 다음 식 (3)~(8)과 같이 계산되며, 자세한 사항은 [7]을 참고할 수 있다.

2.1 Point approximation 방식의 오차계산

레이더 관점에서는 실제 표적의 고도값을 모르기 때문에 [7]의 알고리즘을 오차계산에 바로 사용할 수는 없다. 따라서 레이더 측정거리 R_radar와 [7]의 알고리즘에서 계산된 굴절된 거리 R_b가 일치되도록 반복계산을 수행한다. 이때, 고도 계층은 충분히 잘게 나누어 고도 미소 변위 Δh를 작게 할수록 정확한 결과를 얻을 수 있다.

Fig. 2.

Point approximation method

위 방식을 활용하여 균등한 입력값으로 오차테이블을 구성할 수 있다. 그림 3은 표준대기 굴절지수를 사용하고 레이더 측정거리 20~100km, 측정고각 0~20deg 범위의 거리 및 고각 오차 예시이다.

Fig. 3.

Range and elevation errors for 20~100km range and 0~20deg elevation measurements from point approximation

2.2 Ray spreading 방식의 오차계산

앞선 방식은 레이더가 측정한 거리, 고각에 해당하는 한 포인트에 대해 RT 알고리즘을 반복적으로 수행하는 방식으로, 충분히 고도 미소변위를 작게 하면 정확히 근사 추정할 수 있는 장점이 있다. 그러나 그림 3에서 보듯, 고고도-장거리로 갈수록 굴절각이 작아지기 때문에 고도 미소변위가 작아도 수평거리 이동량은 커지게 되고, 임계치 γ_th를 설정하기 어렵게 된다. 고도 미소 변위를 작게 하면 수평거리 이동량을 어느 정도 줄일 수 있지만, 이는 곧 계산량 증가로 이어진다.

Fig. 4.

Relationship between high altitude and long range (d1 ≪ dK)

하지만 한 포인트를 정확히 근사하는 방식 외에 Ray를 공간상에 펼친 후, Ray 포인트 지점에 대한 오차를 이용하여 레이더 측정치 오차를 계산할 수 있다. 그림 5는 0도부터 20도까지 1도 간격으로 Ray를 방사했을 때의 예시이다.

Fig. 5.

Example of rays(red line) from radar([0,0]), target(red circle), line of sight trajectory between radar and target(green stitch line), points on the rays with same altitude(blue circle), and the earth(blue bottom)

공간상에 뿌려진 m번째 Ray의 고각은 아래와 같이 표현된다.

$\[ {\theta }^m={\theta }^1+\left(m-1\right)\Delta \theta \left(1\le m\le M\right) \]$

(9)

[7]의 알고리즘을 활용하여, m번째 Ray의 j번째 고도계층의 거리오차 $$ \Delta R_j^m$$, 고각오차 $$ \Delta {\theta }_j^m$$는 다음과 같이 계산할 수 있다.

$\[ \Delta R_j^m=R_{b\left(j\right)}^m-R_{LOS\left(j\right)}^m \]$

(10)

$\[ \Delta {\theta }_j^m={\theta }^m-{\theta }_{LOS\left(j\right)}^m \]$

(11)

$\[ R_{b\left(j\right)}^m=\sum _{i=1}^j{n}_i{R}_i \]$

(12)

$\[ \begin{array}{c}{R}_{LOS\left(j\right)}^m=\hfill \\ \sqrt{{\left(l_0^m\right)}^2+{\left(l_j^m\right)}^2-2\left(l_0^m\right){\left(l_j^m\right)}^{2}cos\left(\sum _{i=1}^j{\varphi }_i^m\right)}\hfill \end{array} \]$

(13)

$\[ {\theta }_{LOS\left(j\right)}^m=\text{acos}\left(\frac{l_j^m}{R_{LOS\left(j\right)}^m}\cdot \text{sin}\left(\sum _{i=1}^j{\varphi }_i^m\right)\right) \]$

(14)

이 방식은 Ray가 존재하는 위치에 대해서만 오차를 계산하여 상대적으로 계산이 빠르지만 레이더 측정거리, 즉 굴절거리는 고도에 대한 함수이기 때문에 사용자가 거리축으로 균등하게 오차 테이블을 만들기 어렵다. 따라서 거리축으로 보간(Interpolation)을 수행하여 그림 6과 같이 그림 3의 조건과 동일하게 오차를 계산할 수 있다.

Fig. 6.

Range and Elevation errors for 20~100km range and 0~20deg elevation measurements from Ray spreading

두 방식간의 오차계산 결과를 비교하기 위해 그림 7에 두 방식의 결과의 차이를 계산하였다. 두 방식의 비교 결과 고각이 낮을수록 차이가 존재하지만, 그 수준이 2% 이하로 충분히 무시할 수 있을 정도로 유사한 결과를 보였다.

Fig. 7.

Comparison between two methods

하지만 두 방식의 계산 시간은 큰 차이를 보였는데, 표 1을 통해 간단히 비교하였다. 단순한 실행 시간 비교지만 계산량은 확연히 두 번째 방식이 적음을 확인할 수 있다.

Table 1.

Comparison of simulation time

Ⅲ. 대기굴절 오차보상 알고리즘

앞장에서 구성한 오차테이블은 측정치 보정에 두 가지 방식으로 활용 될 수 있다. 첫 번째는 오차 테이블을 메모리에 저장하여, 매 측정치마다 찾아 사용하는 look-up table 방식이다. 두 번째는 오차 테이블을 2차원 피팅하여 함수화한 뒤, 매 측정치를 함수 입력값으로 사용하여 결과를 도출하는 피팅 방식이다.

3.1 Look-up table 방식

거리범위 R_T와 고각범위 θ_T로 구성된 거리오차 테이블 T_Rerr와 고각오차 테이블 t_θerr을 아래와 같이 정의할 수 있다.

$\[ \begin{array}{c}{R}_T=R_{\text{min}}+\left({idx}_R-1\right)\Delta R{}_T ,\hfill \\ \left(1\le {idx}_R\le \frac{\left(R_{\text{max}}-R_{\text{min}}\right)}{\Delta R_T}+1\right)\hfill \end{array} \]$

(15)

$\[ \begin{array}{c}{\theta }_T={\theta }_{\text{min}}+\left({idx}_{\theta }-1\right)\Delta {\theta }_T ,\hfill \\ \left(1\le {idx}_{\theta }\le \frac{\left({\theta }_{\text{max}}-{\theta }_{\text{min}}\right)}{\Delta {\theta }_T}+1\right)\hfill \end{array} \]$

(16)

$\[ \mathit{\Delta }\mathit{R}=T_{Rerr}\left({idx}_R,{idx}_{\theta }\right) \]$

(17)

$\[ \mathit{\Delta }\mathit{\theta }=T_{\theta err}\left({idx}_R,{\theta }_T{idx}_{\theta }\right) \]$

(18)

위에서 idx_R와 idx_θ는 각각 거리범위와 고각범위의 인덱스값을 의미하고, ΔR_T와 Δθ_T는 각각 거리와 고각 간격을 의미하고, ΔR와 Δθ거리와 고각 오차값들의 집합을 의미한다.

레이다에서 측정된 거리 R_radar을 포함하는 오차 테이블들의 인덱스 idx_R1, idx_R2는 아래 방식으로 찾을 수 있다.

$\[ diff=R_T\left({idx}_R\right)-R_{radar} \]$

(19)

(20)

$\[ {idx}_{R2}={idx}_{R1}+1 \]$

(21)

고각 인덱스 idx_θ1, idx_θ2 역시 동일한 방법으로 찾을 수 있다.

Fig. 8.

Finding error indexes from radar measurements

찾은 인덱스를 통해 고각오차는 아래와 같이 내분점 찾기로 계산할 수 있다.

$\[ \begin{array}{c}\Delta R=T_{Rerr}\left({idx}_{R1},{idx}_{\theta 1}\right)+...\hfill \\ \frac{T_{Rerr}\left({idx}_{R1},{idx}_{\theta 1}\right)-T_{Rerr}\left({idx}_{R2},{idx}_{\theta 2}\right)}{R_T\left({idx}_{R1},{idx}_{\theta 1}\right)-R_T\left({idx}_{R2},{idx}_{\theta 2}\right)}\times ...\hfill \\ \left(R_{radar}-R_T\left({idx}_{R1},{idx}_{\theta 1}\right)\right)\hfill \end{array} \]$

(22)

$\[ \begin{array}{c}\Delta \theta =T_{\theta err}\left({idx}_{R1},{idx}_{\theta 1}\right)+...\hfill \\ \frac{T_{\theta err}\left({idx}_{R1},{idx}_{\theta 1}\right)-T_{\theta err}\left({idx}_{R2},{idx}_{\theta 2}\right)}{{\theta }_T\left({idx}_{R1},{idx}_{\theta 1}\right)-{\theta }_T\left({idx}_{R2},{idx}_{\theta 2}\right)}\times ...\hfill \\ \left({\theta }_{radar}-{\theta }_T\left({idx}_{R1},{idx}_{\theta 1}\right)\right)\hfill \end{array} \]$

(23)

최종 보정된 레이더 측정치(R_compen, θ_compen)는 아래와 같이 계산된다.

$\[ R_{compen}=R_{radar}-\Delta R \]$

(24)

$\[ {\theta }_{compen}={\theta }_{radar}-\Delta \theta \]$

(25)

3.2 피팅 방식

앞선 Look-up table 기반의 오차 보정 방식은 테이블에서 매번 레이더 측정치에 해당하는 인덱스들을 찾아야하기 때문에, 실시간 처리에 영향을 줄 수 있다. 따라서 오차테이블을 2차원 피팅을 통해 함수화하여 사용하면 위 문제를 해결할 수 있다. 아래 그림들은 그림 6 시나리오에서 오차 테이블을 거리축 2차, 고각축 3차의 Poly23 선형모델로 피팅한 예시를 나타낸다.

Fig. 9.

Example of 2D fitting of error table

위 예시에서는 거리오차의 RMSE 결과는 0.22m, 고각오차 피팅의 RMSE 결과는 4.73e-04deg 수준으로 충분히 사용가능한 결과를 얻을 수 있다. 이처럼 오차 테이블을 사용자가 원하는 피팅 모델로 함수화하여, 그 계수를 레이더 메모리에 저장하여 오차를 보상할 수 있다.

Ⅳ. 시뮬레이션

본 장에서는 앞에서 살펴본 오차 테이블 계산 방식과 오차 보상 방식을 실제 해상환경 데이터에 적용하여 확인해본다.

4.1 기상데이터

국제전기통신연합(ITU, International Telecommunication Union)는 고도에 따른 굴절지수를 고도에 따른 기압, 기온, 습도의 함수로 아래와 같이 표현했다[11].

$\[ \begin{array}{c}N\left(h\right)= \hfill \\ 77.6\frac{P\left(h\right)}{T\left(h\right)}-5.6\frac{e\left(h\right)}{T\left(h\right)}+3.75e5\frac{e\left(h\right)}{T{\left(h\right)}^2}\hfill \end{array} \]$

(26)

$\[ n=1+N\cdot {10}^{-6} \]$

(27)

미국 Wyoming 대학에서는 매일 0시, 12시의 고도에 따른 기상데이터를 제공하며, 그림 10은 2021년 10월 21일 포항의 12시 기상데이터를 나타낸다.

Fig. 10.

Atmospheric data at Pohang, KOREA, 2021-10-21

식 (26)을 이용하면 고도에 따른 굴절지수를 그림 11과 같이 계산할 수 있다.

Fig. 11.

Refractivity index with altitude at Pohang, KOREA, 2021-10-21

4.2 오차테이블

그림 12는 그림 11의 굴절지수를 사용하고, Ray spreading 방식을 사용하여 계산한 오차 테이블이다. 거리 범위는 10km 간격으로 1~300km까지, 고각은 1도 간격으로 90도까지 설정하였다. 고도 100km 이하까지만 계산했기 때문에 테이블의 빈 부분을 확인할 수 있다.

Fig. 12.

Error tables for range of 1~300km and elevation of 1~90deg measurements

포항의 해당 시점의 대기굴절에 의한 거리오차는 최대 70m 이상, 고각 오차는 최대 0.25도 이상임을 확인 할 수 있으며, 실제 표적이 존재하는 저고도 구간에 오차가 큰 부분이 집중분포됨을 확인할 수 있다.

4.3 오차보상 결과

위에서 계산한 오차 테이블을 사용하여 실제 포항 해상환경에서 측정한 레이더 측정치를 Look-up table 방식으로 보상하였다. 레이더 측정치는 표 2를 참고한다.

Table 2.

Summary of radar measurements

그림 13(a)는 포항 연안에서 해상방향으로 측정한 레이더 거리, 고각 측정치이고 13(b)는 수평거리-고도 그래프이다. 기회표적은 고각 6도 이하에 존재하며, 수백 km에 해당하는 궤적을 가지고 있다. 굴절오차 보상 이전의 결과는 파란색으로 표현되었고, 보상 후 결과는 붉은색으로 표현되었다. 오차 보상 전-후 차이는 최대 거리 40m, 고각 0.2도, 고도 570m 차이가 나는 것을 확인할 수 있다.

Fig. 13.

Result of refraction compensation

위 데이터에 측정된 측정치(Plot) 수는 탐지방위범위 50도에 대해 12분 동안 28779개였으므로, 전방위 기준으로 계산하면 1초에 280여개의 측정치가 발생할 것으로 예상할 수 있다. 위 해상환경에서는 기회표적 수가 적은 경우였으므로, 실제 표적 수가 많을 경우는 1초에 수천에서 수만개의 측정치가 발생할 수 있다.

표 3에서는 Plot 개수에 따라 계산에 사용된 시간 중 최소값을 표현하였다. 표 3에서 Plot 수 증가량에 비례하여 계산 시간차이 역시 증가함을 확인하여, 직접 오차 테이블을 찾는 것보다 2차원 피팅한 방식의 계산량이 적음을 확인하였다.

Table 3.

Simulation environment

Ⅴ. 결 론

본 논문에서는 해상용 다기능레이더의 대기굴절 오차 보상을 위한 오차계산법을 살펴보고 실제 적용 가능한 보정기법에 대해 소개하였다. 소개하는 오차계산법은 기존방식과 비교할 때, 레이더에서 측정치에 직접적으로 오차를 보정할 수 있도록 개선되었다. 실제 해상환경에서 측정된 기상 데이터로 대기굴절 오차를 확인할 때, 거리 300km 기준 거리오차 약 75미터, 고각 오차 약 0.27도 가량 발생함을 확인하였으며, 기존에 널리 사용되는 반복 방식보다 Ray를 펼쳐 계산하는 방식이 계산 소요시간이 10배 이상 적음을 확인하였다. 오차 보정기법은 테이블 방식과 피팅 방식을 소개하였으며, 레이더 처리 성능에 따라 선택 가능하다.

해상용 다기능레이더는 짧은 운용시간에도 처리해야할 측정치 개수가 매우 많으므로, 제안되는 방식은 그 유용성이 클 것으로 예상된다.

References

• Mark A. Richards, James A. Scheer, and William A. Holm, "Principles of modern radar", SciTech, pp. 130-137, 2010. [https://doi.org/10.1049/SBRA021E]
• B. Salski, et al., "2-D numerical analysis M&S tool of radar wave propagation considering atmospheric density profile", Technical report, AM Technologies, 2019.
• M. Levy, "Parabolic equation methods for electromagnetic wave propagation", The Institution of Engineering and Technology, pp. 4-19, Jun. 2000. [https://doi.org/10.1049/PBEW045E_ch2]
• Hyunwook Moon, et. al., "Development of Exponential Model of Korea for Improved Altitude Estimation Performance of High-Altitude Target at Radar System", KJKIEES, Vol. 23, No. 7, pp. 831-839, Jul. 2012. [https://doi.org/10.5515/KJKIEES.2012.23.7.831]
• BAE system, "RF Refraction correction", unpublished.
• Agee, B., et al., "NASA Ray Trace algorithm', unpublished.
• Bassem R. Mahafza, "Radar signal analysis and processing using MATLAB", CRC Press, pp. 44-47, 2009.
• P. E. Schmid, "Atmospheric tracking errors at S-and C-band frequencies", NASA Technical Note TN D-340, p.9, Aug. 1966.
• Stefano fortunati, et al., "Correction of refracted propagation effects for airborne radar tracking", IEEE Trans. Aerospace and Electronic Systems, Vol. 49, No. 1, pp. 20-41, Jan. 2013. [https://doi.org/10.1109/TAES.2013.6404089]
• H. LEE, et al., "Target positioning with surveillance radar by the estimation of atmospheric refractivity profile", IEEE Trans. AEROSPACE AND ELECTRONIC SYSTEMS, Vol. 56, No. 3, pp. 2245-2255, Jun. 2020. [https://doi.org/10.1109/TAES.2019.2945385]
• ITU-R P.453-11, "The radio refractive index: its formula and refractivity data", International Telecommunication Union, pp. 2-3, Jul. 2015.