Accuracy Improvement of DEM Using Ground Coordinates Package

공공삼각점 위치자료를 이용한 DEM의 위치 정확도 향상

Abstract

In order to correct the provided RPC and DEM generated from the high-resolution satellite images, the acquisition of the ground control point (GCP) must be preceded. This task is a very complicate that requires field surveys, GPS surveying, and image coordinate reading corresponding to GCPs. In addition, since it is difficult to set up and measure a GCP in areas where access is difficult or impossible (tidal flats, polar regions, volcanic regions, etc.), an alternative method is needed. In this paper, we propose a 3D surface matching technique using only the established ground coordinate package, avoiding the ground-image-location survey of the GCP to correct the DEM produced from WorldView-2 satellite images and the provided RPCs. The location data of the public control points were obtained from the National Geographic Information Institute website, and the DEM was corrected by performing 3D surface matching with this package. The accuracy of 3-axis translation and rotation obtained by the matching was evaluated using pre-measured GPS checkpoints. As a result, it was possible to obtain results within 2 m in the plane location and 1 m in height.

고해상도 스테레오 위성영상에서 제공되는 초기 RPC의 보정과 이로부터 생성된 DEM 보정을 위해서는 지상기준점 획득이 선행되어야 한다. 기준점 획득은 현장답사와 GPS 측량, 지상기준점에 대응하는 영상좌표 독취 등을 수행해야 하는 매우 번거로운 작업이다. 그리고 접근이 힘들거나 불가능한 지역(갯벌, 극지대, 화산지대 등)에는 기준점 설치와 측정이 힘들기 때문에 이를 대체할 수 있는 방법이 필요하다. 이 연구에서 WorldView-2 위성영상과 초기 RPC로부터 제작된 DEM 보정을 위해 기준점의 지상좌표와 영상좌표 측량을 피하고 기 구축된 지상좌표 패키지만을 활용한 3차원 표면매칭 기법을 제안한다. 국토지리정보원 국토정보플랫폼에서 공공삼각점의 위치정보를 획득하고 제작된 DEM과 3차원 표면매칭을 실시하여 DEM을 보정하였다. 3차원 표면매칭으로 얻은 3축 이동과 회전량의 정확도를 평가하기 위하여 검사점으로 확인한 결과 평면오차는 2 m, 수직오차는 1 m 이내의 결과를 얻을 수 있었다.

1. 서론

고해상도 위성영상을 사용하여 고품질의 DEM을 확보하기 위해서는 정확한 센서정보 (외부표정요소 또는 RPC(Rational Polynomial Coefficient))가 요구된다. RPC 는 물리적 센서모델인 외부표정요소에 비해 3차원 위치 측정 방법이 간편하고 정확도 또한 떨어지지 않기 때문에 많이 사용되고 있다(Fraser and Hanley, 2005; Tao and Hu, 2002). 현재 운용중인 1 m 이하 급 대부분의 고해상 도 위성영상에서 제공하는 RPC는 해상도에 준하는 위치결정 정확도를 가지고 있지 않다. 이중 WordView-2 위성영상에서 제공되는 RPC는 평면위치 오차(CE90; Circular error at 90% probability)가 약 3 m이다(USGS, 2021). RPC 보정을 위해서는 최소 한 점 이상의 지상기준점이 요구된다. 그러나 지상기준점 측량은 영상에서 독취작업, 기준점의 현장답사와 GPS 측량 등을 수행해 야 하는 매우 번거로운 작업이다. 특히 접근이 곤란하거나 불가능한 지역의 위성영상으로 DEM을 제작하기 위한 기준점 획득은 더욱 힘들다. 기존에는 기준점 대안으로 DEM, 수치지도와 같은 기 구축된 기준자료를 활용하는 연구가 수행되었다(Oh and Jung, 2012; Oh and Lee, 2015). 이들은 영상 전반에 기준점 확보가 가능할 뿐만 아니라 현장 방문에 의한 실측을 하지 않아도 되기 때문에 RPC 정확도 향상을 위한 작업의 효율성과 편리성을 확보할 수 있고, 자동화가 가능하다는 장점이 있다. 그러나 기존 SRTM 또는 TanDEM-X DEM은 고 해상도 위성영상으로 제작되는 DEM에 비해 해상도와 상대정확도가 떨어지며, 수치지도는 접근불능지역이 나 국외지역에서는 확보가 어렵다.

몇몇 연구자들은 고해상도 스테레오 위성영상과 제공되는 RPC로부터 DEM 또는 영상의 기하학적 정확도 향상을 위해, 상대표정 번들조정법을 제안하였다(Toutin et al., 2012; Zhang et al., 2015; Yang et al., 2017). 이 연구들은 주로 접합점의 영상좌표와 RPC로부터 계산된 접합 점의 가상 지상좌표를 이용하는 방식이므로 DEM의 상대적 위치 정확도 향상에만 제한된다.

위치결정 정확도 향상을 위한 또 다른 접근방식은 DEM 매칭이다. 이 방법은 기 확보한 DEM에 위성 초기 센서모델(또는 RPC)로부터 제작된 DEM을 매칭시켜 3축 이동과 회전량을 구하는 것이다. DEM을 일치시키기 위해 가장 많이 활용되는 방법이 ICP(Iterative Closest Point) 방법이다(Bessel and Mckay, 1992). ICP는 두 DEM 에서 임의점들 간 최소거리 조건으로 매칭 쌍을 찾은 다음 두 DEM 간 변환(위치이동, 회전)계수를 구하는 것이다. 그러나 이 방법은 일정 패턴의 밋밋한 형태를 가진 지형이거나 일부분만 중복된 DEM에서는 매칭 왜곡이 발생할 수 있으며, 같은 점을 찾기 위한 반복계산 시간이 많이 걸린다는 단점이 있다. 유사 ICP 중 하나인 LHD(least-squares height-difference) 기법은 같은 평면 위치에서 높이 차가 최소가 되도록 반복계산 하는 방법이다(Rosenholm and Torlegård, 1988; Ebner and Strunz, 1988). Kim and Jeong (2011)은 이 방법을 확장하여 푸 쉬브룸 위성영상 센서모델링에 적용하였다. 그러나 이 방법 또한 기본적으로 좋은 품질의 기준 DEM이 요구되며, DEM경사가 작을수록 매칭의 정확도는 떨어진다. 이 연구는 기준 DEM을 사용하지 않고 기존에 측량해둔 높은 정확도의 기준점을 활용하는 3D 표면매칭 기법을 제안한다. 이 매칭방법을 적용하면 기준점들의 영상독취 작업은 필요로 하지 않기 때문에 DEM 보정을 위한 전체 프로세싱을 자동화할 수 있다. 국토지리정 보원에서는 전국을 대상으로 국가기준점을 설치하고 그 성과를 관리해 오고 있다. 대표적인 국가기준점으로는 삼각점(약 16,000여점)이 있고, 2008년부터 통합기준점(평면위치, 표고, 타원체고, 중력값 등이 정밀 측정되어 있는 점)을 설치하기 시작하여 현재까지 전국적으로 모두 5,500여점이 설치되어 있다(NGII, 2021). 이 점들 의 위치정보는 국토지리정보원 국토정보플랫폼을 통해 무상으로 사용자들에게 제공하고 있다.

Lee and Michael (2019)은 삼각점과 통합기준점, 포인트–표면 매칭기법을 적용하여 KOMPSAT-3 위성영상 제공 RPC의 정확도를 향상시켰다. 이 연구에서 삼각점은 대부분 산 정상에 위치하여 나무들과 숲에 둘러싸여 제작 DEM과의 매칭에서 오차원인이 되고, 탁 트인 공간에 설치된 통합기준점이라 하더라도 위성영상 범위 내에 확보할 수 있는 개수가 적은 문제가 있다. 최근 국토정보플랫폼에서는 주로 시가지 도로주변에 설치된 공공삼각점의 위치정보를 제공하기 시작하였다(NGII, 2021). 공공삼각점의 경우, 주변 장애물이 가리지 않기 때문에 위성영상으로부터 제작된 DEM과의 매칭오차 는 기존 삼각점보다는 작을 가능성이 크며, 통합기준점보다는 많은 수를 확보할 수 있다는 장점이 있는 반면 점들이 도심지에만 몰려 있다는 단점도 있다. 한편, Gruen and Akca(2005)에 의해 제안된 3D 표면매칭 기법은 표면 대표면, 표면 대 곡선, 표면 대 저밀도 점군 매칭에 모두 적용이 가능하다. 따라서 이 연구에서는 3D 표면매칭 기법을 이용하여 WordView-2 위성영상의 초기 RPC로 제작한 DEM 보정을 위한 공공삼각점의 활용 가능성을 테스트하고, 그 정확도를 제시하고자 한다.

2. 3차원 표면매칭

DEM 보정을 위한 연구 흐름도는 Fig. 1의 좌측그림 과 같다. 제공된 초기 RPC로부터 제작된 DEM은 공공 삼각점(Fig. 1에서 Public Control Points, PCPs)을 이용하여 3D 표면매칭으로 DEM과 GCPs 간의 3차원 관계식 즉, 위치이동과 회전 변위량을 구하고 이로부터 DEM 을 보정한다. Fig. 1의 우측 그림은 제작된 DEM을 공공 삼각점 기준으로 매칭하기 위한 개념을 나타낸 것이다.

먼저 확보한 공공삼각점의 위치가 포함된 표면을 f(X, Y, Z)라 가정하고 이와 동일한 DEM 영역의 표면을 g(X, Y, Z)라 하면 두 표면 간의 기하학적 관계는 이상적으로는 동일한 상태가 아니기 때문에 Eq. (1)과 같이 두 표면 사이에 오차벡터 v(X, Y, Z)가 추가된다. 한편 3D표면매칭의 목적은 Eq. (2)을 이용하여 두 표면 사이의 거리를 최소화하기 위한 3축 이동과 회전 변환계수를 구하는 것이다. 참고로 위성영상에서 제공되는 RPC로 제작된 DEM은 실제 기준점 지상좌표들과의 스케일 차 이가 거의 없으므로 이 연구에서는 생략하였다.

\(f_{i}(X, Y, Z)-v_{i}(X, Y, Z)=g_{i}(X, Y, Z)\)       (1)

where i = i^th point

\([X Y Z]_{i}^{T}=M\left[X_{o} Y_{o} Z_{o}\right]_{i}^{T}+t\)       (2)

where M=M(ω, φ, κ) is the 3D rotation matrix; t=[tX tY tZ] T is a 3D translation vector; [XYZ] is the location of the public control point, and [Xo Yo Zo] is the location of the conjugate approximate surface gi o (X, Y, Z) of the DEM.

Eq. (2) 변환계수를 구하기 위해, Eq. (1)에 최소제곱 법 적용을 위한 테일러 선형화를 실시하면 Eq. (3)과 같이 나타낼 수 있다.

\(\begin{gathered} f(X, Y, Z)-\nu(X, Y, Z)=g^{\circ}(X, Y, Z)+ \\ \frac{\partial g^{\theta}(X, Y, Z)}{\partial X} d X+\frac{\partial g^{Q}(X, Y, Z)}{\partial Y} d Y+ \\ \frac{\partial g^{\circ}(X, Y, Z)}{\partial Z} d Z \end{gathered}\)       (3)

where omitted the index i for simplicity

Eq. (3)에서, g^o (X, Y, Z) = aX + bY + cZ + D = 0과 같은 평면방정식으로 표현될 수 있으며, 이로부터 Eq. (3)의 \(\frac{\partial g^{9}(X, Y, Z)}{\partial X} \frac{\partial g^{\prime}(X, Y, Z)}{\partial Y} \frac{\partial g^{\circ}(X, Y, Z)}{\partial X}\)은 DEM의 법선 벡터를 이용하여 Eq. (4)로 표현할 수 있다. 그리고 dX, dY, dZ는 Eq. (5)와 같이 나타낼 수 있다.

Fig. 1. Flow chart for 3D surface matching and correction of the produced DEM. Illustration of the 3D surface matching between PCPs and DEM produced by the provide RPCs

\(\begin{aligned} &\frac{\partial g^{\circ}(X, Y, Z)}{\partial X}=\frac{a}{\sqrt{a^{2}+b^{2}+c^{2}}}=g_{x}, \\ &\frac{\partial g^{\circ}(X, Y, Z)}{\partial Y}=\frac{b}{\sqrt{a^{2}+b^{2}+c^{2}}}=g_{V}, \\ &\frac{\partial g^{\theta}(X, Y, Z)}{\partial Z}=\frac{c}{\sqrt{a^{2}+b^{2}+c^{2}}}=g_{z} \end{aligned}\)       (4)

\(d X=\frac{\partial X}{\partial p_{k}} d p_{k}, \quad \partial Y=\frac{\partial Y}{\partial p_{k}} d p_{k}, d Z=\frac{\partial Z}{\partial p_{k}} d p_{k}\)       (5)

where p_k∈{t_X, t_Y, t_Z, ω, φ, κ} is the k^th transformation parameter in Eq. (2).

\(\begin{aligned} -v(X, Y, Z)=& g_{X} d t_{X}+g_{Y} d t_{Y}+g_{z} d t_{z} \\ &+\left(g_{X} d_{11}+g_{Y} \mathrm{~d}_{21}+g_{z} d_{31}\right) d \omega \\ &+\left(g_{X} d_{12}+g_{Y} d_{22}+g_{z} d_{32}\right) d \varphi \\ &+\left(g_{X} d_{13}+g_{Y} d_{23}+g_{Z} d_{33}\right) d \kappa-L \end{aligned}\)       (6)

where, L=f(X, Y, Z) – go (X, Y, Z) is the normal distances between the points and correspondentsurface elements of the DEM. Surface points (g^o (x, y, z)) is calculated by bi-linear interpolation to determine L. the coefficient terms d_jj are computed by the partial derivatives with respect to the three rotation angles (ω, φ, κ ) in Eq. (5).

Eq. (3)에서 구하고자 하는 미지 변수 항을 Δ, 편미분 항을 B, DEM 표면과 기준점간의 수직거리들로 이루어진 벡터항을 L이라 하면 Eq. (7)과 같이 요약할 수 있고, 이에 대해 최소제곱법을 적용하면 Eq. (8)로 정리할 수 있다.

\(-v=B \Delta-L\)       (7)

\(\Delta=\left(B^{T} W B\right)^{-1}\left(B^{T} W L\right)\)       (8)

where Δ^T = [dtX, dtY, dtZ, dω, dφ, dκ], W = weight

여기서, 3D 표면매칭은 기본적으로 다수의 기준점과 DEM간 동일 평면위치에서의 거리를 이용하는 방법이므로 모든 기준점들과 이에 대응하는 DEM 표면들 간의 거리(di )는 어느 정도 일정해야 한다. 따라서 임의의 기준점에 대응하는 DEM 표면의 과대오차 지점을 제거하기 위해 다음과 같은 조건(W)으로 매칭에 참여하는 점을 제한하였다.

\(W= \begin{cases}1 & \text { if }\left|d_{i}-\operatorname{mean}\left(d_{i}\right)\right|<1.0 \sigma \\ 0 & \text { else }\end{cases}\)       (9)

where σ is the standard deviation of d_i.

Eq. (2)에서 변환계수는 Eq. (8)의 반복계산에 의해 Δ 항이 0에 가까울 때 결정된다. 따라서 매번 반복계산 마다 변환계수는 t_X¹ = t_X⁰ + dt_X, …와 같이 업데이트되므로 제안방법의 수렴여부는 참여하는 기준점과 DEM의 품질에 따라 달라진다.

3. 적용실험과 결과분석

사용된 인공위성 영상자료는 Fig. 2에서 보는 바와 같이 전남 순천시 일대(지상면적 약 20 km × 17 km)를 대상 으로2010년2월에관측수집된WorldView-2Panchromatic 입체영상이다(Table 1). 촬영된 영상의 왼쪽 대부분은 구름이 차지하고 있으므로 DEM 제작은 Fig. 3의 영상에서 박스 부분만을 대상으로 하였다. DEM은 입체영상과 제공 RPC를 ERDAS Imagine Photogrammetry 소프트웨어에 입력하여 제작하였다(Fig. 3에서 오른쪽 그림). DEM보정을 위한 기준점으로 국토정보 플랫폼 웹사이트에서 공공삼각점 32점을 획득하였다(Fig. 3에서 빨간색 점들). 테스트 지역의 공공삼각점들은 대부분 2017년에서 2018년 사이에 설치된 점들이므로 위성영상 촬 영시기와는 다소 차이가 있으나 대부분 고도가 낮은 시가지의 도로 교차지점에 위치하고 있으므로 지형변화가 거의 없다. 따라서 2010년 획득 위성영상으로 제작된 DEM 보정용 기준점으로 사용하기에는 문제가 없다고 볼 수 있다. 결국 기준점과 동일한 평면위치의 DEM 높이 값을 찾고, Eq. (5)을 적용한 결과 최종적으로 15점만 선정되어 매칭에 참여할 수 있었다. 본 기법의 정확도 평가를 위해 위성영상 촬영과 같은 시기인 2010년 7월에 GPS 측량으로 20개 검사점을 획득하였으며(Fig. 3에서 초록색 점들), 이에 대응하는 영상좌표도 부화소 단위로 독취 하였다.

Table 2는 선정된 기준점으로 매칭 후 구한 변환계수를 나타낸 것이다. 여기서, 3축 이동량은 3 m 이내였으며, 회전변환계수는 1분 내외로 미소하였다. WorldView-2스테레오 위성영상의 평면위치와 높이에 대한 공칭오차(CE90, LE90), 평면제곱근오차가 각각 5 m, 2 m 내외인 것을 감안했을 때(USGS, 2021), 매칭 후 3축 이동변 환계수 또한 비슷한 결과를 보여준다.

Fig. 2. WorldView-2 image pair for the Suncheon, Korea.

Table 1. Specifications of tested satellite images

Fig. 3. PCPs and GPS check points location on the image, and the produced DEM based on the provided RPCs.

Table 2. The computed parameters of the produced DEM through 3D surface matching

Fig. 4. Distances between the PCPs and the surface points of the DEM before and after matching, plane distance (1st row), 3D distances (2nd row).

Table 3. Distance errors between PCP locations and DEM surfaces before and after matching

Fig. 4는 구한 변환계수의 신뢰도를 평가하기 위해, PCP와 공간상 가장 가까운 거리에 있는 DEM 표면 지점들 간의 평면상 거리차와 공간상 거리차를 시각적으로 표현한 것으로써 매칭 전에 비해 매칭 후 DEM 오차가 크게 줄어 들었음을 확인할 수 있다. Table 3은 Fig. 4 (2nd row)의 공간상 거리차를 나타낸 것이다. 여기서, RMSE 거리오차와 최대 거리오차는 보정전 3.3 m와 4.0 m에서 보정후 0.1 m와 0.2 m로 크게 향상된 것을 볼 수 있다.

그러나 Table 3의 결과는 PCP와 DEM 표면간 매칭으로 구한 변환계수로 재계산하여 비교한 것이므로 매칭에 대한 절대적인 평가가 될 수 없다. 따라서 검사점들을 이용하여 정확도를 평가하였다. 먼저, 측정한 GPS의 영상좌표와 제공되는 초기 RPC 계수를 이용하여 지상 좌표를 계산하였으며, 측정한 GPS의 지상좌표와 비교 하였다(Fig. 5의 좌측그림, Table 4). 그리고 초기 RPC로 구한 지상좌표를 매칭 후 구한 변환계수로 보정하였다 (Fig. 5의 우측그림, Table 4). Fig. 5는 GPS 포인트들의 보정 전, 후에 대한 평면위치와 공간상 거리오차를 각각 시각적으로 표현한 것이다. 그 결과, 보정 전에 비해 보정 후의 정확도가 다소 향상된 것을 볼 수 있다. Table 4에서도 보는 바와 같이 매칭전에 비해 매칭후의 평면오차는 2 m 이내, 수직오차는 3 m 이내로 감소하였다. 결과적으로 공공삼각점과 3D 표면매칭으로 초기 DEM 을 보정할 수 있는 가능성을 제시할 수 있었다.

Fig. 5. Deviations in GPS checkpoints. Displayed are horizontal (1st row) and spatial differences (2nd row) between measured GPS coordinates (triangle points) and coordinates reconstructed from the image pair (circle points).

Table 4. Errors of GPS points before and after matching

4. 결론

이 연구를 통해, 3D 표면매칭 방법과 기존에 확보한 기준점을 이용하여 스테레오 위성영상과 제공되는 RPC로부터 제작된 DEM을 2 m 정도의 정확도로 보정할 수 있었다. 제안방법은 나아가 스테레오 영상의 접합점을 이용하여 제공되는 RPC 계수 또한 보정이 가능할 것으로 생각한다. 이 방법의 특징은 DEM의 정확도를 향상시키기 위한 기준점의 지상좌표만 필요하고 영상좌표 독취가 필요 없으며, 기존의 물리적 센서모델링과 같은 복잡한 과정이 필요치 않다. 또한 DEM 매칭을 위한 기준 DEM 확보 또한 필요치 않다는 장점이 있다. 아울러 이 방법은 각 나라별 국가기준점이나 기 확보한 GPS 또는 라이다 데이터 등을 이용하여 스테레오 위성 영상과 그 초기 메타자료로부터 제작된 DEM을 보정하는데 범용적으로 활용될 수 있을 것이다.