포스트

(Radar) 12. 좌표계와 좌표변환 — 안테나가 잰 각도가 지도 위 위경도가 되기까지

레이다가 잴 수 있는 것은 자기 정면 기준의 거리 하나와 각도 두 개뿐인데 지휘관은 위경도를 봐야 한다 — 그 사이를 메우는 좌표계들의 원점과 축 규약, 회전행렬이 삼각함수 덧셈정리 한 줄에서 나오는 과정, 안테나 좌표계의 세 얼굴(Polar·Cartesian·UV)과 위상배열이 굳이 UV를 쓰는 사정, INS가 바깥을 전혀 보지 않고도 자세를 아는 방법과 스트랩다운·짐벌의 갈림, NED와 ENU·ECEF와 ECI·타원체와 지오이드의 차이, 함정 고정형 안테나를 놓고 5단계 변환을 설계하고 되돌리는 과정, 그리고 그것을 C로 짜서 왕복 오차 10⁻⁹ m 까지 확인하기.

(Radar) 12. 좌표계와 좌표변환 — 안테나가 잰 각도가 지도 위 위경도가 되기까지

레이다가 “20 km 앞, 왼쪽으로 30도” 라고 말할 때 전시기에는 “북위 36.3610도, 동경 127.5220도” 가 찍혀야 한다. 둘은 같은 표적인데 말하는 사람이 다르다. 그 사이를 메우는 것이 좌표변환이다.

이 글은 그 통역 과정을 처음부터 끝까지 따라간다. 앞부분에서 좌표계 여섯 가지를 하나씩 소개하고, 중간에서 함선 레이다 하나를 놓고 변환 경로를 설계하며, 마지막에 그것을 C로 구현해 숫자로 확인한다. 사전 지식은 삼각함수와 행렬 곱셈 정도면 충분하고, 그마저도 필요한 곳에서 다시 설명한다.


0. 시작하기 전에

0.1 좌표계는 왜 여러 개나 필요한가

레이다 이야기를 하기 전에, 물체 하나를 두 사람이 가리키는 경우를 먼저 보자. 두 사람은 서 있는 자리도 다르고 보고 있는 방향도 다르다.

물체 하나를 두 관측자가 서로 다른 자리와 방향에서 보는 그림 그림 0-1. 물체는 하나인데 A 는 오른쪽, B 는 왼쪽이라고 한다. 서 있는 자리와 보는 방향이 다르면 같은 물체의 위치를 적는 숫자가 달라진다.

관측자 A 는 물체가 자기 앞에서 오른쪽으로 21°, 9 m 거리에 있다고 하고, 관측자 B 는 왼쪽으로 29°, 7 m 거리에 있다고 한다. 둘 다 맞는 말이다. 물체는 한자리에 그대로 있고, 두 사람 사이에 다른 것은 서 있는 자리와 보고 있는 방향 두 가지뿐이다. 서 있는 자리처럼 재는 기준이 되는 점을 원점(origin), 보고 있는 방향처럼 기준이 되는 방향을 축(axis)이라 하고, 이 둘을 미리 정해 놓은 약속이 좌표계(coordinate frame)다. 레이다와 전시기 사이에서도 같은 일이 일어난다.

레이다 안테나는 자기가 지금 어디를 보고 있는지 스스로는 모른다. 안테나가 아는 것은 딱 세 가지뿐이다.

  1. 전파가 갔다 오는 데 걸린 시간 → 거리
  2. 자기 정면에서 좌우로 몇 도
  3. 자기 정면에서 위아래로 몇 도

배가 어느 방향을 향하는지, 지금 위도가 몇 도인지는 안테나의 관심 밖이다.

반대로 지휘관이 보는 전시기는 지도다. 지도에 점을 찍으려면 위도·경도가 필요하다. 그런데 위경도는 지구 전체를 기준으로 한 값이고, 안테나가 준 “왼쪽으로 30도”와는 아무 관계가 없다.

이 간극을 메우려면 중간 단계가 필요하다. 안테나의 “왼쪽 30도”를 배 기준으로 바꾸고, 배 기준을 북쪽·동쪽 기준으로 바꾸고, 그걸 다시 지구 기준으로 바꾼다. 좌표계가 여러 개인 이유는 딱 이것이다.

핵심 좌표계는 “어디”를 말하는 언어다. 표적은 하나뿐인데 언어가 여러 개라서 화자를 바꿀 때마다 통역이 필요하다. 그 통역이 좌표변환이다.

같은 표적, 다섯 가지 대답

3절에서 실제로 계산할 표적 하나를 미리 보자. 아래 다섯 줄은 전부 같은 점을 가리킨다.

말하는 사람표적이 어디 있냐고 물으면
안테나내 정면에서 왼쪽 30°, 20 km
함선선수(뱃머리) 기준 오른쪽 60°, 20 km
바다 위 이 지점남쪽으로 5 198 m, 동쪽으로 19 297 m, 위로 10 m
지구지구 중심에서 (−3 132 054, +4 078 527, +3 760 546) m
지도북위 36.3610°, 동경 127.5220°, 고도 41.3 m

다섯 문장 모두 맞다. 다른 것은 기준점(원점)과 기준 방향(축) 뿐이다.

0.2 좌표계를 읽는 세 가지 질문

새 좌표계를 만나면 아래 셋만 확인하면 된다. 1절의 모든 좌표계를 이 형식으로 설명한다.

질문왜 중요한가
① 원점이 어디인가?“0, 0, 0”이 물리적으로 어느 지점인가. 평행이동 계산의 기준이 된다.
② 축이 어디를 향하는가?특히 z가 위인지 아래인지. 항법 분야는 z를 아래로 두는 관습이 있어 처음엔 거의 반드시 헷갈린다.
③ 무엇에 붙어 있는가?배에 붙어 있으면 배가 흔들릴 때 좌표계도 같이 흔들린다. 표적이 가만히 있어도 좌표값이 변한다.

0.3 변환은 두 동작뿐

좌표변환이라고 하면 복잡한 수학처럼 들리지만 실제로 하는 일은 두 가지뿐이다.

좌표변환의 두 동작인 회전과 평행이동 그림 0-2. 왼쪽은 원점이 같고 축만 θ 만큼 어긋난 경우이고, 오른쪽은 축 방향이 같고 원점만 t 만큼 떨어진 경우다.

  • 회전(rotation) — 기준 방향이 다를 때. “내 정면”과 “배의 정면”이 다르면 회전이 필요하다.
  • 평행이동(translation) — 기준 점이 다를 때. 안테나가 무게중심에서 30 m 떨어져 있으면 그만큼 옮겨준다.

모든 단계는 이 둘 중 하나거나 둘 다다.

회전을 눈으로 보기

표적은 가만히 있다. 움직이는 것은 자(ruler) 다.

0.1절의 표에서 안테나가 “내 정면에서 왼쪽 30°, 20 km” 라고 한 표적을 위에서 내려다본 그림으로 보자. 안테나는 배의 오른쪽 옆면인 우현(starboard)을 보도록 달려 있다.

같은 표적을 안테나 기준과 함선 기준으로 읽은 결과 그림 0-3. 두 그림에서 배와 표적의 위치는 완전히 같다. 달라진 것은 거리를 재는 축뿐이다.

(a) 는 안테나 기준으로 읽은 것이다. 20 km 를 정면 방향과 옆 방향으로 나눠 적으면 표적은 정면 쪽으로 17 320 m 나아가고 거기서 왼쪽으로 10 000 m 벗어난 곳에 있다 (나누는 방법은 1.1.2절에서 다룬다). 옆 방향은 오른쪽을 + 로 적기로 하면 왼쪽 10 000 m 는 −10 000 이고, 읽은 값은 (정면 +17 320, 오른쪽 −10 000) 이다. (b) 는 같은 표적을 함선 기준으로 읽은 것이다. 선수 방향과 우현 방향으로 나눠 적으면 (선수 +10 000, 우현 +17 320) 이다.

안테나의 “정면”이 함선의 “우현”이기 때문에 정면 쪽 17 320 m 는 그대로 우현 쪽 17 320 m 가 된다. 우현을 보고 서면 왼쪽이 선수이고 오른쪽이 배의 뒤쪽인 선미(stern)다. 그래서 안테나의 왼쪽 10 000 m 는 선수 쪽 10 000 m 가 되고, 오른쪽 −10 000 이던 값이 선수 +10 000 으로 적힌다. 두 숫자가 자리를 바꾸고 부호 하나가 뒤집혔다. 축이 90도 돌아가 있으면 좌표는 이렇게 바뀐다.

이 그림은 회전만 보려고 안테나가 무게중심에서 떨어져 있는 만큼은 빼고 그렸다. 그 차이는 아래 「평행이동은 그냥 더하기」에서 더한다.

회전행렬은 어디서 나오는가

90도처럼 딱 떨어지면 숫자를 바꿔치기하면 되지만 45도, 30도는 그렇게 안 된다. 그때 필요한 게 sin·cos이고, 그것을 표로 정리해 놓은 것이 회전행렬이다. 어렵게 외울 것이 아니라 삼각함수 덧셈정리 한 줄에서 바로 나온다.

2차원 회전 유도 그림 0-4. 회전행렬의 유도. 점을 극좌표로 쓰고 각도만 θ 만큼 더하면 끝이다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
① 원래 점 (x, y) 를 극좌표로 쓴다
     x = r·cos α ,   y = r·sin α

② θ 만큼 더 돌리면 각이 (α + θ) 가 된다
     x′ = r·cos(α + θ) = r·cosα·cosθ − r·sinα·sinθ
     y′ = r·sin(α + θ) = r·sinα·cosθ + r·cosα·sinθ
                          (삼각함수 덧셈정리)

③ r·cosα = x,  r·sinα = y  를 되돌려 넣는다
     x′ = x·cos θ − y·sin θ
     y′ = x·sin θ + y·cos θ

④ 이 두 줄을 표로 적은 것이 Rz(θ) 다
     ⎡ x′ ⎤   ⎡ cos θ   −sin θ ⎤   ⎡ x ⎤
     ⎢    ⎥ = ⎢                 ⎥ × ⎢   ⎥
     ⎣ y′ ⎦   ⎣ sin θ    cos θ ⎦   ⎣ y ⎦

3차원에서는 “어느 축 둘레로 도느냐”에 따라 이 2×2 블록이 들어가는 자리만 달라진다. z축 둘레 회전이면 z는 그대로 두고 x·y 만 섞이므로, 위 2×2 를 왼쪽 위에 넣고 나머지를 단위행렬로 채우면 된다.

1
2
3
4
5
6
        ⎡ 1     0        0    ⎤        ⎡  cos β   0   sin β ⎤        ⎡ cos γ  −sin γ   0 ⎤
Rx(α) = ⎢ 0   cos α   −sin α  ⎥ Ry(β)= ⎢    0     1     0   ⎥ Rz(γ)= ⎢ sin γ   cos γ   0 ⎥
        ⎣ 0   sin α    cos α  ⎦        ⎣ −sin β   0   cos β ⎦        ⎣   0       0     1 ⎦

x축 둘레 회전                    y축 둘레 회전                    z축 둘레 회전
(= roll)                        (= pitch)                       (= yaw)

우리 문제의 90도 회전을 넣어 보면 위에서 눈으로 본 결과와 정확히 같다.

1
2
3
4
5
         ⎡ 0  −1   0 ⎤       ⎡ +17 320.5081 ⎤     ⎡ +10 000.0000 ⎤
Rz(90°) =⎢ 1   0   0 ⎥   ×   ⎢ −10 000.0000 ⎥  =  ⎢ +17 320.5081 ⎥
         ⎣ 0   0   1 ⎦       ⎣        0.0000 ⎦     ⎣        0.0000 ⎦

                             정면·오른쪽·아래     선수·우현·아래

안테나가 실제로 쓰는 축은 안테나 면을 기준으로 잡아서 위의 정면·오른쪽·아래와 순서도 방향도 다른데, 그 축으로 적은 값을 정면·오른쪽·아래로 바꿔 읽는 일은 1.2절에서 다룬다.

행렬을 쓰는 이유 세 가지

  1. x·y·z 세 축을 한 번에 처리한다.
  2. 회전을 여러 번 해야 할 때 행렬끼리 곱해서 하나로 합칠 수 있다. roll·pitch·yaw 세 회전도 곱 하나로 뭉쳐진다.
  3. 되돌아갈 때는 전치(행과 열 바꾸기)만 하면 된다. 역행렬을 따로 구할 필요가 없다.

세 번째 성질은 회전행렬이 직교행렬이기 때문에 성립한다. “직교행렬”이란 Cᵀ · C = I(단위행렬)를 만족하는 행렬이고, 회전은 길이와 각도를 바꾸지 않으므로 항상 이 성질을 갖는다. 3절의 코드는 되돌아가는 쪽의 회전행렬을 모두 전치로 만들었고, 갔다가 되돌아온 값이 입력과 같은지를 3.1.6절에서 확인한다.

평행이동은 그냥 더하기

회전이 끝났으면 원점 차이만 더해 주면 된다.

1
2
3
4
5
6
7
p_함선  =  Rz(90°) × p_안테나  +  레버암
                                 └─ (−30, 0, −10) m
                                     x = −30 : 30 m 뒤  (x는 앞이 +)
                                     z = −10 : 10 m 위  (z는 아래가 +)

        =  (+10 000, +17 320.5, 0)  +  (−30, 0, −10)
        =  (9 970, 17 320.5, −10) m

이 오프셋을 레버암(lever arm) 이라고 부른다. 무게중심에서 센서까지 뻗은 “지렛대 팔”이라는 뜻이다.

더한 결과 (9 970, 17 320.5, −10) 은 무게중심에서 잰 표적 위치다. 안테나가 무게중심보다 30 m 뒤에 있으니 무게중심에서 재면 표적은 선수 쪽으로 30 m 덜 나가 있고, 표적은 안테나와 같은 높이에 있는데 그 안테나가 10 m 위에 있으니 표적도 무게중심보다 10 m 높은 곳에 있다.

0.4 이 문서의 표기 약속

좌표변환 버그의 대부분은 수식이 아니라 약속에서 나온다. 같은 기호가 문헌마다 반대 의미로 쓰이기 때문이다. 아래 여섯 가지를 먼저 못 박고 시작한다.

항목약속
z축 방향동체 좌표계(배 기준)와 NED 좌표계(북·동·아래 기준)는 z 가 아래다. “10 m 위”는 z = −10. 안테나 좌표계만 안테나 면을 기준으로 잡아서 x 가 왼쪽, y 가 위, z 가 정면(보어사이트, boresight)이다.
방위각 부호안테나 정면에서 왼쪽으로 잰 각이 +. 위에서 내려다보면 반시계 방향이다.
회전행렬 이름C_ned_from_body = “동체 좌표를 NED 좌표로 바꾸는 행렬”. 이름에 방향을 명시한다.
역변환C_body_from_ned = (C_ned_from_body)ᵀ — 전치하면 된다.
오일러각 순서3-2-1 (yaw → pitch → roll). C_ned_from_body = Rz(yaw)·Ry(pitch)·Rx(roll)
각도 단위계산·저장은 라디안, 화면·인터페이스만 도(degree). 경계에서 한 번만 변환한다.

왜 이름에 방향을 넣는가 C_bn 이라는 변수명은 라이브러리마다 “b→n”이기도 하고 “n→b”이기도 하다. 방향이 반대여도 컴파일은 되고 값도 그럴듯하게 나온다. 그래서 이름에 방향을 박아 넣는 것이 가장 값싼 방어책이다. 3.1.6절의 왕복 확인은 이런 실수가 정변환과 역변환 가운데 한쪽에만 끼어든 경우를 걸러 낸다. 양쪽에 함께 끼어든 경우는 왕복으로 드러나지 않아서 손으로 계산한 값과도 맞춰 본다.


1. 좌표계 여섯 가지

레이다 시스템에서 주로 쓰는 좌표계를 하나씩 살펴본다. 각 좌표계마다 ① 원점 ② 축 ③ 붙은 곳을 먼저 밝히고, 그다음 “왜 필요한가”와 “언제 쓰는가”를 설명한다.


1.1 안테나(Antenna) 좌표계: UV, Polar, Cartesian

안테나 좌표계는 레이다가 표적을 실제로 측정하는 곳이다. 모든 것이 여기서 시작한다.

질문답
① 원점안테나 배열면 한가운데(위상중심)
② 축x = 왼쪽   y = 위   z = 보어사이트(안테나가 정면으로 바라보는 방향)
③ 붙은 곳배 (여기서 다루는 예제는 고정형이라 볼트로 고정 — 안테나가 따로 돌지 않는다)

보어사이트(boresight) 라는 말이 계속 나온다. 안테나 면에 수직으로 뻗은 방향, 즉 “안테나가 똑바로 바라보는 방향”이다. 손전등으로 치면 빛이 가장 센 한가운데 방향이다.

축은 안테나 면을 기준으로 잡았다. 전파를 내보내고 받는 작은 소자들이 늘어선 면(배열면)이 x–y 평면이고, 그 면에 수직인 보어사이트가 z 다. 오른손 좌표계(오른손의 엄지·검지·중지를 서로 직각으로 펴서 차례로 x·y·z 로 삼는 좌표계)에서 y 를 위, z 를 정면으로 두면 x 는 왼쪽으로 정해진다. 여기서 왼쪽·오른쪽은 안테나 뒤에 서서 보어사이트를 따라 내다볼 때를 기준으로 한다. 뒤에 나올 동체 좌표계와 NED 는 z 가 아래인데 안테나 좌표계만 z 가 정면이다. 두 방식의 축을 맞추는 일은 1.2절에서 한다.

같은 방향을 가리키는 데에 세 가지 표현이 쓰인다. 셋 다 같은 정보를 담고 있고, 목적에 따라 골라 쓴다.

1.1.1 Polar (극좌표) — 측정값이 나오는 형태

레이다가 실제로 재는 것은 거리 하나와 각도 둘이다. 그래서 원측정값은 자연스럽게 극좌표다.

성분기호의미
시선거리R (Slant Range)안테나에서 표적까지의 직선 거리. 전파 왕복시간 × 빛의 속도 ÷ 2
방위각Az (Azimuth)보어사이트에서 왼쪽으로 잰 각 (위에서 내려다보면 반시계)
고각El (Elevation)안테나 수평면에서 위로 잰 각

안테나 극좌표 R, Az, El 의 정의 그림 1-1. 방위각 Az 는 위에서 본 그림 (a) 에서 보어사이트를 기준으로 왼쪽으로, 고각 El 은 옆에서 본 그림 (b) 에서 수평면을 기준으로 위쪽으로 잰다.

뒤에서 쓸 예제 입력값은 R = 20 km, Az = 30°, El = 0° 이다. 즉 “안테나 정면에서 왼쪽으로 30도 방향, 20 km 떨어진 곳, 높낮이는 정면과 같음”이라는 뜻이다.

주의 — Az/El 정의는 하나가 아니다 문헌에 따라 방위각을 재는 기준축이나 회전 순서가 다르다. 방위각의 부호를 반대로 잡아 오른쪽을 + 로 두는 규약도 흔하다. 그 규약으로 읽으면 같은 입력 Az = 30° 가 보어사이트의 왼쪽이 아니라 오른쪽을 가리키고, 이 예제에서는 표적이 실제 위치에서 약 20 km 떨어진 다른 점에 찍힌다 (2 × 20 km × sin 30°). 어느 정의를 쓰는지는 문서로 못 박아야 하고, 가장 확실한 방법은 아래 1.1.3의 방향코사인 수식을 함께 적어두는 것이다. 이 문서는 “보어사이트에서 왼쪽으로 Az, 그다음 위로 El” 순서(El-over-Az 계열)를 쓴다.

1.1.2 Cartesian (직교좌표) — 계산하기 좋은 형태

극좌표는 사람이 읽기 좋지만, 회전행렬을 곱하려면 x·y·z 세 성분으로 풀어써야 한다. 그래서 계산의 첫 단계는 항상 극좌표 → 직교좌표다.

바꾸는 식은 직각삼각형의 삼각비를 두 번 쓴 것이다. 삼각비는 빗변의 길이와 각 하나를 알 때 밑변은 빗변 × cos, 높이는 빗변 × sin 이라는 관계다.

먼저 고각으로 위아래 성분을 떼어 낸다 (그림 1-2의 삼각형 ①). 표적까지의 거리 R 을 빗변으로, 고각 El 을 빗변과 밑변 사이의 각으로 하는 직각삼각형을 세우면, 높이 R × sin(El) 이 위 방향 성분 y 다. 밑변 R × cos(El) 은 표적을 안테나 수평면(x–z 평면)에 수직으로 내린 점까지 안테나에서 잰 거리다.

다음으로 그 밑변을 방위각으로 한 번 더 나눈다 (삼각형 ②). 이번에는 R × cos(El) 이 빗변이고 방위각 Az 가 빗변과 보어사이트 사이의 각이다. 보어사이트를 따라 놓인 변이 밑변이므로 R × cos(El) × cos(Az) 가 z 이고, 거기서 왼쪽으로 뻗은 변이 높이이므로 R × cos(El) × sin(Az) 가 x 다.

극좌표를 직교좌표로 푸는 두 직각삼각형 그림 1-2. 삼각형 ① 에서 거리 R 을 고각으로 나눠 높이 y 와 밑변을 얻고, 삼각형 ② 에서 그 밑변을 방위각으로 나눠 z 와 x 를 얻는다.

1
2
3
4
5
6
7
8
x = R × cos(El) × sin(Az)      ← 왼쪽 방향 성분
y = R × sin(El)                ← 위 방향 성분
z = R × cos(El) × cos(Az)      ← 보어사이트 방향 성분

역변환:
R  = √(x² + y² + z²)
Az = atan2(x, z)
El = atan2(y, √(x² + z²))

Az 가 + 이면 x 가 +, El 이 + 이면 y 가 + 다. 각이 커지는 방향과 축의 + 방향을 같게 잡았기 때문에 세 식 어디에도 음수 부호가 붙지 않는다. 역변환은 같은 두 삼각형을 거꾸로 읽은 것이다. Az 는 삼각형 ② 의 두 변 x 와 z 에서, El 은 삼각형 ① 의 높이 y 와 밑변 √(x² + z²) 에서 구한다. atan2(a, b) 는 탄젠트가 a/b 인 각을 구하되 a 와 b 의 부호를 따로 보아 −180° 에서 +180° 까지의 각을 돌려주는 함수다.

예제 값을 넣으면:

1
2
3
x = 20000 × 1 × sin30° = +10 000.0000 m
y = 20000 × sin0°      =        0.0000 m
z = 20000 × 1 × cos30° = +17 320.5081 m

고각이 0° 라 위 방향 성분은 없고, 보어사이트 쪽 성분이 17 320.5 m, 왼쪽 성분이 10 000 m 다. 0.3절에서는 같은 표적을 정면 +17 320 m, 오른쪽 −10 000 m 로 읽었다. 여기서는 왼쪽을 + 로 두어 부호가 바뀌었고 적는 순서가 (왼쪽, 위, 보어사이트) 가 되었을 뿐, 가리키는 위치는 같다.

왜 atan2 이고 왜 asin 이 아닌가 역변환에서 El = asin(y/R) 로 써도 수학적으로는 맞다. 하지만 두 가지 실용적인 문제가 있다. 첫째, 반올림 때문에 인자가 1.0000000000000002 가 되면 asin 은 NaN(Not a Number)을 내놓고, NaN 은 아무 경고 없이 뒤 계산 전체로 번진다. 둘째, 인자가 ±1에 가까울수록(고각이 ±90°에 가까울수록) asin 의 정밀도가 급격히 나빠진다. atan2(분자, 분모) 형태는 두 문제가 모두 없다. 3절 구현도 atan2 를 쓴다.

1.1.3 UV (방향코사인 / 사인공간) — 위상배열이 실제로 다루는 형태

UV 좌표는 표적 방향을 안테나 면에 드리운 그림자로 표현한 것이다.

표적 방향을 가리키는 화살표를 길이 1로 줄인 다음, 그 화살표가 안테나 면 위에 만드는 그림자의 가로 위치가 u, 세로 위치가 v 다. 남은 성분 w 는 보어사이트 방향이다.

식으로 쓰면 1.1.2의 직교좌표 세 성분을 거리 R 로 나눈 것이다. 거리는 지워지고 방향만 남는다. 세 값은 표적 방향이 x·y·z 축과 각각 이루는 각의 코사인이어서 방향코사인(direction cosine)이라고 부른다.

1
2
3
4
5
6
7
8
9
10
11
u = x / R = cos(El) × sin(Az)     ← 안테나 면의 가로(왼쪽) 성분
v = y / R = sin(El)               ← 안테나 면의 세로(위)   성분
w = z / R = cos(El) × cos(Az)     ← 보어사이트 성분

항상 성립:  u² + v² + w² = 1
가시영역 :  u² + v² ≤ 1   ← 원판 밖은 전파가 나가지 않는 영역

역변환:
w  = √(1 − u² − v²)
Az = atan2(u, w)
El = atan2(v, √(u² + w²))

세 성분을 제곱해 더한 값이 늘 1 이라는 것은 표적 방향이 반지름 1 인 구 위의 한 점이라는 뜻이다. 안테나는 면의 앞쪽으로만 전파를 내보내므로 그 구에서 w 가 0 이상인 앞쪽 반구만 쓴다. 그 반구를 안테나 면에 수직으로 내리면 반지름 1 인 원판이 되고, (u, v) 는 늘 이 원판 안에 있다. 이 원판을 가시영역(visible region)이라 한다.

역변환에는 u, v 두 값만 있으면 된다. 길이가 1 이라는 조건에서 w 를 먼저 구하고, 나머지 두 줄은 1.1.2의 역변환에 x, y, z 대신 u, v, w 를 넣은 식이다. w 에 양의 제곱근만 쓰는 것은 앞쪽 반구만 쓰기 때문이다.

예제 값(Az 30°, El 0°)이면 u = 0.5, v = 0, w = 0.866 이다.

반지름 1 인 반구 위의 표적 방향과 안테나 면에 내린 그림자 u, v 그림 1-3. 길이 1 인 화살표의 끝점은 반구 위에 있고, 그 끝점을 안테나 면에 수직으로 내린 자리가 (u, v), 면에서 떨어진 높이가 w 다.

그림은 세 성분이 모두 보이도록 Az 20°, El 35° 인 방향을 그렸고, 이때 u = 0.280, v = 0.574, w = 0.770 이다. 예제 표적은 El 이 0° 라서 v = 0 이고, 그림자가 u 축 위의 u = 0.5 인 자리에 놓인다.

UV를 쓰는 이유 — 위상배열의 사정

현대 레이다는 대부분 위상배열(phased array) 안테나다. 접시가 돌아가는 대신, 안테나 면에 붙은 수백~수천 개의 작은 소자에 신호를 아주 조금씩 시간차를 두고 보내서 빔 방향을 전자적으로 바꾼다. 이 시간차(위상차)를 얼마로 줘야 하는가가 핵심 계산인데,

  • 각도(Az, El) 로는 이 관계가 복잡하다.
  • u, v 로는 정확히 비례한다. 소자 위치에 u, v를 곱하기만 하면 된다.

빔을 만들고 조향하는 계산이 UV 공간에서 덧셈·곱셈 수준으로 단순해지는 것, 이것이 UV를 쓰는 실질적인 이유다. 조금 더 구체적으로 보면,

1
2
3
4
5
소자 (m, n) 에 주어야 할 위상
    ψ_mn = −k · ( m·d_x·u₀  +  n·d_y·v₀ )
             └ k = 2π/λ (파수),  d_x, d_y = x 방향(가로)·y 방향(세로) 소자 간격,  (u₀, v₀) = 빔을 보낼 방향

  → 소자 번호에 u₀, v₀ 를 곱하기만 하면 된다. 각도(Az, El) 로는 이렇게 단순해지지 않는다.

그 밖에도 UV 공간에서만 성립하는 편리한 성질이 있다.

  • 빔의 굵기가 조향각과 무관하게 일정하다. 각도 공간에서는 옆으로 조향할수록 빔이 굵어지고 이득이 떨어진다(배율 1/cos θ, 60°에서 2배·3 dB 손실). 고정형 다면 배열에서 면 하나가 담당할 각도 범위를 정하는 근거가 이 숫자다.
  • grating lobe(격자엽) 위치가 u₀ ± λ/d 로 단순하게 나온다. grating lobe 란 원하지 않는 방향으로 생기는 두 번째 빔이고, 소자 간격 d 가 넓으면 가시영역 안으로 들어와 엉뚱한 방향의 표적을 잡는다. 이를 막는 조건도 UV 로 쓰면 d/λ < 1/(1 + |u|_max) 한 줄이다.
  • 각도 오차를 추정하는 모노펄스 계산이 u, v 에 대해 선형이다.
세 표현의 역할 정리
표현쓰이는 곳형태
Polar신호처리 출력, 운용자 화면, 외부 인터페이스R, Az, El
Cartesian좌표변환 계산 (회전행렬을 곱하는 곳)x, y, z
UV빔 조향·빔 형성, 배열 내부 계산u, v (+ w)

이 설계에서는 신호처리가 Polar로 값을 넘겨주고(2.1.1의 조건 ④), 데이터처리는 이를 Cartesian으로 바꿔서 변환 체인을 태운다. UV는 안테나 내부(빔 제어)에서 쓰이므로 데이터처리 체인에는 직접 등장하지 않는다. 또 UV는 방향만 있고 거리가 없어서 표적 위치를 옮기는 이 체인에는 넣지 않았다.


1.2 동체(Body/Platform) 좌표계

배에 실린 모든 장비의 공통 기준이다. 안테나도, 항법장치도, 무장도 결국 여기로 모인다. 레이다를 실은 배나 항공기를 플랫폼(platform)이라고 불러서 이 좌표계를 플랫폼 좌표계라고도 한다. 1.3절의 INS 설명에 나오는 “플랫폼 좌표계”는 센서를 얹어 늘 수평으로 유지하는 받침대에서 나온 말이라 이것과 뜻이 다르다.

질문답
① 원점배의 무게중심(또는 도면상 정해진 기준점)
② 축x = 선수(앞)   y = 우현(오른쪽)   z = 아래
③ 붙은 곳배 (배와 함께 움직이고 함께 기운다)

축 배치의 머리글자를 따서 FRD(Forward-Right-Down)라고 부른다.

동체 좌표계 FRD 와 자세각 세 가지 그림 1-4. 동체 좌표계의 축과 roll·pitch·yaw. z축은 지면 속으로 들어간다.

자세각 세 가지

동체 좌표계가 로컬 좌표계(1.4)에 대해 얼마나 돌아가 있는지를 나타내는 값이다.

이름기호뜻부호
rollφ좌우로 기우뚱 (횡동요)우현이 내려가면 +
pitchθ앞뒤로 끄덕 (종동요)선수가 들리면 +
yawψ선수가 어느 방위를 향하는가 (heading)진북에서 시계방향 +

세 회전을 합치는 순서는 3-2-1(yaw → pitch → roll)이 표준이다.

1
C_ned_from_body = Rz(yaw) · Ry(pitch) · Rx(roll)

순서가 왜 중요한가 회전은 곱셈 순서를 바꾸면 결과가 달라진다. 책상 위 책을 “옆으로 90도 돌린 뒤 세우기”와 “세운 뒤 옆으로 90도 돌리기”는 다른 자세가 된다. 그래서 순서를 문서로 고정해야 한다.

안테나 설치와 레버암

고정형 안테나는 배에 볼트로 고정되어 있으므로 안테나 ↔ 동체 관계가 항상 같다. 즉 그 회전행렬은 설치 도면이 정하는 상수다.

1
2
C_body_from_antenna = Rz(설치 방위) · Ry(백틸트) · A
p_body = C_body_from_antenna × p_antenna + 레버암
  • 설치 방위: 선수 기준으로 안테나가 몇 도 돌아가 있는가 (예제에서는 90°)
  • 백틸트: 안테나 면을 뒤로 젖혀 보어사이트를 위로 들어올린 각 (예제에서는 0°)
  • 축 맞춤 행렬 A: 안테나 좌표의 세 성분을 정면·오른쪽·아래 순서로 다시 적는 행렬 (각도가 들어가지 않아 어떤 설치에서도 같은 값)
  • 레버암: 무게중심에서 안테나까지의 위치 차이 (예제에서는 30 m 뒤, 10 m 위)

앞의 셋은 축 방향을 맞추는 회전이고, 레버암은 원점을 맞추는 평행이동이다.

축 맞춤, 설치 방위, 레버암과 백틸트 그림 1-5. (a) 안테나 축을 정면·오른쪽·아래로 다시 읽고, (b) 설치 방위만큼 돌린 다음, (b)·(c) 에 표시한 레버암만큼 옮긴다. 백틸트는 (c) 에서 안테나 정면이 수평보다 위로 들린 각이다.

축 맞춤 행렬 A

안테나 좌표계(1.1)는 축을 안테나 면에 맞춰 x = 왼쪽, y = 위, z = 보어사이트로 잡았다. 동체 좌표계는 x = 선수, y = 우현, z = 아래다. 앞을 가리키는 축이 안테나에서는 세 번째(z)이고 동체에서는 첫 번째(x)라서, 두 좌표계는 설치 각도를 따지기 전에 x·y·z 를 어느 방향에 붙였는지부터 다르다. 그래서 안테나 좌표를 동체의 앞·오른쪽·아래와 같은 순서, 곧 안테나의 정면·오른쪽·아래로 먼저 옮겨 적는다. 이 일을 하는 행렬을 축 맞춤 행렬이라 부르고 A 로 적겠다.

A 는 세 성분의 자리를 옮기고 부호를 맞춘다. 정면은 보어사이트이므로 z 를 그대로 가져오고, 오른쪽은 왼쪽의 반대이므로 −x, 아래는 위의 반대이므로 −y 다.

1
2
3
4
5
A × (x, y, z) = ( z, −x, −y ) = ( 정면, 오른쪽, 아래 )

    ⎡  0   0   1 ⎤
A = ⎢ −1   0   0 ⎥
    ⎣  0  −1   0 ⎦

원소가 0 과 ±1 뿐인 상수 행렬이다. 설치 방위나 백틸트 같은 설치 각도가 들어가지 않는다. 표적도 안테나도 움직이지 않고, 같은 위치를 축만 바꿔 다시 읽는다. 90도 회전 두 번을 곱한 Rz(−90°) · Rx(−90°) 도 같은 행렬이 되고, 3절 코드는 이 곱으로 A 를 만든다.

앞의 식 C_body_from_antenna = Rz(설치 방위) · Ry(백틸트) · A 에서 A 는 세 행렬 가운데 맨 오른쪽, 곧 벡터 p_antenna 에 가장 가까운 자리에 있다. 벡터에 가장 먼저 곱해진다는 뜻이다. 먼저 곱하는 까닭은 그다음에 곱하는 두 회전이 정면·오른쪽·아래 순서를 전제로 하기 때문이다. Rz(설치 방위) 는 세 번째 축인 아래 축 둘레로, Ry(백틸트) 는 두 번째 축인 오른쪽 축 둘레로 도는 회전이다. 자세각에서 yaw 가 z축 둘레 회전이고 pitch 가 y축 둘레 회전인 것과 같은 방식이다. A 를 곱하지 않은 안테나 좌표에 Rz 를 바로 곱하면, 그 좌표의 세 번째 축은 보어사이트라서 보어사이트를 축으로 안테나 면을 돌리는 회전이 된다. 설치 방위가 뜻하는 회전은 그것이 아니다.

예제 값이 바뀌는 과정

예제 측정값 R = 20 km, Az = 30°, El = 0° 를 1.1.2의 식으로 풀면 안테나 직교좌표는 왼쪽으로 10 000 m, 위로 0 m, 보어사이트 쪽으로 17 320.5081 m 다. 이 값에 식의 오른쪽 행렬부터 차례로 곱하고 마지막에 레버암을 더한다. 백틸트가 0° 라서 Ry(0°) 는 단위행렬이고 값을 바꾸지 않으므로 표에서는 뺐다.

순서한 일값 [m]세 성분의 뜻
처음안테나 직교좌표 p_antenna(+10 000.0000, 0.0000, +17 320.5081)왼쪽, 위, 보어사이트
①A 를 곱한다(+17 320.5081, −10 000.0000, 0.0000)정면, 오른쪽, 아래
②Rz(90°) 를 곱한다(+10 000.0000, +17 320.5081, 0.0000)선수, 우현, 아래
③레버암 (−30, 0, −10) 을 더한다(+9 970.0000, +17 320.5081, −10.0000)선수, 우현, 아래 (p_body)

①에서는 표적 위치를 정면으로 17 320.5 m, 오른쪽으로 −10 000 m 라고 다시 적었다. 표적이 정면에서 왼쪽에 있으므로 오른쪽을 + 로 두면 음수가 된다.

②에서는 설치 방위 90° 만큼 돌렸다. 안테나는 위에서 내려다볼 때 선수에서 시계방향으로 90° 돌아가 우현을 보고 있으므로, 안테나의 정면이 배의 우현이고 안테나의 왼쪽이 배의 선수 쪽이다. 배 위에서 우현을 바라보고 서면 왼쪽이 선수인 것과 같다. 그래서 정면 17 320.5 m 는 우현 17 320.5 m 가 되고, 오른쪽 −10 000 m, 곧 왼쪽 10 000 m 는 선수 +10 000 m 가 된다. 0.3절에서 Rz(90°) 에 곱해 본 벡터가 ①의 결과다.

③에서는 원점을 안테나에서 무게중심으로 옮겼다. 안테나는 무게중심보다 30 m 뒤, 10 m 위에 있고 동체 좌표로 적으면 (−30, 0, −10) 이다. 이 값을 더하면 선수 방향 값이 30 m 줄고 아래 방향 값은 −10 m 가 된다. 표적이 무게중심보다 10 m 높다는 뜻이다. 레버암은 동체 축으로 잰 값이므로 회전을 모두 마친 뒤에 더한다.

세 행렬을 미리 곱해 두면 예제의 회전은 상수 행렬 하나가 된다.

1
2
3
                                               ⎡ 1   0   0 ⎤
C_body_from_antenna = Rz(90°) · Ry(0°) · A  =  ⎢ 0   0   1 ⎥
                                               ⎣ 0  −1   0 ⎦

행을 위에서부터 읽으면 배의 선수(x)는 안테나의 왼쪽(x)이고, 배의 우현(y)은 안테나의 정면(z)이며, 배의 아래(z)는 안테나의 위(y)의 반대 방향이다.

다른 관습과 충돌 주의 항공·항법 분야는 FRD(앞-오른쪽-아래)를 쓰지만, 로보틱스(ROS 등)는 FLU(앞-왼쪽-위)를 쓴다. 앞을 가리키는 x축이 같아서 roll 은 그대로지만 y축과 z축이 반대라 pitch 와 yaw 의 부호가 뒤집힌다. 외부 시뮬레이터나 시각화 도구를 붙일 때 가장 자주 사고가 나는 경계이므로, 변환 어댑터를 한 군데로 모아두는 것이 좋다.


1.3 INS 좌표계

INS(Inertial Navigation System, 관성항법장치) 는 자이로와 가속도계로 배의 위치·속도·자세를 스스로 계산해 내는 장비다. 좌표변환 관점에서 INS는 회전과 평행이동에 필요한 값을 공급해 주는 곳이다.

질문답
① 원점INS가 설치된 위치
② 축INS 안의 자이로·가속도계 3축이 향하는 방향
③ 붙은 곳배 (센서를 선체에 고정하는 스트랩다운형일 때)

표에 적은 것은 좁은 의미의 INS 좌표계다. INS 문헌에 함께 나오는 다른 프레임은 아래 「”INS 좌표계”는 사실 하나가 아니다」에서 정리한다.

INS는 무엇을 재는가

INS 안에는 두 종류의 센서가 각각 3개씩 들어 있다. 둘 다 자기 자신에게 일어나는 일만 잰다.

센서재는 것여기서 얻는 것
자이로(gyroscope) × 3각속도 — 초당 몇 도씩 돌고 있는가시간에 대해 적분하면 자세(roll·pitch·yaw)
가속도계(accelerometer) × 3비력(specific force) — 중력을 뺀 나머지 가속도두 번 적분하면 속도와 위치

핵심은 INS가 바깥 세상을 전혀 보지 않는다는 점이다. GPS처럼 위성 신호를 받지도 않고 지형을 관측하지도 않는다. 출발할 때의 위치·자세를 알려주면, 그 뒤로는 자기가 느낀 각속도와 가속도만 계속 더해서 “지금 어디 있고 어느 쪽을 보고 있는지”를 추측해 나간다.

그래서 두 가지 성질이 따라온다.

  • 장점: 외부 신호가 필요 없다. 전파방해를 받지 않고, 갱신 속도가 매우 빠르다(보통 100 Hz 이상).
  • 단점: 오차가 시간이 갈수록 쌓인다. 자이로의 미세한 편향(bias)이 적분되어 자세 오차가 되고, 그 자세 오차가 다시 가속도 성분을 잘못 나누어 위치 오차를 키운다. 그래서 실제 시스템은 GPS 같은 외부 기준과 결합(GPS/INS 통합)해 이 누적을 잡아 준다.

초기 정렬(initial alignment)이 필요한 이유 INS는 “출발점”을 스스로 알 수 없다. 그래서 출항 전에 가만히 세워 두고 중력 방향으로 수평을 잡고(leveling), 지구 자전 방향을 감지해 북쪽을 찾는(gyrocompassing) 절차를 거친다. 이 정렬이 끝나야 자세 출력이 의미를 갖는다. 정렬이 부정확하면 그 오차가 그대로 좌표변환 3단계에 실린다.

스트랩다운과 짐벌형

INS에는 역사적으로 두 방식이 있고, 이 차이가 “플랫폼 좌표계”라는 말의 의미를 갈라놓는다.

스트랩다운 INS와 짐벌형 INS 그림 1-6. 짐벌형은 수평 플랫폼이 물리적으로 존재하고, 스트랩다운은 계산상으로만 존재한다.

 짐벌형 (gimbaled)스트랩다운 (strapdown)
구조센서가 짐벌 위에 얹혀 있고, 모터가 짐벌을 돌려 늘 수평을 유지한다센서가 동체에 직접 볼트로 고정되어 배와 함께 기운다
수평 유지기계가 물리적으로 한다계산이 한다 (DCM 또는 쿼터니언을 계속 갱신)
“플랫폼 좌표계”실제로 존재하는 물리적 프레임소프트웨어 안의 상태값일 뿐
특징정밀하지만 크고 무겁고 비싸다작고 튼튼하며 저렴하다. 현대 함정 INS는 사실상 전부 이쪽

문서에서 “플랫폼 좌표계”라는 말을 볼 때 그것이 물리적 장치를 가리키는지 소프트웨어 상태를 가리키는지 확인해야 한다. 스트랩다운에서는 후자다.

“INS 좌표계”는 사실 하나가 아니다

INS 문헌에는 네 개의 프레임이 함께 나온다. INS는 이 넷의 관계를 계속 유지하는 장치다.

기호이름INS에서의 역할
i관성 프레임 (ECI)자이로·가속도계가 실제로 감지하는 기준. 뉴턴 법칙이 그대로 성립하는 곳
e지구고정 (ECEF)위치를 저장하고 GNSS(GPS)와 맞춰보는 기준
n항법 프레임 (NED)속도와 자세를 출력하는 기준
b센서 프레임 (Body)원측정값(각속도·가속도)이 나오는 프레임

즉 좁은 의미의 “INS 좌표계”는 b 프레임 — INS 안의 자이로·가속도계 3축이 실제로 향하는 방향이다.

왜 동체 좌표계와 따로 구분하는가

이상적으로는 INS의 센서 축과 배의 동체 축이 딱 맞아야 하지만, 실제로는

  • 설치 오차: 장착할 때 아주 미세하게 어긋난다 (보통 수 mrad = 수천분의 몇 라디안)
  • 선체 변형: 큰 배는 파도·하중·온도에 따라 선체가 휜다. INS 위치와 안테나 위치의 상대 자세가 시간에 따라 변한다
  • 위치 차이: INS가 무게중심에 있지 않으면 그만큼의 레버암이 또 생긴다

이 어긋남을 보정하는 절차를 정렬(alignment) 이라고 한다.

오차 감각 잡기 1 mrad은 약 0.0573°다. 거리 × 0.001 이 그대로 옆방향 오차가 된다. 10 km → 10 m, 100 km → 100 m, 300 km → 300 m. 즉 “1도도 안 되는 아주 작은 각”이 먼 거리에서는 큰 위치 오차가 된다.

이 상황에서는

INS가 플랫폼 무게중심과 축이 일치하게 설치되어 있다면?

이 조건 덕분에 위의 복잡한 이야기가 전부 사라진다.

1
INS 좌표계  =  동체(Body) 좌표계        (설치 오차 0, 레버암 0)

따라서 이번 설계에서는 INS 좌표계를 독립된 변환 단계로 두지 않는다. INS는 좌표계가 아니라 값의 공급원 역할만 한다.

INS가 주는 값어디에 쓰이는가
roll, pitch, yaw3단계 회전행렬 C_ned_from_body
위도, 경도, 고도4단계 로컬→ECEF 회전과 평행이동
(속도)표적 추적 필터·도플러 보정 (이 글의 범위 밖)

실무에서 확인해야 할 것 INS가 주는 heading이 진북 기준인지 자북·격자북 기준인지, 고도가 타원체고인지 해발고도인지, 자세 데이터에 유효시각(언제 측정한 값인가) 이 붙어 있는지. 셋 다 좌표변환 결과를 통째로 틀리게 만들 수 있다.


1.4 NED, ENU 좌표계

배가 지금 떠 있는 바다 위 한 지점에서 북쪽·동쪽 방향을 축으로 삼은 좌표계다. “로컬 좌표계” 또는 “국지수평 좌표계”라고 부르는 것이 보통 이것이다.

질문답
① 원점배의 현재 위치 (또는 지정한 고정 기준점)
② 축NED: N=북, E=동, D=아래  /  ENU: E=동, N=북, U=위
③ 붙은 곳바다 위의 한 지점 — 배와 함께 이동하지만 함께 회전하지는 않는다

③번이 이 좌표계의 존재 이유다. 배가 아무리 흔들려도 북쪽은 계속 북쪽이다.

NED와 ENU의 차이

두 좌표계는 같은 개념의 축 순서·부호만 다른 형제다.

 첫째 축둘째 축셋째 축주로 쓰는 분야
NED북 (N)동 (E)아래 (D)항공우주·항법·무기체계
ENU동 (E)북 (N)위 (U)측지·측량·GIS·로보틱스

NED 와 ENU 축 비교 그림 1-7. 두 좌표계는 같은 정보를 담는다. 축 순서와 z 부호만 다르다.

1
2
NED 값이 (N, E, D) 일 때  →  ENU 값은 (E, N, −D)
ENU 값이 (E, N, U) 일 때  →  NED 값은 (N, E, −U)

NED를 쓰는 이유: 동체 좌표계(FRD: 앞-오른쪽-아래)와 축의 의미가 맞아떨어져서 자세 회전행렬이 자연스럽게 연결된다. 이번 설계도 NED를 쓴다. ENU를 쓰는 이유: 고도가 +z라 사람이 읽기 직관적이다. 지도·측량 쪽 데이터와 붙일 때 편하다.

어느 쪽이든 하나로 통일하는 것이 중요하다. 섞이면 남북이 뒤집히거나 고도 부호가 반대가 된다.

“아래”의 정확한 의미

D축(아래)은 지구 중심 방향이 아니다. 지구는 완전한 구가 아니라 적도 쪽이 부푼 타원체라서, 지표면에 세운 수직선(법선)은 지구 중심을 지나지 않는다.

1
2
3
4
5
6
7
8
       북극
        │
        │      ● 내가 서 있는 지점 P
        │     ↙ 법선 방향(= 진짜 '아래'의 반대)
        │   ↙
        ●───────────  적도
      지구중심
              ↖ 지구 중심 방향 (P에서 보면 법선과 살짝 다르다)

두 방향의 차이는 위도 45° 부근에서 최대 약 0.19° 다. 작아 보이지만 이는 지표에서 약 21 km에 해당하는 각도이므로, 구형 지구로 근사한 코드에 측지 위도를 그대로 넣으면 그만큼의 계통 오차가 전 체인에 실린다.

로컬 평면의 한계 — 지구는 둥글다

로컬 좌표계는 평평한 판이다. 원점 근처에서는 문제가 없지만, 멀어질수록 실제 지구 표면이 그 판 아래로 휘어 내려간다.

원점에서의 거리지구 표면이 평면보다 내려간 양
5 km약 2 m
10 km약 8 m
20 km약 31 m
50 km약 196 m
100 km약 783 m

(3절에서 실제로 확인한다. 20 km 표적의 고도가 로컬 평면 기준 10 m인데 측지 고도로는 41.3 m가 나오는 이유가 바로 이 31.3 m다.)


1.5 ECEF 좌표계, ECI 좌표계

여기서부터는 지구 전체를 다루는 좌표계다. 원점이 지구 중심으로 옮겨간다.

1.5.1 ECEF — Earth-Centered, Earth-Fixed (지구중심 지구고정)

질문답
① 원점지구의 질량중심
② 축Z = 북극 방향   X = 적도면과 그리니치 자오선이 만나는 방향   Y = 오른손계를 완성하는 방향
③ 붙은 곳지구 — 지구와 함께 하루 한 바퀴 자전한다
1
2
3
4
5
6
7
8
9
        Z (북극)
        │
        │      ● 어떤 지점
        │    /
        │  /
        ●────────── Y (동경 90°)
       /
     /
   X (적도 ∩ 그리니치 자오선, 즉 경도 0°)

왜 필요한가. 로컬 좌표계에서 위경도로 곧장 넘어갈 수는 없다. 로컬 좌표계의 원점(배의 위치)이 지구 어디인지는 지구 기준 좌표로만 표현되기 때문이다. ECEF는 그 다리 역할을 한다. 또 여러 배·여러 레이다의 결과를 합칠 때 서로 자세가 다르니 공통 언어가 필요한데, ECEF가 그 역할을 한다.

특징. 미터 단위의 평범한 직교좌표라서 더하기·빼기가 그냥 된다. 다만 값이 600만 m 규모라 소수점 아래를 다루기엔 정밀도 여유가 적다. (가까운 두 점의 상대 위치는 ECEF가 아니라 로컬 좌표계에서 계산해야 한다.)

실현. 실제로 쓰는 ECEF는 WGS-84라는 국제 표준으로 정의되어 있고, GPS가 방송하는 좌표도 이것이다.

로컬에서 ECEF 로 가는 회전

로컬 좌표 (N, E, D) 는 원점(배의 위치)에서 북으로 N, 동으로 E, 아래로 D 만큼 간다는 뜻이다. 그런데 북·동·아래가 지구 축 X·Y·Z 로 어느 쪽인지는 원점이 지구 어디에 있느냐에 따라 다르다. 적도에서는 북이 Z축과 나란하지만 북극 가까이에서는 북이 Z축과 거의 수직이다. 그래서 로컬 벡터를 ECEF 로 옮기려면 원점에서의 북·동·아래를 먼저 지구 축 기준으로 적어 두어야 한다.

지구 위 한 점 P 의 북·동·아래 방향과 ECEF 축 그림 1-8. 점 P 에서 북은 자오선을 따라 북극 쪽, 동은 위도선을 따라 동쪽, 아래는 타원체 면에 수직으로 안쪽을 가리킨다. P 의 위도와 경도가 바뀌면 세 화살표가 X·Y·Z 축에 대해 놓이는 방향도 바뀐다. 그림은 지구를 구로 그려서 아래 화살표가 지구 중심을 향하지만 실제 타원체에서는 중심을 조금 비껴 간다.

세 방향을 길이가 1 인 벡터, 곧 단위벡터(unit vector)로 잡아 ECEF 성분으로 적으면 아래와 같다. 그림 1-8 의 점 P 가 로컬 원점이고, 그 위도를 φ, 경도를 λ 로 적었다. 그림에서 N, E, D 라고 이름 붙인 세 화살표가 이 단위벡터다. 식에서는 좌표값 (N, E, D) 와 섞이지 않게 북, 동, 아래로 적는다.

1
2
3
북   = ( −sin φ·cos λ,   −sin φ·sin λ,    cos φ )
동   = ( −sin λ,          cos λ,          0     )
아래 = ( −cos φ·cos λ,   −cos φ·sin λ,   −sin φ )

식이 맞는지는 적도와 그리니치 자오선이 만나는 점(φ = 0, λ = 0)을 넣어 보면 바로 확인된다. X축이 지표를 뚫고 나오는 자리라서 세 방향이 좌표축과 나란해진다. 적도 위에서는 타원체 면의 수직선이 지구 중심을 지나기 때문에, 이 점에서는 아래가 지구 중심 쪽이기도 하다.

1
2
3
북   = (  0,  0,  1 )      ← Z축 방향. 북극 쪽
동   = (  0,  1,  0 )      ← Y축 방향. 동경 90° 쪽
아래 = ( −1,  0,  0 )      ← X축의 반대 방향. 지구 중심 쪽

로컬 벡터는 북·동·아래 단위벡터에 좌표값 N, E, D 를 각각 곱해 더한 것이다. 세 단위벡터를 열로 세운 행렬을 만들어 두면 그 계산이 행렬 곱 한 번이 된다.

1
2
3
4
N × 북 + E × 동 + D × 아래  =  C_ecef_from_ned × (N, E, D)

C_ecef_from_ned = [ 북  동  아래 ]        ← 세 벡터를 열로 세운 3×3 행렬
C_ned_from_ecef = (C_ecef_from_ned)ᵀ      ← 되돌아갈 때는 전치. 세 벡터가 행으로 눕는다

같은 행렬을 회전 두 번으로도 만들 수 있다.

1
C_ecef_from_ned = Rz(λ) · Ry(−90° − φ)

북·동·아래가 X·Y·Z축과 겹쳐 있는 상태에서 시작한다고 생각하면 된다. Y축 둘레로 −90° 돌리면 북이 Z축을, 아래가 X축의 반대쪽을 향해서 위에서 확인한 적도 위의 점과 같아진다. 같은 방향으로 φ 만큼 더 돌리면 위도 φ 인 점의 방향이 되고, 이어서 Z축 둘레로 λ 만큼 돌리면 경도 λ 인 자오선 위로 옮겨 간다. 3절 코드는 이 방식으로 행렬을 만든다. 원소 아홉 개를 손으로 적지 않아도 되기 때문이다.

이 회전은 방향만 맞춘다. 원점을 배의 위치에서 지구 중심으로 옮기려면 회전한 벡터에 배의 ECEF 위치를 더해야 하고, 그 위치를 위경도에서 구하는 식은 1.6에 있다.

1.5.2 ECI — Earth-Centered Inertial (지구중심 관성)

질문답
① 원점지구의 질량중심 (ECEF와 같다)
② 축특정 시각(보통 J2000.0)의 별자리 방향에 고정
③ 붙은 곳아무것도 아님 — 지구가 돌아도 축은 그대로 있다

ECEF와 원점은 같고 축이 도느냐 안 도느냐만 다르다.

ECEF 와 ECI 의 차이 그림 1-9. 원점은 같고 축이 도느냐 마느냐만 다르다.

1
2
ECEF :  지구와 함께 돈다   →  지상의 건물은 좌표가 안 변한다
ECI  :  돌지 않는다        →  지상의 건물은 좌표가 하루 한 바퀴 원을 그린다

왜 이런 게 필요한가. 물리 법칙 때문이다. 뉴턴의 운동법칙(F = ma)은 회전하지 않는 좌표계에서만 그대로 성립한다. 회전하는 좌표계에서 물체의 운동을 계산하려면 원심력·코리올리력 같은 실제로 존재하지 않는 가짜 힘을 추가로 넣어야 한다.

그래서 중력만 받고 날아가는 물체 — 인공위성, 탄도미사일 — 를 계산할 때는 ECI를 쓴다. ECI에서는 중력 하나만 넣으면 되지만, ECEF에서는 가짜 힘들을 전부 모델링해야 한다. (참고로 INS 내부 계산도 ECI를 기준으로 출발한다. 1.3의 i 프레임이 이것이다.)

ECEF ↔ ECI 변환. 지구 자전만 고려하면 z축 둘레 회전 하나로 끝난다.

1
2
p_ECI  = Rz(θ)  × p_ECEF        θ = 그리니치 항성시(GMST), 시각으로 계산
p_ECEF = Rz(θ)ᵀ × p_ECI

θ 는 그 시각에 지구가 돌아가 있는 각이다. ECEF 의 X축(그리니치 자오선 쪽)이 ECI 의 X축에서 동쪽으로 얼마나 돌아가 있는지를 잰 값이고, 그리니치 항성시(GMST, Greenwich Mean Sidereal Time)라고 부른다. 지구는 거의 일정한 빠르기로 돌기 때문에, 기준 시각의 각에 그 뒤로 지난 날수만큼 돈 각을 더하면 어림값이 나온다.

1
2
3
4
5
6
θ ≈ 280.46061837° + 360.98564736629° × (JD − 2451545.0)      ← 360° 로 나눈 나머지를 쓴다

  JD                = 율리우스일(Julian Date). 날짜와 시각을 하루 단위의 수 하나로 적은 것
  2451545.0         = 2000년 1월 1일 12시(세계시)의 율리우스일
  280.46061837°     = 그 시각의 θ
  360.98564736629°  = 지구가 하루 동안 도는 각

하루 동안 도는 각이 360° 보다 1° 가까이 큰 것은 하루를 해 기준으로 세기 때문이다. 지구가 한 바퀴 도는 사이 공전 궤도를 따라 조금 나아가 있어서, 해가 다시 같은 자리에 오려면 그만큼 더 돌아야 한다.

뒤에서 쓸 예제의 함선 위치(위도 36.408°, 경도 127.307°, 고도 0 m)가 2026년 9월 1일 12시(한국 시각)에 ECI 로는 어디인지 계산해 보겠다. 한국 시각 12시는 세계시 3시이고, 율리우스일로는 2461284.625 다. 아래의 함선 ECEF 좌표는 이 위도·경도·고도를 1.6의 LLA → ECEF 식에 넣어 구한 값이고, 거기에 Rz(θ) 를 곱한 것이 함선 ECI 좌표다.

1
2
3
4
5
JD − 2451545.0 = 9 739.625 일
θ = 25.296348272°                      ← 360° 를 9 767 번 덜어 내고 남은 각

함선 ECEF   (−3 114 830.1107,  +4 087 762.8525,  +3 764 723.0978) m
함선 ECI    (−4 562 850.4352,  +2 364 818.7377,  +3 764 723.0978) m

z 성분이 그대로다. z축 둘레로 돌렸기 때문이다. 원점에서의 거리도 그대로다. 회전은 길이를 바꾸지 않는다. θ 와 좌표의 자릿수를 길게 적은 것은 계산을 따라 해 볼 수 있게 하려는 것이고, 어림식이 그 자리까지 정확하다는 뜻은 아니다.

정밀도가 더 필요하면 세차·장동(지구 자전축이 천천히 흔들리는 현상)과 극운동까지 넣어야 하지만, 입문 단계에서는 “시각을 알면 각도가 나오고, 그 각도만큼 z축 둘레로 돌리면 된다“로 충분하다.

시각 정밀도 감각 지구는 적도에서 초속 약 465 m로 돈다. 즉 시각이 1 ms 어긋나면 적도에서 위치가 약 0.47 m 어긋난다. 여러 센서의 결과를 합칠 때 시각 동기가 왜 중요한지 보여주는 숫자다.

1.5.3 이 설계에 ECI가 필요한가

필요 없다. 여기서 다루는 표적은 항공기/수상함 수준이고, 결과를 지구 기준 위경도로 내면 된다. ECI는 표적을 관성 동역학으로 전파해야 할 때(탄도탄 방어, 위성 추적) 필요하다.


1.6 LLA 좌표계

우리가 아는 그 지도 좌표다. Latitude(위도) · Longitude(경도) · Altitude(고도).

질문답
① 원점(직교좌표가 아니라 각도 좌표라 원점 개념이 다르다) 기준은 적도면과 그리니치 자오선
② 축위도 = 적도에서 남북으로 잰 각 / 경도 = 그리니치에서 동서로 잰 각 / 고도 = 타원체 표면에서 잰 높이
③ 붙은 곳지구

지구는 구가 아니라 타원체다

정확한 위경도를 다루려면 지구를 어떤 모양으로 볼 것인가를 먼저 정해야 한다. 그 표준이 WGS-84이고, GPS와 해도가 쓰는 것도 이것이다.

상수값의미
a (장반경)6 378 137.0 m적도 방향 반지름
1/f (편평률의 역수)298.257223563얼마나 납작한가
b (단반경)6 356 752.314245 m극 방향 반지름 (계산으로 나옴)
e² (제1이심률 제곱)0.00669437999014타원 공식에 계속 등장

적도 반지름과 극 반지름의 차이가 약 21 km다. 지구 크기에 비하면 0.3 % 정도지만, 이걸 무시하고 구로 계산하면 위치 오차가 수십 km까지 벌어진다.

지구 타원체, 측지위도와 지심위도, 고도의 두 기준 그림 1-10. 타원체 법선은 지구 중심을 지나지 않는다. 그래서 측지위도와 지심위도가 다르고, 고도의 기준도 타원체(LLA)와 지오이드(해발)로 갈린다.

LLA ↔ ECEF 변환

LLA → ECEF (한 번에 계산됨)

1
2
3
4
N = a / √(1 − e²·sin²(위도))        ← 그 위도에서의 곡률반경
X = (N + 고도) × cos(위도) × cos(경도)
Y = (N + 고도) × cos(위도) × sin(경도)
Z = (N·(1 − e²) + 고도) × sin(위도)   ← (1 − e²) 빼먹는 실수가 가장 흔하다

ECEF → LLA (반복이 필요하다)

이 방향은 딱 떨어지는 공식이 없다. 위도를 알아야 N을 구하는데, N을 알아야 위도를 구할 수 있는 닭-달걀 관계이기 때문이다. 그래서 위도를 어림값으로 놓고 시작해, 고도와 위도를 번갈아 다시 계산하면서 값이 더 변하지 않을 때까지 되풀이한다. 경도는 이 문제와 무관해서 X 와 Y 만으로 바로 나온다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ρ    = √(X² + Y²)                      ← 지구 자전축에서 떨어진 거리
경도 = atan2(Y, X)                      ← 반복 없이 바로 나온다

초기값
  위도₀ = atan2(Z, ρ)                   ← 지구를 구로 보고 잡은 위도
  고도₀ = 0

반복
  N    = a / √(1 − e²·sin²(위도₀))
  고도 = ρ / cos(위도₀) − N
  위도 = atan2( Z·(N + 고도),  ρ·(N·(1 − e²) + 고도) )

  위도 변화가 1e−12 rad 보다 작고 고도 변화가 1e−6 m 보다 작으면 멈춘다.
  아니면 방금 구한 위도와 고도를 위도₀, 고도₀ 로 삼아 다시 계산한다. 최대 10회

반복 안의 고도 식과 위도 식은 바로 위 LLA → ECEF 식을 거꾸로 푼 것이다. X 식과 Y 식을 제곱해 더하고 제곱근을 잡으면 ρ = (N + 고도) × cos(위도) 이고, 이것을 고도에 대해 풀면 고도 식이 된다. Z 식을 이 ρ 식으로 나누어 정리하면 tan(위도) = (Z / ρ) × (N + 고도) / (N·(1 − e²) + 고도) 가 되고, 이것이 위도 식이다. 고도 식에는 cos(위도) 와 N 이, 위도 식에는 N 과 고도가 들어 있어서 두 식 모두 위도를 알아야 계산된다. 그래서 지금 가진 위도로 먼저 계산하고, 거기서 나온 위도를 다시 넣는다.

0.1절 표의 표적을 넣어 보겠다. ECEF 좌표는 (−3 132 053.7872, +4 078 527.1919, +3 760 545.9526) m 다.

회차위도 [°]고도 [m]
초기값36.1774113620
136.361399161−14 886.2737
236.36096615776.5738
336.36096717841.1992
436.36096717641.2826
536.36096717641.2824
636.36096717641.2824

6회차에서 위도와 고도의 변화가 둘 다 허용치 아래로 내려가 멈춘다. 위도는 4회차, 고도는 5회차에 이미 표에 적은 자리까지 값이 정해져 있고, 6회차는 그것을 확인하는 회차다. 결과는 위도 36.360967176°, 경도 127.522008441°, 고도 41.2824 m 다.

초기값 36.177° 는 수렴한 위도 36.361° 보다 0.18° 작다. 두 각을 재는 기준이 다르기 때문이다. atan2(Z, ρ) 는 지구 중심과 그 점을 이은 선이 적도면과 이루는 각(지심위도)이고, 구하려는 위도는 타원체 면에 세운 수직선이 적도면과 이루는 각(측지위도)이다. 그림 1-10 에 그린 두 각이다. 둘의 차이는 위도 45° 부근에서 가장 커서 0.192° 이고, 이 표적의 위도에서는 0.18° 다. 1회차 고도가 −14 886 m 로 크게 빗나간 것도 이 어긋난 위도로 고도를 먼저 계산했기 때문이다. 위도가 한 번 고쳐진 2회차부터는 수십 m 안으로 들어온다.

극점에서는 이 식을 쓸 수 없다. ρ 와 cos(위도) 가 함께 0 이 되어 고도 식의 ρ / cos(위도) 가 0/0 이 된다. 3절 코드는 이 순서를 그대로 옮긴 것이고 극점은 따로 처리하지 않았다.

고도의 함정 — 타원체고 vs 해발고도

LLA의 “고도”는 타원체 표면에서 잰 높이(타원체고, HAE)다. 그런데 해도나 고도계가 말하는 “고도”는 보통 평균해면 기준 높이(해발고도, MSL)다.

1
해발고도  =  타원체고  −  지오이드고

지오이드란 “바다가 육지 안까지 이어져 있다면 만들어질 평균 해수면”이다. 중력이 균일하지 않아서 이 면은 타원체와 일치하지 않고, 전 지구적으로 약 −107 m ~ +86 m 차이가 난다. 한국 근해에서는 대략 +20 ~ +25 m 다.

여기에 항공기가 쓰는 기압고도까지 끼어들면 기준이 셋이 된다. 시스템 내부에서는 한 가지로 통일하고, 표시할 때만 변환하는 것이 안전하다. 이번 설계는 내부적으로 타원체고를 쓴다.


1.7 한눈에 보는 정리

#좌표계원점축붙은 곳주 용도
1.1안테나안테나 면 중심x=왼쪽, y=위, z=보어사이트배표적 측정 (R, Az, El) / 빔 조향 (u, v)
1.2동체(Body)배 무게중심x=선수, y=우현, z=아래배탑재 장비의 공통 기준
1.3INSINS 설치 위치자이로·가속도계 3축배위치·자세·속도 공급
1.4NED / ENU배의 현재 위치N,E,D / E,N,U바다 위 지점 (회전 안 함)데이터 처리(추적·필터링)
1.5ECEF지구 중심Z=북극, X=적도∩그리니치지구 (자전함)지구 기준 공통 좌표, 위경도로 가는 다리
1.5ECI지구 중심별자리에 고정없음 (회전 안 함)위성·탄도 표적의 운동 계산
1.6LLA(적도면 / 그리니치)위도·경도·고도지구지도 표시, 외부 전달

변환 체인 한 장 요약

좌표변환 체인 전체도 그림 1-11. 각 단계의 회전과 평행이동을 무엇이 공급하는지가 오차 예산의 항목이 된다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 신호처리
    │  R, Az, El  (안테나 극좌표)
    ▼
[1.1] 안테나 좌표계  ── 극좌표 → 직교좌표
    │  회전: C_body_from_antenna = Rz(설치 방위) · Ry(백틸트) · A   (설치 도면, 상수. A 는 축 맞춤)
    │  이동: 레버암
    ▼
[1.2] 동체 좌표계
    │  회전: C_ned_from_body (INS 자세 — roll, pitch, yaw)     ← [1.3] INS가 공급
    ▼
[1.4] 로컬(NED) 좌표계  ── ★ 데이터 처리를 여기서 수행
    │  회전: C_ecef_from_ned (INS 위치 — 위도, 경도)           ← [1.3] INS가 공급
    │  이동: 배의 ECEF 위치
    ▼
[1.5] ECEF 좌표계
    │  위도·고도 반복 계산
    ▼
[1.6] LLA 좌표계  ── 전시기·통제제어로 전달

2. 실전 설계 — 함선 고정형 레이다

좌표계를 하나씩 봤으니 이제 실제 시스템을 하나 놓고 설계해 보자. 신호처리가 넘겨준 측정값을 받아 데이터 처리를 거쳐 운용 전시기까지 보내는 경로다.

2.1 상황

  • 플랫폼 — 해상에서 운용하는 함선. 레이다는 고정형 안테나(기계적으로 돌지 않는다)
  • 안테나 설치 — 선수 기준 시계방향 90도를 바라보고, 무게중심에서 30 m 뒤 / 10 m 위
  • INS 설치 — 플랫폼 무게중심과 축이 일치
  • 입력 — 신호처리가 주는 값은 안테나 기준 거리·방위각·고각
  • 처리 — 데이터 처리는 로컬 좌표계에서 수행
  • 출력 — 결과를 위/경/고도와 안테나 극좌표(극좌표) 두 가지로 통제제어에 전달

2.1.1 각 조건이 설계를 어떻게 묶는가

조건내용설계에 미치는 영향
① 플랫폼해상 함선, 고정형 안테나안테나가 기계적으로 돌지 않는다 → 안테나↔동체 회전행렬이 상수. 대신 배가 흔들리는 것을 전부 계산으로 보상해야 한다
② 안테나 설치선수 기준 시계방향 90°, 무게중심에서 30 m 뒤 / 10 m 위설치 방위 = 90°, 백틸트 = 0°, 레버암 = (−30, 0, −10) m
③ INS 설치무게중심과 축 일치INS 좌표계 = 동체 좌표계 → 변환 단계 하나가 통째로 생략된다
④ 입력안테나 기준 거리·방위각·고각체인의 출발점이 안테나 극좌표로 확정
⑤ 처리로컬 좌표계에서 수행체인의 중간 목적지가 NED로 확정
⑥ 출력위/경/고도 + 안테나 극좌표정변환뿐 아니라 역변환도 필요하다

즉 ④·⑤·⑥ 이 체인의 출발점·중간점·도착점을 이미 지정하고 있다. 설계자가 결정할 것은 그 사이를 어떤 좌표계로 어떻게 이을 것인가다.

2.1.2 전체 데이터 흐름

1
2
3
4
5
6
7
8
9
┌──────────┐   R, Az, El     ┌────────────────────────┐   위/경/고도   ┌────────────┐
│  신호처리  │ ─────────────> │       데이터 처리        │ ────────────> │  통제제어   │
│          │  (안테나 극좌표)  │  (좌표변환 + 추적·필터링)  │   R, Az, El   │   전시기    │
└──────────┘                 └────────────────────────┘ ────────────> └────────────┘
                                     ▲
                                     │  위도·경도·고도, roll·pitch·yaw
                                ┌─────────┐
                                │   INS    │
                                └─────────┘

데이터 처리 내부는 다시 세 부분으로 나뉜다.

1
2
3
4
 [입력 변환]  안테나 극좌표 ──> 로컬(NED)        ... 2.1.4의 1~3단계
 [처리]       로컬 좌표계에서 추적·필터링         ... ⑤
 [출력 변환]  로컬 ──> 위/경/고도  (정변환)       ... 2.1.4의 4~5단계
              로컬 ──> 안테나 극좌표 (역변환)     ... 2.1.5

2.1.3 단계별 좌표계 선정과 그 이유

단계선정한 좌표계표현선정 이유
0안테나 좌표계 — 극좌표(Polar)R, Az, El④. 레이다가 물리적으로 잴 수 있는 것이 거리와 두 각도뿐이므로 원측정값은 필연적으로 극좌표다
1안테나 좌표계 — 직교좌표(Cartesian)x, y, z (왼쪽, 위, 보어사이트)회전행렬을 곱하려면 벡터 형태여야 한다. 극좌표끼리는 회전을 합성할 수 없다
2동체(Body) 좌표계x, y, z (선수, 우현, 아래)안테나·INS·무장이 만나는 유일한 공통 기준. 안테나는 배에 붙어 있으므로 배를 거치지 않고 바깥 세상으로 나갈 방법이 없다. 또 여러 안테나·센서를 붙일 때 이 단계가 확장점이 된다
3로컬(NED) 좌표계N, E, D⑤. 배가 흔들려도 이 좌표계는 흔들리지 않는다. 표적이 실제로 움직인 것과 배가 흔들린 것을 구분할 수 있어야 추적 필터가 성립한다. ENU가 아니라 NED를 고른 것은 동체 FRD와 축 의미가 맞아 자세 회전이 자연스럽게 연결되기 때문
4ECEF 좌표계X, Y, Z로컬 좌표계에서 위경도로 곧장 갈 수 없다. 로컬 원점(배 위치)이 지구 어디인지는 지구 기준 좌표로만 표현되므로, 지구중심 직교좌표를 다리로 삼는다
5LLA 좌표계위도, 경도, 고도⑥. 전시기·통제제어가 지도 위에 표시하는 형식
역안테나 좌표계 — 극좌표R, Az, El⑥의 두 번째 출력. 처리 결과를 안테나 기준 각으로 되돌려야 다음 빔을 보낼 방향을 정할 수 있고 운용자도 안테나 기준으로 확인할 수 있다
선정하지 않은 좌표계와 그 이유
좌표계왜 쓰지 않는가
INS 좌표계 (별도 단계로)③에 의해 동체 좌표계와 일치. 별도 변환 단계를 두면 항등행렬을 곱하는 셈이라 계산만 늘고 오류 가능성만 커진다. INS는 좌표계가 아니라 값의 공급원으로 설계에 반영한다
UV 좌표계빔 조향·빔 형성은 안테나 내부(신호처리 이전)의 일이다. 또 UV는 방향만 있고 거리가 없어 표적 위치를 옮기는 이 체인에는 넣을 값이 없다
ECI 좌표계표적을 관성 동역학으로 전파할 때 필요하다(위성·탄도탄). 이 설계는 지구 기준 위치 산출이 목적이므로 불필요
ENU 좌표계NED와 정보량이 같다. 둘 다 쓰면 축 순서·부호 혼동만 생기므로 하나로 통일한다

2.1.4 단계별 좌표변환 과정

각 단계에서 무엇이 회전을 공급하고, 무엇이 평행이동을 공급하는지를 함께 적는다. 식 아래에는 예제값이 그 단계에서 어떤 값이 되는지도 붙였다. 예제의 측정값은 R = 20 km, Az = 30°, El = 0° 이고, 배는 위도 36.408°, 경도 127.307°, 고도 0 m 에 선수를 진북에서 시계방향 45° 로 두고 기울지 않은 채 멈춰 있다 (roll = pitch = 0°, yaw = 45°). 플랫폼이 이동하지 않고 고정되어 있으므로 배의 위치와 자세는 계산하는 동안 바뀌지 않는다.


1단계 — 안테나 극좌표 → 안테나 직교좌표

변환 종류: 표현 변경 (회전·이동 없음)

1
2
3
x_a = R × cos(El) × sin(Az)      ← 왼쪽 성분
y_a = R × sin(El)                ← 위 성분
z_a = R × cos(El) × cos(Az)      ← 보어사이트 성분

이유: 회전행렬은 벡터에만 곱할 수 있다. 각도로 표현된 방향을 세 성분으로 풀어쓰는 준비 단계다.

예제값을 넣으면 (x_a, y_a, z_a) = (10 000.0000, 0.0000, 17 320.5081) m 다. 고각이 0° 라 위 성분은 없고, 20 km 가 보어사이트 쪽 약 17.3 km 와 왼쪽 10 km 로 나뉜다.


2단계 — 안테나 좌표계 → 동체 좌표계

변환 종류: 회전 + 평행이동

1
2
3
4
5
6
7
8
9
10
C_body_from_antenna = Rz(ψ_mount) · Ry(θ_mount) · A
p_body = C_body_from_antenna × p_antenna + r_lever

  A       = 축 맞춤 행렬   (각도가 들어가지 않는 상수)
            A × (x, y, z) = (z, −x, −y) = (정면, 오른쪽, 아래)
  ψ_mount = 90°           (조건 ②: 선수 기준 시계방향 90도)
  θ_mount = 0°            (조건 ②에 백틸트 언급 없음)
  r_lever = (−30, 0, −10) m
            └ x = −30 : 30 m 뒤 (동체 x축은 앞이 +)
              z = −10 : 10 m 위 (동체 z축은 아래가 +)

세 행렬은 벡터에 가까운 오른쪽 것부터 차례로 곱해진다. 축 맞춤 행렬 A 가 가장 먼저이고 그다음이 백틸트 Ry, 설치 방위 Rz 이며, 마지막에 레버암을 더한다. A 를 가장 먼저 곱해야 하는 까닭은 1.2의 「축 맞춤 행렬 A」에 있다.

회전의 공급원: 설치 도면 (고정형이므로 상수, 조건 ①·②) 평행이동의 공급원: 설치 도면 (레버암, 조건 ②)

예제의 설치값을 넣으면 세 행렬의 곱은 상수 행렬 하나가 된다.

1
2
3
4
5
6
7
8
                                             ⎡ 1   0   0 ⎤
C_body_from_antenna = Rz(90°) · Ry(0°) · A = ⎢ 0   0   1 ⎥
                                             ⎣ 0  −1   0 ⎦

p_antenna                       = (10 000.0000,      0.0000, 17 320.5081) m     (왼쪽, 위, 보어사이트)
C_body_from_antenna × p_antenna = (10 000.0000, 17 320.5081,      0.0000) m     (선수, 우현, 아래)
r_lever                         = (   −30.0000,      0.0000,    −10.0000) m     (선수, 우현, 아래)
p_body                          = ( 9 970.0000, 17 320.5081,    −10.0000) m     (위 두 줄의 합)

안테나가 우현을 보고 있어서 안테나의 왼쪽(x)이 배의 선수 쪽이고 보어사이트(z)가 배의 우현 쪽이다. 그래서 정면에서 왼쪽으로 30° 에 있는 표적의 선수 성분이 +10 000 m 로 나온다. 여기에 레버암을 더하면 원점이 안테나에서 무게중심으로 옮겨 간다. 값이 줄마다 왜 이렇게 바뀌는지는 1.2에 풀어 두었다.

왜 레버암을 반영해야 하는가: 20 km 표적에 30 m는 무시해도 될 것 같지만, 이 오차는 거리가 멀어져도 줄지 않는다. 근거리 표적이나 정밀 사격통제에서 그대로 문제가 되고, 안테나가 여러 개일 때는 안테나마다 다른 레버암 때문에 같은 표적의 트랙이 갈라진다. (3.2절에서 30 m 어긋남을 실제로 확인한다.)


3단계 — 동체 좌표계 → 로컬(NED) 좌표계

변환 종류: 회전만 (원점을 배 위치에 두므로 이동 없음)

1
2
C_ned_from_body = Rz(yaw) · Ry(pitch) · Rx(roll)
p_ned = C_ned_from_body × p_body

회전의 공급원: INS의 자세각 (조건 ③에 의해 INS 좌표계 = 동체 좌표계)

이 단계가 설계의 핵심인 이유: 고정형 안테나(①)이기 때문에, 배의 흔들림을 보상할 기계적 수단이 없다. 전부 이 회전행렬이 감당한다. 따라서 이 단계의 정확도가 시스템 전체의 지향 정확도를 결정한다.

예제에서는 roll 과 pitch 가 0° 라 yaw 45° 회전 하나만 남는다. z축 둘레 회전이므로 선수·우현 성분만 섞이고 아래 성분은 그대로다.

1
2
3
4
5
C_ned_from_body = Rz(45°) · Ry(0°) · Rx(0°) = Rz(45°)

N = cos45° × 9 970.0000 − sin45° × 17 320.5081 =  −5 197.5941 m
E = sin45° × 9 970.0000 + cos45° × 17 320.5081 = +19 297.3033 m
D =                                                  −10.0000 m

배가 있는 지점에서 남쪽으로 5 198 m, 동쪽으로 19 297 m, 위로 10 m 라는 뜻이다. 수평 성분을 진북에서 시계방향으로 잰 방위로 바꿔 읽으면 atan2(E, N) = 105.07° 다.

예제 배치도. 진북, 선수, 안테나 정면, 표적 사이의 각도 그림 2-1. 배가 기울지 않았을 때는 진북에서 선수까지 45°, 선수에서 안테나 정면까지 90° 를 더하고, 표적이 정면에서 선수 쪽으로 되돌아와 있는 30° 를 뺀다: 45° + 90° − 30° = 105°.

왜 각도 덧셈으로는 안 되는가 배가 기울지 않았다면 진북 기준 방위 = yaw + 설치방위 − Az 로 간단히 계산된다 (45° + 90° − 30° = 105°). Az 를 빼는 것은 방위각을 왼쪽으로 + 로 재기 때문이다. 진북 기준 방위는 위에서 내려다볼 때 시계방향으로 재는데, 안테나 정면에서 왼쪽은 반시계방향이다. 위에서 구한 105.07° 와의 차이 0.07° 는 레버암에서 나온다. 그러나 roll이나 pitch가 조금이라도 있으면 이 덧셈은 즉시 깨진다. 세 축의 회전이 서로 섞이기 때문이다. 3.2절에서 roll 20°만으로 고각이 0°에서 −17°로 바뀌는 것을 확인한다. 해상 플랫폼에서 회전행렬은 선택이 아니라 필수다.

이 단계에서 반드시 챙겨야 할 것 — 자세 데이터의 시각

회전행렬을 정확히 만들어도, 언제 측정한 자세인가가 틀리면 소용이 없다. 배는 파도에 계속 흔들리고 있으므로 자세는 매 순간 변한다.

1
2
3
4
5
함정 횡동요 각속도 10°/s 일 때 지연에 따른 지향 오차

    지연  2 ms  →  0.020°  (0.35 mrad)
    지연 10 ms  →  0.100°  (1.75 mrad)     ← 20 km 에서 약 35 m
    지연 50 ms  →  0.500°  (8.73 mrad)     ← 20 km 에서 약 175 m

이 값은 대개 INS 자체의 자세 정확도보다 크다. 즉 좋은 INS를 쓰고도 시각 처리를 잘못하면 그 성능을 그대로 버리는 셈이다. 설계에서 챙길 것은 세 가지다.

  1. INS가 보내는 자세에 유효시각(언제 측정한 값인가) 타임스탬프를 함께 받는다.
  2. 각속도를 이용해 실제 측정이 일어난 시각으로 자세를 외삽한다.
  3. 레이다와 INS를 같은 시각원(PTP/IRIG-B 등)에 물려 시계를 맞춘다.

(이 글의 예제는 플랫폼이 고정되어 있어 이 문제가 드러나지 않지만, 실제 해상 운용에서는 3단계의 정확도를 좌우하는 가장 큰 요인이다.)


【 데이터 처리 구간 】 조건 ⑤

로컬(NED) 좌표계에서 추적·필터링을 수행한다. 이 좌표계를 고른 이유는 위에서 설명한 대로 배의 흔들림과 무관하기 때문이다. 안테나 좌표계에서 추적하면 배가 롤링할 때마다 표적이 좌우로 튀는 것처럼 보여 필터가 발산한다.

(이 글은 변환 검증이 목적이므로 필터를 두지 않고 값을 그대로 통과시킨다. 그 덕분에 역변환 결과가 입력과 정확히 같아야 한다는 강력한 검증 조건이 생긴다.)


4단계 — 로컬(NED) 좌표계 → ECEF 좌표계

변환 종류: 회전 + 평행이동

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
① 배의 위치를 ECEF로 변환 (LLA → ECEF 공식)
   N = a / √(1 − e²·sin²(lat))
   p_ecef_platform = ( (N+h)·cos(lat)·cos(lon),
                       (N+h)·cos(lat)·sin(lon),
                       (N·(1−e²)+h)·sin(lat) )

② 로컬 벡터를 지구 방향으로 회전
                     ⎡ −sin(lat)·cos(lon)  −sin(lat)·sin(lon)   cos(lat) ⎤
   C_ned_from_ecef = ⎢       −sin(lon)            cos(lon)          0     ⎥
                     ⎣ −cos(lat)·cos(lon)  −cos(lat)·sin(lon)  −sin(lat) ⎦

   C_ecef_from_ned = (C_ned_from_ecef)ᵀ

③ 더한다
   p_ecef_target = p_ecef_platform + C_ecef_from_ned × p_ned

회전의 공급원: INS의 위치(위도·경도). 지구 위 어디에 있느냐에 따라 “북쪽”이 지구 기준으로 어느 방향인지가 완전히 달라지기 때문이다. 평행이동의 공급원: 배의 ECEF 위치

①의 N 은 1.6의 곡률반경이다. 3단계에서 북쪽 좌표로 쓴 N 과 글자만 같다. ②에 풀어 쓴 C_ned_from_ecef 의 세 행은 배가 있는 지점의 북·동·아래 방향을 ECEF 성분으로 적은 것이다. 이 세 방향이 위도와 경도에서 어떻게 정해지는지는 1.5.1에 있다. 거기서 φ, λ 로 적은 위도와 경도가 여기의 lat, lon 이다.

예제값으로 따라가면, 배의 위치(위도 36.408°, 경도 127.307°, 고도 0 m)를 ①에 넣어 함선의 ECEF 위치를 얻고, ②에서 3단계의 p_ned 를 지구 축 기준으로 돌린 다음, ③에서 둘을 더한다.

1
2
3
① p_ecef_platform         = (−3 114 830.1107, +4 087 762.8525, +3 764 723.0978) m
② C_ecef_from_ned × p_ned = (   −17 223.6765,     −9 235.6606,     −4 177.1452) m
③ p_ecef_target           = (−3 132 053.7872, +4 078 527.1919, +3 760 545.9526) m

②의 값은 p_ned 와 같은 벡터를 지구 축 X, Y, Z 로 다시 읽은 것이라 길이는 그대로다. 여기에 지구 중심에서 잰 배의 위치를 더하면 표적의 위치가 지구 중심 기준으로 나온다.


5단계 — ECEF 좌표계 → LLA 좌표계

변환 종류: 좌표계 종류 변경 (직교 → 각도)

경도는 X 와 Y 만으로 바로 나온다. 위도와 고도는 한 번에 풀리지 않는다. 위도를 알아야 곡률반경 N 이 나오고, N 을 알아야 위도가 나오기 때문이다. 그래서 지구를 구로 본 위도에서 출발해 N, 고도, 위도를 차례로 다시 계산하고, 값이 더 바뀌지 않을 때까지 되풀이한다. 회차마다 값이 어떻게 수렴하는지는 1.6에 표로 있다.

1
2
3
4
5
6
7
8
9
10
ρ    = √(X² + Y²)
lon  = atan2(Y, X)

초기값   lat₀ = atan2(Z, ρ)                           ← 지구를 구로 본 위도
         h₀   = 0
반복     N    = a / √(1 − e²·sin²(lat₀))
         h    = ρ / cos(lat₀) − N
         lat  = atan2( Z·(N + h),  ρ·(N·(1 − e²) + h) )
         lat 과 h 의 변화가 허용치 아래면 멈춘다.
         아니면 lat₀·h₀ 를 새 값으로 바꿔 다시 계산한다

예제 표적은 여섯 번째 계산에서 멈춘다.

1
lat = 36.360967176°,   lon = 127.522008441°,   h = 41.2824 m

3단계에서 위로 10 m 였던 표적이 여기서는 고도 41.3 m 로 나온다. 로컬 좌표계의 수평면은 평평한데 지구 표면은 20 km 앞에서 그 아래로 31.3 m 휘어 내려가 있어서, 그만큼이 더해진 것이다(1.4).

주의점: 여기서 나온 고도는 타원체고다. 전시기가 해발고도를 요구하면 지오이드고를 빼서 변환해야 한다(1.6 참조). 시스템 내부는 타원체고로 통일한다.

2.1.5 출력 2종과 역변환

조건 ⑥ 은 결과를 두 가지 형태로 요구한다.

출력얻는 방법
위/경/고도위의 정변환 1→5단계 결과 그대로
안테나 극좌표처리 결과를 역변환해서 안테나 기준으로 되돌린다

안테나 극좌표를 따로 내보내는 까닭은, 처리한 결과를 다시 안테나 기준 각으로 되돌려야 다음 빔을 어디로 보낼지 정할 수 있기 때문이다.

역변환은 새로 배울 것이 없다. 세 가지만 지키면 된다.

1
2
3
① 회전은 전치(transpose)     C_ned_from_body  →  (C_ned_from_body)ᵀ
② 평행이동은 빼기            + 레버암          →  − 레버암
③ 순서는 거꾸로
1
2
3
4
5
6
7
8
표적 LLA
  → ECEF                                    (정변환과 같은 공식)
  → 로컬 NED    C_ned_from_ecef × (표적_ecef − 배_ecef)
  → 동체        (C_ned_from_body)ᵀ × p_ned
  → 안테나      (C_body_from_antenna)ᵀ × (p_body − 레버암)
  → 극좌표      R  = √(x² + y² + z²)
                Az = atan2(x, z)                ← 보어사이트(z)에서 왼쪽(x)으로
                El = atan2(y, √(x² + z²))       ← 안테나 수평면에서 위(y)로

C_body_from_antenna 안에 축 맞춤 행렬 A 가 들어 있으므로, 전치 하나로 설치 방위·백틸트와 함께 축도 안테나 면 기준(x 왼쪽, y 위, z 보어사이트)으로 되돌아온다. 마지막 극좌표 식은 1단계 식을 거꾸로 푼 것이다.

위 체인은 정변환을 끝에서부터 되짚도록 표적의 위/경/고도에서 시작해 적었다. 2.1.2의 흐름대로 로컬(NED)에 있는 처리 결과를 안테나 극좌표로 내보낼 때는 p_ned 가 이미 있으므로 → ECEF 와 → 로컬 NED 두 줄을 건너뛰고 → 동체 줄부터 계산한다. 3절에서는 다섯 단계를 모두 확인하려고 위/경/고도에서 출발한다.

역변환이 곧 검증 수단이 된다 필터를 거치지 않았다면 역변환 결과는 입력값과 정확히 같아야 한다. 정변환과 역변환 가운데 한쪽에만 부호나 전치 방향이 틀린 곳이 있으면 왕복한 값이 입력과 어긋난다. 3.1.6절의 왕복 확인에서 실제로 이것을 돌려 보고, 왕복만으로는 드러나지 않는 실수가 무엇인지도 거기에 적었다. (실제 운용에서는 필터가 값을 다듬으므로 필터가 보정한 만큼만 달라진다. 그 차이가 예상 범위를 벗어나면 필터나 변환 중 하나가 이상하다는 신호다.)

2.1.6 설계 결정 요약

결정 항목선택근거
데이터 처리 좌표계NED (로컬)조건 ⑤ + 배의 흔들림과 무관
로컬 좌표계 종류NED (ENU 아님)동체 FRD와 축 의미 일치 → 자세 회전이 자연스럽게 연결
로컬 원점배의 현재 위치배가 고정되어 있으므로 고정 원점과 동일. 이동 플랫폼이면 원점 갱신 정책 필요
z축 방향동체·NED 는 아래(Down) 로 통일. 안테나는 면 기준 축을 쓰고 축 맞춤 행렬 한 곳에서만 바꾼다섞이면 고도 부호가 뒤집힌다. 축을 바꾸는 자리를 하나로 모아 두면 실수를 찾기 쉽다
오일러각 순서3-2-1 (yaw→pitch→roll)항법 분야 표준
고도 기준타원체고(HAE) 내부 통일해발고도는 표시 단계에서만 변환
INS 좌표계 처리별도 단계 없음조건 ③ (무게중심·축 일치)
레버암반영거리와 무관한 계통 오차이므로
각도 역계산atan2 사용 (asin 아님)NaN 방지 + 극단 각도에서의 정밀도
검증 방법왕복 확인 (roll·pitch 가 0 인 대칭 자세와 세 각이 서로 다른 비대칭 자세 모두) 과 손으로 계산한 값 대조왕복은 정변환과 역변환이 서로 어긋난 부호·전치 오류를 찾고, 양쪽이 함께 틀린 경우는 손으로 계산한 값(60°, 105°)과 맞춰 찾는다

3. C 로 짜서 확인하기

3.1 구현과 실행

3.1.1 파일 구성과 빌드

1
2
3
4
code/
├── coord_frames.h    구조체·상수·함수 선언        (112줄)
├── coord_frames.c    변환 구현                    (300줄)
└── main.c            예제 입력값 실행 프로그램    (106줄)

Visual Studio 2022(MSVC)에서 경고 수준 /W4 로 빌드했고 경고 없이 빌드된다. 소스 파일을 UTF-8 로 저장해서 컴파일 옵션에 /utf-8 을 함께 주었다. 변환 코드가 쓰는 라이브러리는 표준 C 의 math.h 하나다.

3.1.2 구현 원칙 세 가지

좌표변환 버그의 특징은 컴파일도 되고 값도 그럴듯해 보인다는 것이다. 그래서 애초에 실수하기 어렵게 만드는 것이 중요하다.

① 함수 이름에 변환 방향을 넣는다

어느 좌표에서 어느 좌표로 가는 변환인지를 함수 이름에 적었다.

1
2
3
// 3) 동체 <-> NED (x 북, y 동, z 아래)
ST_Vec3     f_BodyToNed(const ST_Attitude stAtt, const ST_Vec3 stBody);
ST_Vec3     f_NedToBody(const ST_Attitude stAtt, const ST_Vec3 stNed);

f_BodyToNed 는 동체 좌표를 받아 NED 좌표를 돌려주고 f_NedToBody 는 그 반대다. 이름에 방향이 있어서 호출하는 쪽에서 헷갈릴 일이 없다. 변환 함수는 2.1.4의 다섯 단계마다 이런 쌍이 하나씩 있어 모두 열 개다.

인자로 오가는 구조체는 다섯 가지다. ST_Vec3 는 x·y·z 세 성분이고, ST_Polar 는 거리·방위각·고각, ST_Lla 는 위도·경도·고도, ST_Attitude 는 roll·pitch·yaw, ST_Mount 는 설치 방위·백틸트·레버암을 묶은 것이다.

② 회전행렬은 기본 회전 셋의 곱으로만 만든다

원소를 손으로 적는 행렬은 x·y·z 축 둘레의 기본 회전 f_RotX, f_RotY, f_RotZ 셋뿐이다. 안테나에서 동체로, 동체에서 NED 로, NED 에서 ECEF 로 가는 회전행렬은 전부 이 셋을 곱해서 만든다. 아홉 칸을 삼각함수로 풀어 적는 곳이 없으니 부호를 잘못 옮길 자리가 기본 회전 셋으로 줄어든다. 되돌아갈 때는 같은 행렬을 전치해서 쓰고, 역방향 행렬을 따로 만들지 않는다.

③ 단위는 안에서 라디안과 미터로 통일한다

구조체에 담는 각도는 전부 라디안이고 길이는 전부 미터다. 도(degree)는 main.c 에서 입력값을 넣을 때 DEG2RAD 로, 화면에 찍을 때 RAD2DEG 로 한 번씩만 바꾼다. 라디안과 도를 섞는 실수는 좌표변환에서 가장 흔한 버그다.

3.1.3 공식과 코드

2.1.4와 2.1.5의 식을 옮긴 함수를 단계별로 식과 나란히 놓았다. 회전행렬은 3×3 배열 e[3][3] 하나를 담은 구조체 ST_Mat3 로 두었다. 행렬에 하는 연산은 3×3 곱(f_MatMul), 전치(f_MatT), 행렬과 벡터의 곱(f_MatVec) 세 가지뿐이고 역행렬을 구하는 함수는 없다. f_VecAdd 와 f_VecSub 는 성분별 덧셈과 뺄셈이고, FLOAT64·INT32·VOID 는 double·int·void 에 붙인 이름이다.

회전행렬
1
2
3
         ⎡ cos γ  −sin γ   0 ⎤
Rz(γ) =  ⎢ sin γ   cos γ   0 ⎥
         ⎣   0       0     1 ⎦
1
2
3
4
5
6
7
8
9
10
11
static ST_Mat3 f_RotZ(const FLOAT64 dAng)
{
    ST_Mat3 stR = f_MatZero();

    stR.e[0][0] =  cos(dAng);
    stR.e[0][1] = -sin(dAng);
    stR.e[1][0] =  sin(dAng);
    stR.e[1][1] =  cos(dAng);
    stR.e[2][2] = 1.0;
    return stR;
}

0.3절의 Rz(γ) 를 배열에 칸 그대로 옮겼다. 1행 1열이 e[0][0], 1행 2열이 e[0][1] 이다. 전체를 0 으로 채운 뒤 0 이 아닌 다섯 칸만 넣는다. f_RotX 와 f_RotY 도 같은 방식이다.

1단계 — 안테나 극좌표 ↔ 직교좌표
1
2
3
x = R × cos(El) × sin(Az)          R  = √(x² + y² + z²)
y = R × sin(El)                    Az = atan2(x, z)
z = R × cos(El) × cos(Az)          El = atan2(y, √(x² + z²))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
ST_Vec3 f_PolarToXyz(const ST_Polar stPolar)
{
    ST_Vec3 stXyz;

    stXyz.dX = stPolar.dRange * cos(stPolar.dEl) * sin(stPolar.dAz);   // x 는 왼쪽. Az 가 반시계(+) 라 sin
    stXyz.dY = stPolar.dRange * sin(stPolar.dEl);                      // y 는 위
    stXyz.dZ = stPolar.dRange * cos(stPolar.dEl) * cos(stPolar.dAz);   // z 는 보어사이트
    return stXyz;
}

ST_Polar f_XyzToPolar(const ST_Vec3 stXyz)
{
    ST_Polar stPolar;
    FLOAT64  dXz = sqrt(stXyz.dX * stXyz.dX + stXyz.dZ * stXyz.dZ);    // 안테나 면의 좌우-보어사이트 평면 거리

    stPolar.dRange = sqrt(dXz * dXz + stXyz.dY * stXyz.dY);
    stPolar.dAz    = atan2(stXyz.dX, stXyz.dZ);                        // 보어사이트(z) 에서 왼쪽(x) 으로
    stPolar.dEl    = atan2(stXyz.dY, dXz);                             // 그 평면에서 위(y) 로
    return stPolar;
}

역변환의 dXz 는 √(x² + z²) 로, 표적까지의 벡터를 x축(좌우)과 z축(보어사이트)이 만드는 평면에 내린 길이다. 방위각은 그 평면 안에서 보어사이트(z)로부터 왼쪽(x)으로 잰 각이고, 고각은 그 평면에서 위(y)로 올라간 각이다. 두 각 모두 1.1.2절에서 본 대로 asin 이 아니라 atan2 로 되찾는다.

2단계 — 안테나 ↔ 동체
1
2
3
C_body_from_antenna = Rz(설치 방위) · Ry(백틸트) · A          A = Rz(−90°) · Rx(−90°)
p_body    = C_body_from_antenna × p_antenna + 레버암
p_antenna = (C_body_from_antenna)ᵀ × (p_body − 레버암)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
// 안테나 면 기준 축(x 왼쪽, y 위, z 보어사이트) 을 동체 축(x 선수, y 우현, z 아래) 으로 돌려 놓는 상수 행렬.
// A = Rz(-90) * Rx(-90)
static ST_Mat3 f_RotAntAxisToBody(VOID)
{
    return f_MatMul(f_RotZ(-PI / 2.0), f_RotX(-PI / 2.0));
}

// 안테나 -> 동체 회전. 축을 먼저 맞추고 설치 방위, 백틸트 순으로 돌림. 순서를 바꾸면 값이 달라짐
static ST_Mat3 f_RotAntToBody(const ST_Mount stMount)
{
    ST_Mat3 stMnt = f_MatMul(f_RotZ(stMount.dAz), f_RotY(stMount.dTilt));

    return f_MatMul(stMnt, f_RotAntAxisToBody());
}

// 회전하고 나서 설치 위치만큼 옮김
ST_Vec3 f_AntToBody(const ST_Mount stMount, const ST_Vec3 stAnt)
{
    ST_Vec3 stRot = f_MatVec(f_RotAntToBody(stMount), stAnt);

    return f_VecAdd(stRot, stMount.stOffset);
}

// 정변환 반대로. 먼저 빼고 전치한 행렬로 회전
ST_Vec3 f_BodyToAnt(const ST_Mount stMount, const ST_Vec3 stBody)
{
    ST_Vec3 stDiff = f_VecSub(stBody, stMount.stOffset);

    return f_MatVec(f_MatT(f_RotAntToBody(stMount)), stDiff);
}

A 는 1.2절의 축 맞춤 행렬이다. 식에서 맨 오른쪽에 있으므로 벡터에 가장 먼저 곱해지고, 코드에서도 f_MatMul(stMnt, f_RotAntAxisToBody()) 로 맨 오른쪽에 놓았다. f_RotAntAxisToBody 의 주석은 A 가 안테나 축을 동체 축(x 선수, y 우현, z 아래)으로 돌려 놓는다고 적었는데, 설치 방위와 백틸트가 0° 일 때 그렇다. 예제처럼 설치 방위가 있으면 A 를 곱한 값은 안테나의 정면·오른쪽·아래이고, 백틸트와 설치 방위까지 곱해야 배의 선수·우현·아래가 된다. stMount.dAz 가 설치 방위, stMount.dTilt 가 백틸트, stMount.stOffset 이 레버암이다. 주석의 “설치 방위, 백틸트 순”은 행렬을 왼쪽부터 곱해 적는 순서다. 벡터에는 오른쪽 행렬부터 곱해지므로 축 맞춤, 백틸트, 설치 방위 차례가 된다.

역변환은 더하기를 빼기로, 행렬을 전치로 바꾸고 순서도 뒤집는다. 레버암은 동체 좌표로 적은 값이라서 동체 좌표일 때 먼저 빼고, 그다음에 전치한 행렬로 안테나 축으로 되돌린다. 두 줄의 순서를 바꾸면 레버암을 엉뚱한 좌표계에서 빼게 된다.

3단계 — 동체 ↔ NED
1
2
C_ned_from_body = Rz(yaw) · Ry(pitch) · Rx(roll)
p_ned = C_ned_from_body × p_body          p_body = (C_ned_from_body)ᵀ × p_ned
1
2
3
4
5
6
7
// 동체 -> NED 회전. yaw, pitch, roll 순서로 곱함. 순서를 바꾸면 값이 달라짐
static ST_Mat3 f_RotBodyToNed(const ST_Attitude stAtt)
{
    ST_Mat3 stZy = f_MatMul(f_RotZ(stAtt.dYaw), f_RotY(stAtt.dPitch));

    return f_MatMul(stZy, f_RotX(stAtt.dRoll));
}

두 좌표계는 원점이 같아서 더하고 뺄 것이 없다. f_BodyToNed 는 이 행렬을 곱하기만 하고 f_NedToBody 는 전치한 행렬을 곱한다.

4단계 — NED ↔ ECEF
1
2
3
C_ecef_from_ned = Rz(경도) · Ry(−90° − 위도)
p_ecef_target = C_ecef_from_ned × p_ned + p_ecef_platform
p_ned         = (C_ecef_from_ned)ᵀ × (p_ecef_target − p_ecef_platform)
1
2
3
4
5
6
7
8
9
10
11
12
13
// NED -> ECEF 회전. Rz(경도) * Ry(-90 - 위도)
static ST_Mat3 f_RotNedToEcef(const ST_Lla stOrigin)
{
    return f_MatMul(f_RotZ(stOrigin.dLon), f_RotY(-PI / 2.0 - stOrigin.dLat));
}

// NED 벡터를 지구 방향으로 돌리고 함선 ECEF 위치를 더함
ST_Vec3 f_NedToEcef(const ST_Lla stOrigin, const ST_Vec3 stNed)
{
    ST_Vec3 stRot = f_MatVec(f_RotNedToEcef(stOrigin), stNed);

    return f_VecAdd(stRot, f_LlaToEcef(stOrigin));
}

회전 두 번의 곱을 풀어 쓰면 1.5.1절에서 본, 북·동·아래 단위벡터를 열로 세운 행렬이 된다. 함선의 ECEF 위치는 함수 안에서 f_LlaToEcef 로 구하므로 호출하는 쪽은 함선의 위도·경도·고도만 넘기면 된다. f_EcefToNed 는 반대로 함선 위치를 먼저 빼고 전치한 행렬로 돌린다.

5단계 — ECEF ↔ LLA

위도·경도·고도에서 ECEF 로 가는 쪽은 1.6절의 식을 그대로 옮기면 된다.

1
2
3
4
N = a / √(1 − e²·sin²(위도))
X = (N + 고도)·cos(위도)·cos(경도)
Y = (N + 고도)·cos(위도)·sin(경도)
Z = (N·(1 − e²) + 고도)·sin(위도)
1
2
3
4
5
6
7
8
9
10
11
12
ST_Vec3 f_LlaToEcef(const ST_Lla stLla)
{
    ST_Vec3 stEcef;
    FLOAT64 dN;     // 그 위도에서 본 지구 곡률반경 (동서 방향). 적도에서 가장 작고 극에서 가장 큼

    dN = WGS84_A / sqrt(1.0 - WGS84_E2 * sin(stLla.dLat) * sin(stLla.dLat));

    stEcef.dX = (dN + stLla.dAlt) * cos(stLla.dLat) * cos(stLla.dLon);
    stEcef.dY = (dN + stLla.dAlt) * cos(stLla.dLat) * sin(stLla.dLon);
    stEcef.dZ = (dN * (1.0 - WGS84_E2) + stLla.dAlt) * sin(stLla.dLat);   // z 만 (1 - e^2) 가 붙음. 지구가 극 쪽으로 눌린 만큼 줄이는 것
    return stEcef;
}

반대 방향은 한 번에 풀리지 않아서 1.6절의 반복 계산을 그대로 옮겼다.

1
2
3
ρ = √(X² + Y²)                          경도 = atan2(Y, X)
N = a / √(1 − e²·sin²(위도₀))          고도 = ρ / cos(위도₀) − N          (위도₀ 는 직전 회차의 위도)
위도 = atan2( Z·(N + 고도),  ρ·(N·(1 − e²) + 고도) )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
ST_Lla f_EcefToLla(const ST_Vec3 stEcef)
{
    ST_Lla  stLla;
    FLOAT64 dP = sqrt(stEcef.dX * stEcef.dX + stEcef.dY * stEcef.dY);   // 지구 자전축에서 떨어진 거리
    FLOAT64 dLat0;
    FLOAT64 dAlt0;
    FLOAT64 dLat;
    FLOAT64 dAlt;
    FLOAT64 dN;
    INT32   i;

    stLla.dLon = atan2(stEcef.dY, stEcef.dX);

    // 초기값은 지구를 그냥 공으로 봤을 때 위도, 고도는 0
    dLat0 = atan2(stEcef.dZ, dP);
    dAlt0 = 0.0;
    dLat  = dLat0;
    dAlt  = dAlt0;

    for (i = 0; i < LLA_ITER_MAX; i++)
    {
        dN   = WGS84_A / sqrt(1.0 - WGS84_E2 * sin(dLat0) * sin(dLat0));
        dAlt = dP / cos(dLat0) - dN;
        dLat = atan2(stEcef.dZ * (dN + dAlt), dP * (dN * (1.0 - WGS84_E2) + dAlt));

        if ((fabs(dLat - dLat0) < LLA_TOL_LAT) && (fabs(dAlt - dAlt0) < LLA_TOL_ALT))
        {
            break;
        }
        dLat0 = dLat;
        dAlt0 = dAlt;
    }

    stLla.dLat = dLat;
    stLla.dAlt = dAlt;
    return stLla;
}

위도 변화가 LLA_TOL_LAT(1e−12 rad), 고도 변화가 LLA_TOL_ALT(1e−6 m) 아래로 내려가면 멈추고, 그렇지 않더라도 LLA_ITER_MAX(10 회)에서 끝난다. 예제의 표적 점에서는 1.6절의 표처럼 여섯 번째에 멈춘다. dP / cos(dLat0) 는 극점에서 0/0 이 되는데 이 코드는 극점을 따로 처리하지 않는다.

3.1.4 실행 결과

입력

구분값
Slant Range20 km
Azimuth30° (보어사이트에서 왼쪽)
Elevation0°
플랫폼 위도 / 경도 / 고도36.408° / 127.307° / 0 m
플랫폼 자세 (roll, yaw, pitch)0° / 45° / 0°
안테나 설치선수 기준 시계방향 90°, 30 m 뒤 / 10 m 위

전체 실행 출력

main.c 는 측정값을 f_PolarToXyz, f_AntToBody, f_BodyToNed, f_NedToEcef, f_EcefToLla 에 차례로 넣어 위도·경도·고도를 얻고, 그 결과를 f_LlaToEcef, f_EcefToNed, f_NedToBody, f_BodyToAnt, f_XyzToPolar 에 차례로 넣어 극좌표로 되돌린다. 아래는 단계마다 나온 값을 그대로 찍은 것이다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
[입력]
  측정값       R = 20000.0 m, Az = 30.000 deg, El = 0.000 deg
  플랫폼 위치  lat = 36.408000 deg, lon = 127.307000 deg, alt = 0.0 m
  플랫폼 자세  (r, y, p) = (0.000, 45.000, 0.000) deg
  안테나 설치  az = 90.000 deg, tilt = 0.000 deg, offset = (-30.0, 0.0, -10.0) m

[정변환]
  1) 극좌표 -> 안테나 직교 (x, y, z)
            10000.0000           0.0000       17320.5081  [m]
  2) 안테나 -> 동체 (x, y, z)
             9970.0000       17320.5081         -10.0000  [m]
  3) 동체 -> NED (N, E, D)
            -5197.5941       19297.3033         -10.0000  [m]
  4) NED -> ECEF (X, Y, Z)
         -3132053.7872     4078527.1919     3760545.9526  [m]
     함선 ECEF (X, Y, Z)
         -3114830.1107     4087762.8525     3764723.0978  [m]
  5) ECEF -> LLA
      lat = 36.360967176 deg, lon = 127.522008441 deg, alt = 41.2824 m

[역변환]
  1) LLA -> ECEF (X, Y, Z)
         -3132053.7872     4078527.1919     3760545.9526  [m]
  2) ECEF -> NED (N, E, D)
            -5197.5941       19297.3033         -10.0000  [m]
  3) NED -> 동체 (x, y, z)
             9970.0000       17320.5081         -10.0000  [m]
  4) 동체 -> 안테나 직교 (x, y, z)
            10000.0000           0.0000       17320.5081  [m]
  5) 안테나 직교 -> 극좌표
      R = 20000.000000 m, Az = 30.000000000 deg, El = 0.000000000 deg
      입력과 차이  dR = 6.8e-10 m, dAz = 6.8e-13 deg, dEl = 3.8e-12 deg

[최종 출력]
  1) 표적 위/경/고도     lat = 36.360967 deg, lon = 127.522008 deg, alt = 41.3 m
  2) 표적 안테나 극좌표  R = 20000.0 m, Az = 30.000 deg, El = 0.000 deg

3.1.5 단계별 결과 해설

1단계 결과
1
p_antenna = ( +10 000.0000,  0.0000,  +17 320.5081 ) m        (x 왼쪽, y 위, z 보어사이트)

거리 20 km 가 왼쪽 성분 20000 × sin30° = 10000 과 보어사이트 성분 20000 × cos30° = 17320.5 로 나뉜다. 고각이 0° 라 위 성분은 0 이다. 손으로 확인할 수 있는 값이다.

2단계 결과
1
2
3
4
축 맞춤 A  : ( +17 320.5081,  −10 000.0000,    0.0000 )   ← 정면, 오른쪽, 아래 순서로 다시 읽은 값
Rz(90°)    : ( +10 000.0000,  +17 320.5081,    0.0000 )   ← 두 숫자가 자리를 바꾸고 부호 하나가 뒤집혔다
+ 레버암   : (  +9 970.0000,  +17 320.5081,  −10.0000 )
선수 기준 방위 60.074485°,  수평거리 19 985.0169 m

실행 출력에는 셋째 줄만 찍힌다. 둘째 줄은 f_AntToBody 안의 stRot 이고, 첫째 줄은 축 맞춤 행렬만 따로 곱해 본 값이다. 방위와 수평거리는 셋째 줄의 선수 성분과 우현 성분에서 atan2(우현, 선수) 와 √(선수² + 우현²) 로 계산했다.

안테나 정면은 선수에서 오른쪽으로 90° 이고 표적은 그 정면에서 선수 쪽으로 30° 되돌아온 방향에 있다 (그림 2-1). 그래서 90° − 30° = 60° 가 예상값인데 실제로는 60.074°가 나왔다. 차이 0.074°는 레버암 30 m가 만든 것이다. 20 km 거리에서 30 m는 각도로 atan(30/20000) ≈ 0.086° 이다. 레버암의 수평 성분은 선미 쪽을 향하고 표적은 선수에서 60° 방향이라, 30 m 가운데 표적 방향에 수직인 성분만 방위를 바꾼다. 그래서 0.086° 에 sin 60° 를 곱한 약 0.074° 가 된다.

3단계 결과
1
2
p_ned = ( N −5 197.5941,  E +19 297.3033,  D −10.0000 ) m
진북 기준 방위 105.074485°,  수평거리 19 985.0169 m

45°(yaw) + 90°(설치) − 30°(측정) = 105° 와 일치한다(레버암분 0.074° 제외). 방위각은 왼쪽이 + 라서 더하지 않고 뺀다. D가 −10 이므로 표적은 로컬 평면보다 10 m 위에 있다.

4·5단계 결과
1
2
3
4
함선 ECEF           = ( −3 114 830.1107,  +4 087 762.8525,  +3 764 723.0978 ) m
로컬 벡터를 돌린 값 = (    −17 223.6765,      −9 235.6606,      −4 177.1452 ) m
표적 ECEF (둘의 합) = ( −3 132 053.7872,  +4 078 527.1919,  +3 760 545.9526 ) m
표적 LLA            = 36.360967176°N,  127.522008441°E,  41.2824 m

둘째 줄은 실행 출력에 없는 값으로, f_NedToEcef 안의 stRot 이다.

고도가 왜 10 m가 아니라 41.3 m인가. 레이다는 수평(El = 0°)으로 쐈고 안테나는 10 m 높이에 있으니 “표적 고도 10 m”일 것 같다. 그런데 빔은 직진하는 반면 바다 표면은 멀어질수록 아래로 휘어 내려간다. 20 km 지점에서 그 차이가 31.3 m이고, 10 + 31.3 = 41.3 m가 된다. 로컬 평면 좌표를 그대로 고도로 쓰면 안 되는 이유가 바로 이것이다.

지구 곡률 때문에 생기는 고도 차이 그림 3-1. 빔은 직진하지만 바다 표면은 멀어질수록 내려간다. 그 차이가 거리의 제곱에 비례해 커진다.

(실제로는 대기 굴절 때문에 빔도 아래로 살짝 휘어서 이보다 조금 작아진다. 정밀 고도 산출에서는 4/3 지구 모델 등으로 굴절을 함께 고려한다.)

3.1.6 왕복 확인

좌표변환은 “적당히 그럴듯한 값”이 나오기 때문에 눈으로 봐서는 검증되지 않는다. 그래서 갔다가 되돌아왔을 때 처음 값이 나오는지를 본다. main.c 는 정변환으로 얻은 표적의 위도·경도·고도를 역변환 함수 다섯 개에 차례로 넣어 안테나 극좌표까지 되돌린 다음, 그 거리·방위각·고각을 입력값과 비교한다. 실행 출력의 [역변환] 블록에는 [정변환] 블록의 단계별 값이 반대 순서로 그대로 다시 나온다.

자세 (roll, pitch, yaw)거리 차이방위각 차이고각 차이
(0°, 0°, 45°) 예제 자세6.8e−10 m6.8e−13°3.8e−12°
(20°, 10°, 45°)5.9e−10 m2.2e−12°1.1e−12°

첫 행이 실행 출력의 값이다. 이 차이는 계산 도중의 반올림에서 나온다. ECEF 좌표는 6.4e6 m 규모이고 배정밀도 실수(double)의 상대정밀도는 2.2e−16 이라서, 둘을 곱한 약 1.4e−9 m 가 ECEF 좌표 하나에 실리는 반올림 오차의 규모다. 거리 차이 6.8e−10 m 는 그 안에 있고, 고각 차이 3.8e−12° 도 20 km 거리에서 길이로 바꾸면 약 1.3e−9 m 로 같은 수준이다.

ECEF에서 미터 이하를 다루면 배정밀도의 유효 자릿수 절반을 이미 써 버린 상태가 된다. 그래서 가까운 두 점의 상대 위치는 반드시 로컬 좌표계에서 계산해야 한다.

둘째 행은 자세만 바꿔 같은 왕복을 다시 돌린 결과다. 왕복 차이는 예제 자세와 같은 수준이다.

왜 자세를 전부 다른 값으로 두는가 예제 자세는 roll 과 pitch 가 0 이라 Rx 와 Ry 가 단위행렬이고, 세 회전 가운데 yaw 하나만 실제로 돈다. 되돌아가는 행렬을 따로 풀어 적으면서 곱하는 순서를 뒤집지 않은 실수는 이런 자세에서 드러나지 않는다. 예를 들어 (Rz·Ry·Rx)ᵀ 는 Rxᵀ·Ryᵀ·Rzᵀ 인데, 이것을 Rzᵀ·Ryᵀ·Rxᵀ 로 잘못 풀어도 roll = pitch = 0 에서는 값이 같다. roll·pitch·yaw를 전부 다른 값으로 두어야 왕복 차이로 나타난다.

왕복이 확인해 주는 것은 정변환과 역변환이 서로 역이라는 것까지다. 이 코드는 되돌아가는 회전행렬을 가는 쪽 행렬의 전치로 만들기 때문에, 가는 쪽 행렬을 곱하는 순서가 틀렸거나 회전 방향이 반대여도 왕복은 그대로 통과한다. 그래서 3.1.5절의 60° 와 105° 처럼 손으로 계산할 수 있는 값과도 맞춰 봐야 한다.

3.1.7 최종 출력값

조건 ⑥ 이 요구한 두 가지 출력

[1] 표적 위도 / 경도 / 고도
항목값
위도36.360967° N
경도127.522008° E
고도41.3 m (WGS-84 타원체고)

소수 아홉 자리까지의 값은 3.1.4절의 실행 출력에 있다.

[2] 표적 안테나 기준 시선거리 / 방위각 / 고각
항목값입력값과의 차이
시선거리20 000.0 m6.8e−10 m
방위각30.000°6.8e−13°
고각0.000°3.8e−12°

추적 필터를 넣지 않았으므로 역변환 결과는 입력값과 같아야 하고, 실제로 컴퓨터 반올림 수준까지 일치한다.


3.2 추가 확인 실험

설계에서 강조한 세 가지 판단이 실제로 근거가 있는지 코드로 확인했다. 세 실험 모두 3.1절의 변환 함수에 입력값만 바꿔 넣은 결과다.

실험 A — 레버암을 빼먹으면

 위도경도고도
레버암 반영36.360967176°127.522008441°41.28 m
레버암 무시36.361157843°127.522245658°31.33 m
차이약 21 m약 21 m약 10 m

레버암을 반영한 점에서 무시한 점을 로컬 좌표로 빼면 (N, E, D) = (−21.213, −21.213, −10.000) m 다. 수평으로 정확히 30 m, 높이로 10 m 이니 레버암 그 자체다. 이 오차는 거리가 멀어져도 줄지 않는다.

실험 B — 배가 기울면 각도 덧셈이 깨진다

배가 기울면 빔이 실제로 향하는 방향 그림 3-2. 안테나는 배에 고정되어 있어서 우현이 내려가도록 배가 기울면 우현을 보는 빔도 바다 쪽으로 내려간다. 레이다는 여전히 “고각 0도”라고 보고하지만 그건 안테나 기준일 뿐이다.

같은 측정값(Az 30°, El 0°)에 자세만 바꿔 넣은 결과다. 단순 덧셈은 yaw + 90° − 30° 이고, 실제 진북 방위 atan2(E, N) 과 실제 고각 atan2(−D, √(N² + E²)) 은 변환해서 얻은 NED 좌표에서 구했다.

배 자세 (roll, pitch, yaw)단순 덧셈실제 진북 방위실제 고각표적 고도
(0°, 0°, 45°)105°105.07°+0.03°41.3 m
(0°, 0°, 60°)120°120.07°+0.03°41.3 m
(10°, 0°, 45°)105°104.70°−8.63°−2 967 m
(20°, 0°, 45°)105°103.52°−17.21°−5 886 m
(0°, 10°, 45°)105°105.46°+5.00°+1 772 m
(20°, 10°, 45°)105°101.33°−11.82°−4 063 m

roll이 20°만 되어도 고각이 0°에서 −17.21° 로 바뀐다. 레이다는 여전히 “고각 0도”라고 보고하지만 그건 안테나 기준일 뿐이고, 세상 기준으로는 표적 방향이 바다 쪽으로 수평보다 17도 아래라는 뜻이다. 20 km 거리에서 이는 고도로 약 5.9 km 차이다(표의 표적 고도가 음수인 이유). 기운 각 20° 가 전부 실리지 않고 17° 남짓인 것은 표적이 안테나 정면에서 30° 벗어나 있어서다.

pitch +10° 에서는 반대로 고각이 +5.00° 로 올라간다. 선수가 들리는 자세인데 표적이 선수에서 오른쪽으로 60° 방향, 곧 배의 앞쪽 절반에 있어서 선수와 함께 들리기 때문이다. yaw 만 바꾼 둘째 행은 단순 덧셈과의 차이가 첫 행과 같은 0.07° 그대로이고, roll 과 pitch 가 함께 들어간 마지막 행은 방위도 단순 덧셈에서 약 3.7° 벗어난다. 첫 두 행의 +0.03° 는 안테나가 무게중심보다 10 m 높아서 생긴 값이다.

각도를 더하는 방식으로는 이 효과를 만들어낼 수 없다. 세 축의 회전이 서로 섞이기 때문이다. 행렬 곱셈이 이 섞임을 자동으로 처리해 준다. 이것이 2.1.4의 3단계가 설계의 핵심인 이유다.

실험 C — 지구 곡률

거리로컬 평면 기준 고도측지 고도차이
5 km10.00 m11.95 m1.95 m
10 km10.00 m17.81 m7.81 m
20 km10.00 m41.28 m31.28 m
50 km10.00 m205.69 m195.69 m
100 km10.00 m792.96 m782.96 m

차이가 거리의 제곱에 비례해서 커진다. 1.4절의 표에 어림값으로 적어 둔 것이 이 값이다. 로컬 평면의 D 성분을 고도로 그대로 쓸 수 없고, 고도는 ECEF 를 거쳐 LLA 로 계산한 값을 써야 한다.


부록 A. 자주 하는 실수

#실수증상예방책
1z 방향 혼동 (위 vs 아래)고도 부호가 뒤집힌다. 표적이 땅속에 있다동체·NED는 z를 아래로 통일. “10 m 위 = z −10”. 안테나 면 기준 축은 축 맞춤 행렬 한 곳에서만 바꾼다
2전치 방향 뒤바뀜자세가 0일 때는 정상, 기울면 틀림함수 이름에 방향 명시 (f_BodyToNed). roll·pitch·yaw를 전부 다른 값으로 시험
3라디안/도 혼용결과가 57배 또는 1/57배로 엉뚱계산·저장은 라디안으로 통일. 도는 입출력 경계에서 한 번만 변환 (DEG2RAD, RAD2DEG)
4Z = (N·(1−e²)+h)·sin(lat) 에서 (1−e²) 누락위도가 조금씩 틀림. 적도에선 안 보이고 위도 45° 부근에서 가장 큼 (약 0.19°)LLA→ECEF→LLA 왕복 시험
5atan(a/b) 사용분모 b가 음수일 때 방위가 180° 틀림항상 atan2(a, b)
6asin 인자 클램프 누락특정 입력에서 NaN이 나오고 조용히 번짐atan2 형태로 바꾸거나 클램프
7레버암 누락결과가 항상 일정량 어긋남 (거리 무관)레버암을 넣은 결과와 뺀 결과를 비교해 차이가 레버암만큼인지 확인 (3.2절 실험 A)
8로컬 평면 고도를 그대로 사용먼 표적일수록 고도가 낮게 나옴반드시 LLA로 환산
9타원체고와 해발고도 혼동고도가 20 m 남짓 계통 오차내부는 타원체고로 통일, 표시만 변환
10오일러각 순서 불일치큰 자세각에서만 틀림3-2-1로 고정하고 문서에 명시
11방위각 부호 규약 혼동 (왼쪽이 + vs 오른쪽이 +)표적이 보어사이트를 사이에 두고 반대편에 찍힌다. 거리는 그대로라 값이 그럴듯해 보인다어느 쪽이 +인지를 식과 함께 문서에 적는다. 이 글에서는 왼쪽이 + (x = R·cos(El)·sin(Az), x는 왼쪽)

부록 B. 용어집

용어뜻
보어사이트 (Boresight)안테나가 정면으로 바라보는 방향. 안테나 면의 수직 방향
시선거리 (Slant Range)안테나에서 표적까지의 직선 거리. 지표를 따라 잰 거리가 아니다
방위각 (Azimuth)기준 방향에서 좌우로 잰 각. 기준이 보어사이트인지 진북인지, 어느 쪽이 +인지 항상 명시할 것. 이 글에서 안테나 방위각은 보어사이트에서 왼쪽이 +, 진북 기준 방위는 시계방향이 +다
고각 (Elevation)수평면에서 위아래로 잰 각
레버암 (Lever arm)무게중심에서 센서까지의 위치 차이 벡터
백틸트 (Back tilt)안테나 면을 뒤로 젖혀 보어사이트를 위로 들어올린 각
축 맞춤 행렬 (A)안테나 면 기준 축(x 왼쪽, y 위, z 보어사이트)으로 적은 좌표를 정면·오른쪽·아래 순서로 바꿔 적는 상수 행렬. A × (x, y, z) = (z, −x, −y). 각도가 들어가지 않고, 설치 방위·백틸트보다 먼저 곱한다
DCM (Direction Cosine Matrix)방향코사인 행렬 = 회전행렬. 좌표계 사이의 회전 관계를 담은 3×3 행렬
전치 (Transpose)행렬의 행과 열을 바꾸는 것. 회전행렬에서는 곧 역회전
직교행렬Cᵀ·C = I 를 만족하는 행렬. 회전행렬은 항상 직교행렬이다
FRDForward-Right-Down. 앞-오른쪽-아래 축 배치
NED / ENUNorth-East-Down / East-North-Up. 국지수평 좌표계의 두 관습
ECEFEarth-Centered Earth-Fixed. 지구중심 지구고정 직교좌표계
ECIEarth-Centered Inertial. 지구중심 관성좌표계 (자전하지 않음)
LLALatitude-Longitude-Altitude. 위경도 고도 좌표계
WGS-84지구 타원체의 국제 표준. GPS·해도가 쓴다
타원체고 (HAE)타원체 표면에서 잰 높이. LLA의 고도가 이것
지오이드 (Geoid)중력이 만드는 평균 해수면. 해발고도의 기준
INSInertial Navigation System. 관성항법장치
정렬 (Alignment)센서 축과 동체 축의 어긋남을 추정·보정하는 절차
GMSTGreenwich Mean Sidereal Time. ECEF↔ECI 회전각을 주는 시각
방향코사인 (Direction cosine, u·v·w)단위 방향벡터를 각 축에 정사영한 성분. UV 좌표가 이것
가시영역 (Visible region)u² + v² ≤ 1 인 원판. 이 밖으로는 빔을 만들 수 없다

부록 C. 오차는 어디서 오는가

좌표변환을 “정확히” 한다는 것이 무슨 뜻인지 감을 잡으려면, 각 단계의 오차가 최종 표적 위치에 얼마나 영향을 주는지 알아야 한다. 각도 오차는 거리에 비례해서 커지고, 위치 오차는 거리와 무관하게 그대로 남는다 — 이 두 문장이 핵심이다.

오차원어느 단계최종 위치에 미치는 영향20 km 표적 기준
안테나 각도 측정 오차 1 mrad1단계거리 × 1e−320 m
안테나 설치 정렬 오차 1 mrad2단계거리 × 1e−320 m
레버암 누락 30 m2단계거리와 무관하게 30 m30 m
INS 자세 오차 1 mrad3단계거리 × 1e−320 m
자세 데이터 지연 10 ms (roll 10°/s)3단계거리 × 1.75e−335 m
INS 위치 오차 10 m4단계거리와 무관하게 10 m10 m
지구 곡률 무시5단계거리² / (2R) — 고도 성분31 m (고도)
지오이드고 무시5단계고도에 계통 오차20~25 m (고도)

여기서 읽어야 할 것

  1. 1 mrad = 0.0573° 이고, 이것이 거리 20 km에서 20 m가 된다. “1도도 안 되는 작은 각”이 실제로는 큰 위치 오차라는 감각을 갖는 것이 중요하다.
  2. 자세 지연이 자세 정확도보다 클 수 있다. 표에서 지연 10 ms의 기여(35 m)가 INS 자세 오차 1 mrad(20 m)보다 크다. 좋은 INS를 쓰고도 시각 처리를 잘못하면 그 성능이 그대로 사라진다.
  3. 레버암과 INS 위치 오차는 거리가 멀어져도 줄지 않는다. 원거리 표적에서는 각도 성분에 묻히지만, 근거리 표적이나 정밀 사격통제에서는 지배적인 항이 된다.
  4. 고도 오차는 별도로 관리해야 한다. 곡률과 지오이드는 수평 위치에는 거의 영향을 주지 않지만 고도에는 수십 미터씩 실린다.

이 표는 규모 감각을 잡기 위한 것이고, 실제 오차 예산은 각 항을 제곱합해서(RSS) 합산하고 상관관계를 고려해야 한다. 다만 어느 항이 지배적인지를 먼저 아는 것이 설계에서 훨씬 중요하다.


부록 D. 코드

3절에서 쓴 코드는 아래 세 파일이다. 3.2절의 실험은 같은 함수에 입력값만 바꿔 넣어 돌렸다. 변환 코드(coord_frames.c)가 쓰는 외부 함수는 표준 C 수학 함수(<math.h>)뿐이고, Visual Studio 2022(MSVC)에서 /utf-8 /W4 로 경고 없이 빌드된다.

파일내용
coord_frames.h구조체 · WGS-84 상수 · 함수 선언
coord_frames.c3×3 행렬 연산, 단계별 정변환·역변환 함수 구현
main.c예제 입력값으로 정변환과 역변환을 돌려 단계별 값을 출력

변환 함수는 coord_frames.h 와 coord_frames.c 두 파일에 다 들어 있어서, 이 둘만 가져가면 다른 프로젝트에도 그대로 붙일 수 있다. 좌표계 이름을 함수 이름과 인자 이름에 넣어서(f_BodyToNed, stBody 같은 식) 이름만 봐도 어느 좌표계에서 어느 좌표계로 가는지 보이도록 했다.


요약

  1. 좌표계는 “누구 입장에서 말하는가” 이고, 좌표변환은 그 통역이다.
  2. 변환의 재료는 회전과 평행이동 둘뿐이다.
  3. 이번 설계의 체인은 안테나 → 동체 → 로컬(NED) → ECEF → LLA 이고, 각 단계는 앞 단계 없이는 성립하지 않는다. 안테나에서 동체로 넘어갈 때는 안테나 면 기준 축을 정면·오른쪽·아래로 바꿔 적는 축 맞춤 행렬을 가장 먼저 곱한다.
  4. 회전은 설치 도면(설치 방위·백틸트)과 INS의 자세·위치가, 평행이동은 설치 도면(레버암)과 배의 위치가 공급한다.
  5. 조건 ③(INS 축 일치) 덕분에 INS 좌표계는 별도 단계로 두지 않는다.
  6. 돌아올 때는 전치하고, 빼고, 순서를 뒤집는다. 왕복 결과가 원래 값이면 정변환과 역변환이 서로 역이라는 것까지 확인된다. 가는 쪽이 맞는지는 손으로 계산한 값(60°, 105°)과 따로 맞춰 본다.
  7. 고정형 안테나는 배의 흔들림을 전부 계산으로 보상한다. 그래서 회전행렬이 필수다 — roll 20°에서 고각이 0°→ −17°로 바뀐다.
  8. 로컬 평면 좌표를 고도로 그대로 쓰면 안 된다. 20 km에서 31.3 m가 어긋난다.

최종 결과

1
2
3
4
5
6
7
8
9
[1] 표적 위/경/고도
      위도  36.360967176 °N
      경도  127.522008441 °E
      고도  41.2824 m  (WGS-84 타원체고)

[2] 표적 안테나 기준 극좌표
      시선거리  20 000.0000 m
      방위각    30.000000 °  (보어사이트에서 왼쪽)
      고각      0.000000 °
이 기사는 저작권자의 CC BY 4.0 라이센스를 따릅니다.