# 기본 측위 알고리즘(삼변측량 원리)

#### GNSS에서의 거리 기반 측위 개념

GNSS 위성으로부터의 신호를 수신해 사용자의 위치를 결정하는 핵심 과정은 기본적으로 ‘삼변측량(Trilateration)’에 기반한다. 삼각측량(Triangulation)이 각도 측정을 주로 이용하는 반면, 삼변측량은 서로 다른 세 지점(또는 그 이상)과의 **거리**를 측정하여 2차원 혹은 3차원 공간에서 수신기의 위치를 추정한다. GNSS에서는 여러 위성으로부터 전파가 도달하는 시간을 측정하고, 전파의 전파 속도(진공 중 빛의 속도)를 고려해 거리(또는 의사거리, pseudorange)를 추정한다. 이 거리를 정확히 측정하면, 하나의 위성으로부터는 구(球) 형태의 위치 후보가, 두 위성으로부터는 두 구의 교집합인 원(圓) 형태가, 세 위성으로부터는 그 교집합인 두 개의 점이 도출된다. 이를 적절히 해석하고, 네 번째 위성을 추가적으로 사용하여 시계(Clock) 바이어스 등을 보정하면 3차원 공간에서 사용자 위치를 결정할 수 있다.

아래 예시는 간단히 3차원 공간에서 네 개 이상의 위성으로부터의 거리 측정 정보를 이용하는 구도를 나타낸 것이다:

{% @mermaid/diagram content="graph LR
S1(위성 1) --> U(사용자)
S2(위성 2) --> U
S3(위성 3) --> U
S4(위성 4) --> U" %}

#### 기본 방정식 정립

사용자 위치를 $\mathbf{x}$라 하고, $i$번째 위성의 위치를 $\mathbf{p}\_i$라 하자. 이상적으로는 $i$번째 위성으로부터 사용자의 거리 $\rho\_i$가 다음과 같이 정의된다.

$$
\rho\_i = | \mathbf{x} - \mathbf{p}\_i |
$$

여기서 $\mathbf{x} = \begin{bmatrix} x \ y \ z \end{bmatrix}$ 이고, $\mathbf{p}\_i = \begin{bmatrix} x\_i \ y\_i \ z\_i \end{bmatrix}$ 이다. 또한 $|\cdot|$는 유클리드 노름(Euclidean Norm)을 의미한다. 그러나 실제 GNSS에서는 사용자의 수신기 시계 오차 때문에 측정된 거리 $\rho\_i$가 참값과 완전히 일치하지 않고, 다음과 같은 의사거리(pseudorange) $P\_i$를 측정한다고 본다.

$$
P\_i = | \mathbf{x} - \mathbf{p}\_i | + c,\delta t
$$

여기서

* $c$는 빛의 속도(전파 속도),
* $\delta t$는 수신기의 시계 오차를 의미한다. 시계 오차를 포함하여 삼변측량을 수행하기 위해서는 최소 4개의 위성 측정값이 필요하다. 각각의 의사거리 방정식은 아래와 같이 표현된다.

$$
P\_i - | \mathbf{x} - \mathbf{p}\_i | = c,\delta t, \quad i = 1,2,3,4
$$

위 네 개 방정식을 해 $\mathbf{x}$와 $\delta t$에 대해 풀면, 이론적으로는 $\mathbf{x}$를 고유하게 구할 수 있다.

#### 비선형 방정식 구조

삼변측량 기본 방정식은 아래와 같은 비선형 형태를 갖는다.

$$
f\_i(\mathbf{x}, \delta t) = P\_i - | \mathbf{x} - \mathbf{p}\_i | - c, \delta t = 0, \quad i = 1,2,3,4
$$

이를 직접 해석적으로 풀면 단순하지 않기 때문에, 보통은 뉴턴-랩슨(Newton-Raphson)이나 가우스-뉴턴(Gauss-Newton) 알고리즘과 같은 수치해석 기법을 이용한다. 이때 초기값 설정이 중요한 역할을 하며, 관측된 위성들의 위치 정보와 대략적인 사용자 위치 추정값을 바탕으로 점진적으로 해를 개선한다. 해를 구하기 위한 반복 과정에서 야코비(Jacobian) 행렬은 다음과 유사하게 구성될 수 있다. $\mathbf{x} = \[x ; y ; z]^\mathsf{T}$라 할 때,

$$
\mathbf{J}\_{(i,:)} = \frac{\partial f\_i(\mathbf{x}, \delta t)}{\partial \[x ; y ; z ; \delta t]}
$$

와 같이 정의된다. 구체적으로 $| \mathbf{x} - \mathbf{p}\_i |$의 편도함수들은 다음과 같은 꼴을 갖는다.

$$
\frac{\partial}{\partial x} | \mathbf{x} - \mathbf{p}\_i | =  \frac{x - x\_i}{\sqrt{(x - x\_i)^2 + (y - y\_i)^2 + (z - z\_i)^2}}
$$

( $y$, $z$ 방향 편도함수도 유사한 구조 )

#### 보정 항과 지표면 좌표 변환

추가로, GNSS 시스템에서는 전리층 및 대류권 지연, 다중경로(multipath), 안테나 오프셋, 위성 궤도력(력학 모델) 등에 의해 측정된 의사거리 $P\_i$가 편차를 갖는다. 따라서 실제 측정방정식은

$$
P\_i = | \mathbf{x} - \mathbf{p}*i | + c,\delta t + d*{\text{ionosphere},i} + d\_{\text{troposphere},i} + \varepsilon\_i
$$

와 같이 더 복잡해진다. 여기서 $d\_{\text{ionosphere},i}$와 $d\_{\text{troposphere},i}$는 각각 전리층 지연, 대류권 지연을 의미하고, $\varepsilon\_i$는 기타 잡음(noise) 및 오차 항을 포괄한다. 수신기 내부 알고리즘이나 보정 서비스(SBAS, GBAS, PPP 등)를 통해 이러한 오차 항들을 추정ㆍ보정하여, 이상 방정식에 최대한 가깝게 맞춰 나가는 과정을 거친다.

#### 선형화 기법

반복 알고리즘 구현 시에는, 목표 함수를 선형 근사로 단순화하여(테일러 전개 등 활용) 해를 갱신한다. 예를 들어, 현재 추정치 $(\mathbf{x}^{(k)}, \delta t^{(k)})$에서의 잔차(residual) 벡터 $\mathbf{r}^{(k)}$를 정의하면,

$$
\mathbf{r}^{(k)} = \begin{bmatrix} r\_1^{(k)} \ r\_2^{(k)} \ r\_3^{(k)} \ r\_4^{(k)} \end{bmatrix} =  \begin{bmatrix} P\_1 - | \mathbf{x}^{(k)} - \mathbf{p}\_1 | - c,\delta t^{(k)} \ P\_2 - | \mathbf{x}^{(k)} - \mathbf{p}\_2 | - c,\delta t^{(k)} \ P\_3 - | \mathbf{x}^{(k)} - \mathbf{p}\_3 | - c,\delta t^{(k)} \ P\_4 - | \mathbf{x}^{(k)} - \mathbf{p}\_4 | - c,\delta t^{(k)} \end{bmatrix}
$$

그리고 야코비 행렬 $\mathbf{J}^{(k)}$를 구한 뒤, 선형 근사를 통해

$$
\Delta \mathbf{u}^{(k)} = \begin{bmatrix} \Delta \mathbf{x}^{(k)} \ \Delta \delta t^{(k)} \end{bmatrix} = (\mathbf{J}^{(k)})^\dagger , \mathbf{r}^{(k)}
$$

와 같이 추정치 보정 벡터를 계산할 수 있다. 여기서 $(\cdot)^\dagger$는 의사역행렬(pseudoinverse)을 의미한다. 그다음

$$
\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \Delta \mathbf{x}^{(k)}, \quad \delta t^{(k+1)} = \delta t^{(k)} + \Delta \delta t^{(k)}
$$

와 같이 추정치를 갱신한다. 이 과정을 잔차 벡터가 충분히 작아질 때까지 반복한다.

#### 반복 알고리즘 구현 예시 (가우스-뉴턴 기법)

가우스-뉴턴(Gauss-Newton) 방식은 비선형 최소제곱(Nonlinear Least Squares) 문제를 다룰 때 자주 사용된다. GNSS의 삼변측량 문제에서도 측정 방정식이 비선형이므로 가우스-뉴턴 알고리즘을 활용해 사용자 위치와 시계 오차를 추정할 수 있다. 다음과 같은 최소제곱 문제를 고려한다.

$$
\min\_{\mathbf{x},,\delta t} \sum\_{i=1}^{N} \Bigl\[P\_i - |\mathbf{x} - \mathbf{p}\_i| - c,\delta t \Bigr]^2
$$

여기서 $N \ge 4$로, 사용 가능한 위성(측정값) 개수를 의미한다. 위 과정에서 측정 방정식의 오차가 최소화되도록 $(\mathbf{x}, \delta t)$를 찾으면 된다.

1. **초기값 설정**
   * $\mathbf{x}^{(0)}$를 대략적으로 설정한다. 일반적으로는 이전 에포크(epoch)의 해 또는 지상국에서 제공되는 초기 정보 등을 사용할 수 있다.
   * $\delta t^{(0)}$ 역시 대략적인 값을 부여한다(예: 0 혹은 다른 추정값).
2. **반복 과정**
   * (a) **잔차(residual) 계산**

     $$
     \mathbf{r}^{(k)} = \begin{bmatrix} r\_1^{(k)} \ r\_2^{(k)} \ \vdots \ r\_N^{(k)} \end{bmatrix} = \begin{bmatrix} P\_1 - |\mathbf{x}^{(k)} - \mathbf{p}\_1| - c,\delta t^{(k)} \ P\_2 - |\mathbf{x}^{(k)} - \mathbf{p}\_2| - c,\delta t^{(k)} \ \vdots \ P\_N - |\mathbf{x}^{(k)} - \mathbf{p}\_N| - c,\delta t^{(k)} \end{bmatrix}
     $$

     (b) **야코비(Jacobian) 행렬 계산**

     $$
     $\mathbf{J}^{(k)} \equiv \frac{\partial \mathbf{r}^{(k)}}{\partial \[x ; y ; z ; \delta t]}
     $$

     구체적으로 $i$번째 행 ($i \in {1, \dots, N}$)에 대해서,

     $$
     \frac{\partial}{\partial x} \Bigl\[P\_i - |\mathbf{x}^{(k)} - \mathbf{p}\_i| - c,\delta t^{(k)}\Bigr] = -\frac{(x^{(k)} - x\_i)}{|\mathbf{x}^{(k)} - \mathbf{p}\_i|}
     $$

     그리고 $y,, z$ 방향도 동일한 형태로 편미분할 수 있다. 마지막으로 $\delta t$에 대한 편미분은 $-c$가 된다.
   * (c) **업데이트 벡터 계산** 가우스-뉴턴에서는 선형 근사를 이용해 다음의 보정값(증분값)을 구한다.

     $$
     \Delta \mathbf{u}^{(k)} = \begin{bmatrix} \Delta \mathbf{x}^{(k)} \ \Delta \delta t^{(k)} \end{bmatrix} = \bigl((\mathbf{J}^{(k)})^\mathsf{T}\mathbf{J}^{(k)}\bigr)^{-1} , (\mathbf{J}^{(k)})^\mathsf{T} , \mathbf{r}^{(k)}
     $$

     여기서 $(\mathbf{J}^{(k)})^\mathsf{T}\mathbf{J}^{(k)}$이 역행렬이 가능하다는 전제(또는 유사역행렬 활용)가 필요하다.
   * (d) **추정치 갱신**

     $$
     \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \Delta \mathbf{x}^{(k)}, \quad \delta t^{(k+1)} = \delta t^{(k)} + \Delta \delta t^{(k)}
     $$

     (e) **수렴 판정** $|\Delta \mathbf{u}^{(k)}|$ 혹은 $|\mathbf{r}^{(k)}|$이 충분히 작으면 반복을 멈춘다.

#### 가중 최소제곱 (Weighted Least Squares)

여러 위성에서 측정되는 $\mathrm{C/N\_0}$ (carrier-to-noise ratio)나 예상 오차 분산 등을 고려해, 각각의 측정값에 차등적인 가중(weight)을 부여할 수 있다. 이 경우 최소화해야 할 목적함수는 다음과 유사하게 표현된다.

$$
\min\_{\mathbf{x},,\delta t} \sum\_{i=1}^{N} w\_i \Bigl\[P\_i - |\mathbf{x} - \mathbf{p}\_i| - c,\delta t \Bigr]^2
$$

여기서 $w\_i$는 해당 위성 측정값의 신뢰도를 반영하는 가중 계수다. 이 경우, 야코비 기반 반복식에서 $\mathbf{J}^{(k)}$와 $\mathbf{r}^{(k)}$가 가중 행렬 $\mathbf{W}$에 의해 조정되며,

$$
\Delta \mathbf{u}^{(k)}  = \bigl((\mathbf{J}^{(k)})^\mathsf{T}\mathbf{W},\mathbf{J}^{(k)}\bigr)^{-1} , (\mathbf{J}^{(k)})^\mathsf{T}\mathbf{W},\mathbf{r}^{(k)}.
$$

이와 같은 구조로 오차가 큰 위성 측정값은 자동으로 상대적으로 낮은 가중이 부여되고, 신뢰도가 높은 위성으로부터의 측정값이 해에 더 큰 영향을 미치게 된다.

#### 기하학적 해석 (구(球)의 교점 문제)

3차원 공간에서 위성 위치가 $\mathbf{p}\_1, \mathbf{p}\_2, \mathbf{p}\_3, \mathbf{p}\_4$라고 가정하자. 각 위성으로부터의 거리 측정값이 $\rho\_1, \rho\_2, \rho\_3, \rho\_4$라면, 이는 사실상 다음과 같은 4개의 구로 나타낼 수 있다.

1. 구 1: 중심 $\mathbf{p}\_1$, 반지름 $\rho\_1$
2. 구 2: 중심 $\mathbf{p}\_2$, 반지름 $\rho\_2$
3. 구 3: 중심 $\mathbf{p}\_3$, 반지름 $\rho\_3$
4. 구 4: 중심 $\mathbf{p}\_4$, 반지름 $\rho\_4$

이 4개의 구가 이상적으로 교차하는 지점은 일반적으로 0개, 1개, 2개 등 다양한 가능성이 존재한다. GNSS 시스템에서는 수신기의 시계 오차가 함께 존재하므로, 실제로는 하나의 후보점만 물리적으로 의미를 가질 수 있도록 추가 제약이 들어간다(시계 바이어스, 대략적 사용자 높이 등). 수치해석 알고리즘(예: 가우스-뉴턴, 뉴턴-랩슨)을 수행할 때는 위에서 언급한 ‘구들의 교점’을 찾는 과정을 반복적으로 근사하는 형태라고 볼 수 있다.

#### 추가 고려사항

* **Rank 조건**: $\mathbf{J}^{(k)}$의 칼럼들이 선형종속인 경우(또는 그에 준하는 상황) 역행렬 연산이 불가하므로, 알고리즘이 불안정해질 수 있다. 위성의 배치(기하학적 배치)가 좋지 않으면(Dilution of Precision, DOP이 높으면) 상태방정식 해석이 악화된다.
* **측정치 이상치(Outlier)**: 실측 데이터에서 일부 위성이 이상치(오류가 매우 큰 값)를 보낼 경우, 단순 최소제곱 해가 크게 왜곡될 위험이 있으므로 RANSAC 등 강인(Robust) 추정 기법을 적용하기도 한다.

#### 동적 사용자 시나리오에서의 필터링 기법 (확장 칼만 필터)

GNSS 수신기가 움직이는 물체(예: 차량, 항공기, 드론 등)와 같이 동적인 환경에서 사용될 때는, 매 에포크(epoch)마다 독립적으로 삼변측량 방정식을 풀어 위치를 구하는 것뿐만 아니라 ‘시간적으로 연속된 정보’를 적극 활용하는 방법이 더 효율적일 수 있다. 이때 자주 사용되는 방법 중 하나가 확장 칼만 필터(EKF, Extended Kalman Filter)다.

**상태벡터 설정**

EKF 적용 시, 추정해야 할 상태(state)를 예시적으로 다음과 같이 설정할 수 있다.

$$
\mathbf{x}\_k =  \begin{bmatrix} x\_k \ y\_k \ z\_k \ \dot{x}\_k \ \dot{y}\_k \ \dot{z}\_k \ \delta t\_k \ \dot{\delta t}\_k \end{bmatrix}
$$

* $x\_k,, y\_k,, z\_k$: $k$번째 시점에서의 위치
* $\dot{x}\_k,, \dot{y}\_k,, \dot{z}\_k$: 속도
* $\delta t\_k$: 시계 오차
* $\dot{\delta t}\_k$: 시계 오차 변화율(드리프트)

추가적으로 가속도나 자세(Attitude) 파라미터 등을 포함하여 상태를 확장할 수도 있다.

**프로세스 모델**

시스템(프로세스) 방정식은 일반적으로 다음과 같은 형태의 상태천이(transition)를 가정한다. 예를 들어, 짧은 시간간격 $\Delta t$ 후의 상태는

$$
\mathbf{x}\_{k+1} = \mathbf{F}\_k \mathbf{x}\_k + \mathbf{w}\_k
$$

로 나타낼 수 있다. 여기서

$$
\mathbf{F}\_k \approx \begin{bmatrix} 1 & 0 & 0 & \Delta t & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & \Delta t & 0 & 0 & 0 \ 0 & 0 & 1 & 0 & 0 & \Delta t & 0 & 0 \ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 1 & \Delta t \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}
$$

와 같이 위치-속도, 시계 오차-드리프트 사이의 선형 관계를 나타내고, $\mathbf{w}\_k$는 시스템 노이즈(process noise)에 해당한다.

**측정 방정식**

GNSS로부터 얻어지는 측정값(의사거리)은 앞서 언급된

$$
P\_i^k = |\mathbf{p}\_i^k - \mathbf{x}\_k^{pos}| + c ,\delta t\_k + \epsilon\_i^k
$$

와 유사한 형태를 갖는다. 여기서 $\mathbf{x}\_k^{pos} = \[,x\_k ; y\_k ; z\_k,]^\mathsf{T}$는 상태벡터 중 위치 성분만 떼어낸 것이고, $\mathbf{p}\_i^k$는 $i$번째 위성의 위치(알려진 값), $\epsilon\_i^k$는 오차 항을 의미한다. EKF 측면에서 측정 벡터 $\mathbf{z}\_k$를 다음과 같이 두고,

$$
\mathbf{z}\_k =  \begin{bmatrix} P\_1^k \ P\_2^k \ \vdots \ P\_N^k \end{bmatrix}
$$

이를

$$
\mathbf{z}\_k = h(\mathbf{x}\_k) + \mathbf{v}\_k
$$

로 표현할 수 있다. 비선형 함수 $h(\cdot)$는 각 위성까지의 거리와 시계 오차 항을 반영하는 역할을 하며, $\mathbf{v}\_k$는 측정 잡음(measurement noise)다.

**EKF 선형화 과정**

EKF는 칼만 필터의 선형 모델 가정이 깨지는 비선형 문제에서, 측정식과 상태천이를 1계 테일러 전개로 선형 근사하여 적용한다.

1. **예측 단계 (Prediction)**
   * 상태 예측:

     $$
     \hat{\mathbf{x}}*{k \mid k-1} = \mathbf{F}*{k-1} ,\hat{\mathbf{x}}\_{k-1 \mid k-1}
     $$
   * 공분산 예측:

     $$
     \mathbf{P}*{k \mid k-1} = \mathbf{F}*{k-1} ,\mathbf{P}*{k-1 \mid k-1},\mathbf{F}*{k-1}^\mathsf{T} + \mathbf{Q}\_{k-1}
     $$

     여기서 $\mathbf{Q}\_{k-1}$는 시스템 노이즈 공분산.
2. **갱신 단계 (Update)**
   * 야코비 계산:

     $$
     \mathbf{H}*{k}  = \frac{\partial h(\mathbf{x})}{\partial \mathbf{x}} \Biggl\rvert*{\mathbf{x}=\hat{\mathbf{x}}\_{k \mid k-1}}
     $$
   * 칼만 이득(Kalman Gain):

     $$
     \mathbf{K}*k  = \mathbf{P}*{k \mid k-1},\mathbf{H}\_k^\mathsf{T}  \Bigl\[\mathbf{H}*k,\mathbf{P}*{k \mid k-1},\mathbf{H}\_k^\mathsf{T} + \mathbf{R}\_k\Bigr]^{-1}
     $$
   * 상태 갱신:

     $$
     \hat{\mathbf{x}}*{k \mid k}  = \hat{\mathbf{x}}*{k \mid k-1}  + \mathbf{K}\_k   \Bigl(\mathbf{z}*k - h(\hat{\mathbf{x}}*{k \mid k-1})\Bigr)
     $$
   * 공분산 갱신:

     $$
     \mathbf{P}\_{k \mid k}  = \Bigl(\mathbf{I} - \mathbf{K}\_k,\mathbf{H}*k\Bigr),\mathbf{P}*{k \mid k-1}
     $$

이 과정을 통해, 직전 시점 정보(예측치)와 현재 GNSS 측정치를 종합하여 최적의 위치, 속도, 시계 관련 추정치를 구한다.

**장점**

* **시간연속성 활용**: 독립적인 에포크별 위치 추정에 비해 잡음이나 일시적인 위성 관측 이상치에 대한 강인성이 높아진다.
* **속도 정보 추정**: EKF 과정에서 속도(및 가속도까지 확장 가능)를 명시적으로 모델링하므로, 동적 정확도를 높일 수 있다.
* **추가 센서 융합 용이**: 관성측정장치(IMU), 전자나침반 등과의 센서 융합(Integration)에 EKF가 자주 사용된다.

**유의사항**

* **모델링 정확도**: 잘못된 동역학 모델(가속도 예측 등)을 사용하면 오히려 필터 오차가 누적될 수 있다.
* **초기 공분산 설정**: 초기 추정치의 신뢰도(공분산)에 따라 필터의 수렴 속도와 안정성이 좌우된다.
* **선형화 오차**: 비선형성이 큰 환경(예: 큰 기동, 과도한 속도 변동)에서는 EKF의 선형 근사가 충분치 않을 수 있다. 이 경우 UKF(Unscented Kalman Filter) 같은 대안을 고려하기도 한다.

#### RAIM (Receiver Autonomous Integrity Monitoring)

GNSS 수신기는 보통 여러 위성으로부터 측정된 의사거리 정보를 가지고 사용자 위치를 추정한다. 그러나 일부 위성 측정값에 이상치(outlier)가 존재할 경우, 단순 최소제곱이나 EKF만으로는 이를 검출하고 제거하기 어려울 수 있다. RAIM은 이러한 GNSS 측정값의 무결성(Integrity)을 사용자가 스스로 모니터링하는 기법으로, 삼변측량 기본 알고리즘과 연계되어 동작한다.

**RAIM 기본 아이디어**

* **초과측정(Redundancy)**: 3차원 측위와 시계 바이어스 추정에는 최소 4개 위성이 필요하다. 하지만 보통은 5개 이상의 위성을 관측 가능하므로 초과측정이 발생한다.
* **잔차 기반 판단**: 측정방정식을 풀어 낸 뒤, 남는 잔차(residual)가 이상적으로는 오차통계와 부합해야 한다. 특정 위성 측정치가 크게 틀리면 잔차 분포가 비정상적으로 커진다.
* **위성 제외 테스트**: RAIM은 가능한 한 위성 조합별로 (혹은 부분집합별로) 잔차 검사(Chi-square test, 혹은 다른 통계적 방법)를 수행해, 문제를 일으킨 것으로 추정되는 위성 측정값을 식별·배제한다.

**잔차 기반 RAIM 수식 예시**

단순화하여, 5개의 측정값(위성 1\~5)이 있다고 하자. 3차원 위치와 시계 바이어스($\delta t$) 추정에는 4개의 방정식이 필요하므로, 5개의 방정식을 최소제곱으로 풀면 1차원 여유도가 생긴다(자유도: $5 - 4 = 1$). 이때 각 에포크별로 최소제곱 해 $\hat{\mathbf{u}}$를 구한 후,

$$
\mathbf{r} = \mathbf{z} - \mathbf{H},\hat{\mathbf{u}}
$$

로 계산되는 잔차 벡터 $\mathbf{r}$를 이용해 통계량(예: $\chi^2$ 분포 기반)을 구한다:

$$
rT = \mathbf{r}^\mathsf{T} \mathbf{r}
$$

(가중 행렬 $\mathbf{W}$가 있을 경우, $\mathbf{r}^\mathsf{T}\mathbf{W}\mathbf{r}$ 형태가 됨) 이 $T$가 특정 임계값(Threshold)보다 크면 이상치가 있다고 판단하여, RAIM 절차에 따라 의심 위성을 찾고 제외하거나 경고를 발생시킨다.

**위성 제외 절차 예시**

1. 모든 위성(5개)을 사용해 해를 구한다.
2. 잔차 $\mathbf{r}$를 이용해 통계량 $T$를 계산한다.
3. 임계값을 초과하면, 각 위성을 하나씩 제외한 뒤(즉, 4개 위성만 사용하는 경우 5가지) 각각 새로운 해를 구하고, 잔차를 재검토한다.
4. 특정 위성을 제외했을 때 잔차가 정상 범위로 돌아오면, 그 위성이 이상값을 낸 것이라고 추정할 수 있다.

**RAIM 한계**

* **초과측정 필요**: 적어도 5개 이상의 위성을 안정적으로 관측 가능해야 한다.
* **검출 대 기각 사이의 Trade-off**: 너무 엄격하게 설정하면 정상 측정도 제외될 수 있고, 너무 느슨하면 진짜 이상치가 검출되지 않을 수 있다.
* **수치적 복잡도**: 많은 위성을 관측할수록(예: 8기 이상), 모든 부분집합을 검사하는 것은 계산량이 커진다.

#### 다중 GNSS 활용

과거에는 GPS 단일 위성군만 사용했으나, 현재는 GLONASS, Galileo, BeiDou 등 여러 GNSS를 동시에 활용할 수 있다. 이렇게 다중 위성을 병합하여 측위를 수행하면,

* 보이는 위성 수가 대폭 증가한다(초과측정 증가).
* RAIM이나 DOP 개선 효과가 커진다.
* 다양한 주파수대( L1, L2, L5 등 )를 이용하면 전리층 및 다중경로 오차 보정이 향상된다.

#### 정밀측위(추가 기법)

기본 삼변측량 알고리즘은 의사거리 기반으로 위치를 구하지만, 오차 소스를 더 정교하게 보정하거나 위성-수신기 간 위상(Phase) 측정을 사용하면 cm\~mm 급 정밀도를 얻을 수 있다. 예를 들어,

* **RTK (Real-Time Kinematic)**: 기준국과 이동국 간 위상 관측치 차이(Differential)로 정밀도를 확보.
* **PPP (Precise Point Positioning)**: 광범위한 오차 모델(정밀 위성 궤도·시계 오차 등)을 적용하여 단독수신기 상태에서도 cm 급 정밀도 가능.

다만 이러한 정밀측위 기법은 기본 의사거리 측정 기반 삼변측량에서 더 세분화된 관측량(Phase Range)을 추가로 활용하고, 다양한 보정 매개변수를 엄밀하게 추정한다는 점에서 부차적인 알고리즘적 복잡도가 높다.
