1. 서론
최근 TerraSAR-X, Cosmo-SkyMed, Kompsat-5 등 사용가능한 고해상도 SAR 위성 영상의 수가 늘어남에 따라 다양한 분야에서 이를 활용한 연구가 진행중에 있다. 특히, 군사, 재해, 재난 분야에서 변화탐지 등의 다양한 적용 분야에 대한 요구가 커지고 있으며, 이를 위해서는 다양한 환경에서의 SAR영상간 정밀 영상 정합 연구가 필수적으로 필요하다. 그러나 SAR 영상은 다량의 스펙클 노이즈와 영상 촬영 각도, 모드 등 다양한 촬영 환경에 따라 영상밝기 값이 변화하는 특성이 있기 때문에, 자동화된 정밀정합 기법의 적용에 어려움이 있다 (Sansoti et al., 2006). 다중관측각 환경에서의 SAR영상은 고층빌딩이 위치하는 도심지와 산악지형에서 촬영각도에 따른 변화가 크게 나타나는지역이거나, 군사, 보안목적의 이유로 지상정밀 DEM(Digital Elevation Model), DSM(Digital Surface Model)을 취득하지 못하는 경우 이러한 기하학적 왜곡의 보정은 더욱 힘든일이다(Thiele et al., 2007). 따라서 SAR영상의 기하학적 변화와 스펙클노이즈에 강건한 자동영상정합기술이 필요하다.
대표적인 정합기법으로는 크게 특징 기반 기법과 밝기값 기반 기법이 존재한다. 특징 기반 기법은 주로 코너, 엣지 등과 같은 이미지내에서 크게 두드러지는 특징점을 찾아 서술자를 추출하여 매칭하는 기법이고 밝기값 기반 기법은 두 이미지 밝기값을 통계적 비교를 통해 매정합을 수행하는 기법이다.
대표적인 특징점 기반 기법으로는 SIFT(Scale Invariant Feature Transform) 알고리즘이 존재한다(Lowe2004). SIFT알고리즘은 컴퓨터 비전분야에서 주로 활용되어지는 기법이지만 위성영상 정합에서도 다양하게 활용되어 지고있다(Goncalves et al., 2011). SIFT 알고리즘은 크게 특징점 추출, 특징점의 서술자 추출, 특징점 매칭 과정을 통한 3단계로 나눠진다. 이러한 전통적인 SIFT 알고리즘을 위성영상간 정합에 사용하고자 개선하려는 연구들이 지속되어 왔다. Schiwind et al. (2010)은 SIFT 알고리즘의 첫번째 스케일에서 추출되는 키포인트는 비재현성이 가장 높기 때문에 이를 배제하는 SIFT-OCT 기법을 제안하였다. Fanet al. (2013)은 기준영상과 대상영상 사이에 회전에 대한 영향이 거의 없다면 서술자 추출단계의 방향성 분지정단계(Orientation Assignment)를 생략해도 좋음을 보였다. 그 외에도 Wang et al. (2012)는 Bilateral filter를 사용한BFSIFT(Bilateral Filter SIFT)를 제안하였고, Wang et al. (2014)는 가우시안 필터를 개선하여적용한 AAG-SIFT(Adapted anistropic Gaussian SIFT)를 제안했다. SongandZhang(2010)은 SIFT알고리즘의 속도를 개선하기 위해 특징점 추출단계를 간소화한 SURF 알고리즘을 정합에 사용했다. Delinger et al. (2015)는 SAR 영상간정합을 위하여 SIFT 알고리즘의 차분 기울기 연산자를 영역비 기울기 연산자로 대체하여 스펙클노이즈를 줄인 SAR-SIFT 알고리즘을 제안했다. JeonandKim(2018)은 SAR-SIFT와 DLSS(DenseLocal Self-Similarity)를 결합하여 광학영상과SAR영상간 정합기법을 제안하였다. 그 외에도 다양한 SAR-SIFT기반 알고리즘이 연구되어 왔지만, 여전히 SAR-SIFT는 영상사이즈가 커짐에 따라 연산량이 많아지고 추출되는 특징점의 수가 많아지면서 오경보율이 높아지는 단점이 존재한다.
밝기값 기반 정합기법에 대한 연구도 진행되어 왔다. MI(Mutual Information)와NCC(Normalized Cross Corre - lation) 기법은 밝기값 기반 정합기법의 대표적인 방법으로 원격탐사분야에서 사용되어져 왔다(Chen et al., 2003). Suriand Reinartz(2010)은 도심지역에서의 MI를 이용한 SAR위성과 광학위성 영상정합을 수행했다. MI 기법은 비선형 영상에 대해 정합이 가능하지만 처리시간이 오래 걸린다는 단점이 존재한다. Xiang et al. (2018) 은 엣지를 강조해주는PC(Phase Congruency) 영상에 (Kovesi, 2000) SAR영상의 스펙클노이즈에 강건하게 만들어주는 2-DGabor filter를 적용 후 NCC기법을 적용하여 중국 Gaofen-3 위성영상의 다중 모드 간 정합에 사용하였다. 개략정합 후 정밀정합을 수행하는 방식으로, 정밀정합 단계에서 PC영상을 생성하여 상호 상관계수를 통한 정합방법을 제안하였다.
그러나 이러한 기존 정합방법들은 변환식을 통해 정합을 수행하는 강성정합(Rigid)구조이기 때문에 지형효과에 의해 왜곡이 발생하는 다중관측각 영상에는 적합하지 않다. 조밀 비강성 정합(DenseNon-rigid) 모델은 모든 픽셀에 대하여 변화량을 찾기 때문에 변환식을 사용하지 않고 정합을 수행할 수 있다. 대표적인 조밀 비강성 정합 모델 중 하나인 광학흐름(Opticalflow)는 주로 컴퓨터 비전분야에서 객체를 추적하기 위하여 사용된다. Plyer et al. (2015)는 기존Folki(FlowOpticalFlow Lucas-KanadeIterative) 기법(Lucas and Kanade, 1981)을 위성영상에적합하도록 개선하여 다중모드 SAR영상 간 정합과 간섭도 생성을 위한 상호정합을 수행하였다. Brigot et al. (2016)은 Folki알고리즘을 개선한 Gefolki (Geo science Extended Folki) 알고리즘을 제안하여 SAR 영상과 Lidar(Light Detection and Ranging) 영상 간 정합에 성공적으로 사용했다. 개략정합이 수행 되어져야 하는 선행조건이 있지만 Gefolki 알고리즘은 해당지역의 정밀DEM, DSM 이 없이 정합이 수행이 가능하다는 장점이 있다.
본 논문에서는 다중 관측각 영상에 대해 coarsetofine 단계의 정합기법을 수행한다. 개략 정합 단계는 Xiang et al. (2018)에서 제안한 기법에 착안하여, 기준 영상과 대상 영상 사이즈를 1000 픽셀이 하가 되도록 다운 샘플링 후 SAR-SIFT를 수행한다. 이를 통하여 연산량을 줄일 수 있고, SAR 영상에서 나타나는 기하학적 왜곡을 효과적으로 제어할 수 있다. 이때, SAR-SIFT 단계에서 추출한매칭점의 개수가 10개 이하일 경우 샘플링 횟수를 늘려 반복한다.
정밀정합 단계에서는 성능비교를 위하여 NCC 기법과 PC맵에 NCC기법을 적용한 PC-NCC기법, MI기법, Gefolki기법을 사용하여 비교 분석한다. NCC, PC-NCC, MI기법은 기준영상과 개략 정합된 영상간에 패치별로 오프셋을 찾아 정합을 수행하는 방법으로 소요시간을 고려하여 적절하게 패치 사이즈와 window 사이즈를 조절하여 수행한다. 이때, PCmap은 엣지를 강조한다고 알려져 있는maximum moment 값을 사용한다(Kovesi, 2000). NCC, PC-NCC, MI기법은 개략정합 모델과 정밀정합단계의 변환 모델을 합쳐 최종변환 모델을 추출한다. Gefolki는 조밀비 강성정합 모델로 모든 픽셀에서 오프셋을 추정하기 때문에, 변환모델 없이 개략정합된 이미지에서 정밀정합된 최종 이미지를반환한다.
최종적으로 정합기법의 비교분석을 모자이크 이미지를 통한 정성적 분석과, 각 영상에서두드러지는 점에서 취득한 포인트에 대한 RMSE(RootMeanSquareError) 와 이미지 퀄리티 평가지표중 하나인 FSIM(Feature Similarity) 지수를 사용하여 정량적 분석을 수행한다.
2. 연구 지역 및 데이터
연구 지역은 고해상도 SAR 영상의 특성인 기하학적 왜곡을 관측할 수 있도록 수계, 도심지, 산악지형을 포함하는 대한민국 인천 석남동 지역을선정했다(Fig. 1).
Fig. 1. Study area and the coverages of TerraSAR-X image: 2018/04/03 (Yellow), 2018/04/08 (Blue), 2018/04/13 (Red).
1) TerraSAR-X
TerraSAR-X는 독일 항공우주국에서 개발된 X밴드 SAR 위성으로 2007년 6월 발사되어 현재까지 운영 중이다. 촬영 모드로는 Wide ScanSAR모드, StripMap모드 Spotlight모드 High Resolution Spotlight모드, Staring Spotlight모드가 있으며, 11일의 재방문 주기를 갖는다. 본 연구에서는 다중관측각 영상간 정합 성능비교를 위하여 고해상도 모드인 TerraSAR-XStaring Spotlight모드 영상 3장을 사용했다. 취득한 영상의 비행방향 해상도는 약 16 cm이며 거리 방향 해상도는 50 cm에서 1 m까지 관측각에 따라 상이하다(Table1). 영상은 2018년 4월 3일, 2018년 4월 8일, 2018년 4월 13일에 취득한 영상으로 차례로 약 24.1도, 40.4도, 52.4도의 관측각을갖는다. 3장의 영상 모두 Ascending Track, Right looking, HH편광모드로 촬영되었다.
Table 1. Summary of experimental data of TerraSAR-X staring spotlight mode
3. 연구방법
본 연구에서는 고해상도 다중관측각 SAR영상간 정합 기법을 비교 분석하고자 한다. 비록 이를 위해Terra SAR-X Staring Spotlight모드 L1A영상을 사용하여 개략정합 단계와 정밀정합 단계로 구성된2단계의 정합기법을 수행한다(Fig.2). 개략정합 단계에서는 기준영상과 대상영상 사이의 상대적으로 큰회전, 이동, 스케일에 대한 차이를 처리하기 위하여 적응형 샘플링 기법과 SAR- SIFT 기법을 사용하여 정합을 수행한다. SAR-SIFT는 SIFT 알고리즘을 SAR영상의 특징인 스펙클노이즈의 영향을줄이도록 개선한 알고리즘이다. 본 연구에서는 SAR-SIFT의 이미지 사이즈가 커질수록 매칭점이 늘어나고 그에 따른 시간과 오탐지율이 높아지는 SAR-SIFT 기법을 보완하기 위하여 매회마다 절반의이미지 사이즈로 평균다운 샘플링을 수행한다. 샘플링 과정을 통해 처리 시간이 감소하고, 기준영상과 대상영상 사이의 기하학적 차이가 줄어들면서 정합성능이 늘어나게 된다.
정밀정합 단계에서는 총 4가지 방법을 통해 비교분석을 수행한다. 변환모델을 통한 NCC와 PC 맵을 활용한 PC-NCC, MI기법을수행하며, 모든 픽셀에서 오프셋을 추정하는 Gefolki 기법을 수행한다. 변환모델을 사용하는 기법들은 일정패치 간격으로 윈도우 기반 지역적 이동량을 추정한다. 해당 이동량을 포함한 변환 최종변환 모델을 통해 정합을 수행하게 된다. Gefolki 기법은 모든 픽셀에 대한 오프셋이 계산되기 때문에, 변환모델 없이 정합을 수행한다. 전체적인 알고리즘흐름도는 Fig. 2와같다.
Fig. 2. Data processing flow chart for coarse-to-fine registration.
1) Coarse registration
SAR-SIFT는 기존 SIFT 기법을 SAR 영상에 적합하도록 변형한 알고리즘으로 SIFT 기법의 차분기울기 연산자를 영역비 기울기 연산자로 대체하여 스펙클노이즈의 영향을 줄였다. SAR-SIFT는 총 4단계의 구성으로, 키포인트 탐지(Key-point detection), 방향성분 지정(Orientation Assignment), 서술자 추출(Descript or extraction), 키포인트 매칭(Key-point matching)으로나눠진다. 키포인트탐지단계에서 logscale의 ROWEA(Ratio of Exponential Weighted Average) 연산자를 통한 기울기성분을 추출한다. 수평성분(식1), 수직성분(식2)을 추출하며, 그 식은 아래와 같다.
\(\begin{gathered} R_{\text {horisontal, } \alpha}= \\ \log \left(\frac{\sum_{i=-M / 2}^{M / 2} \sum_{j=1}^{N / 2} I(x+i, y+j) e^{-\frac{|i|+|j|}{\alpha}}}{\sum_{i=-M / 2}^{M / 2} \sum_{j=-N / 2}^{-1} I(x+i, y+j) e^{-\frac{|i|+|j|}{\alpha}}}\right. \end{gathered}\) (1)
\(\begin{gathered} R_{\text {vertical, } \alpha}= \\ \log \left(\frac{\sum_{i=1}^{M / 2} \sum_{j=-N / 2}^{N / 2} I(x+i, y+j) e^{-\frac{|i|+|j|}{\alpha}}}{\sum_{i=-M / 2}^{-1} \sum_{j=-N / 2}^{N / 2} I(x+i, y+j) e^{-\frac{|i|+|j|}{\alpha}}}\right. \end{gathered}\) (2)
해당 식에서 M과 N은 ROWEA 필터의 사이즈, x, y는 픽셀 위치, α는 스케일 가중 파라미터, I는 픽셀 Intensity 값을 의미한다. 이후 각 스케일별로 SAR Harris map이 생성되며 각 스케일별로 극점을 찾아 키포인트 후보군으로 추출한다. SAR Harris map 생성식은 아래와 같다 (식3).
\(\begin{aligned} &C_{S H}=\operatorname{Gaussian}_{2 \alpha}^{-} *\left[\begin{array}{cc} (R(h, \alpha))^{2} & R h, \alpha \cdot R_{v, \alpha} \\ R h, a \cdot R v, a & (R(h, \infty))^{2} \end{array}\right]\\ &R_{s H}(x, y, \alpha)=\operatorname{det}\left(C_{s H}(x, y, \alpha)\right)-\mathrm{d} \cdot \operatorname{tr}(C s H(x, y, \infty)) \end{aligned}\) (3)
여기서 CSH는 SAR-Harris행렬을 의미하며 Rh, α는 식 1에 의해 생성된 수평성분· Rv, α는 식 2에의해 생성된 수직 성분을 의미한다. RSH(x, y, α)는 SAR-Harris 행렬에 대한 응답함수로 d는 임의값을 의미한다.
방향 성분지정과 서술자 추출단계에서는 스케일별로 기울기 성분의 방향(Orientation) 히스토그램 중 가장 피크 값을 선택함으로써, 회전 불변성을 획득한다. 서술자 추출단계에서는 원형의 로그폴라좌표변환을 통해 서술자를 구성한다. 마지막으로 키포인트 매칭단계에서는NNDR(Nearest Neighbor Distance Ratio)기법을 통해 매칭한다(Delinger et al., 2015).
개략 정합 단계에서는 SAR-SIFT를 적용하기에 앞서, 적응형 샘플링 기법을 적용한다. 개략 정합 단계의 처리시간을 감소시키고 지형적 효과를 완화하기 위하여 영상의 사이즈를 1000픽셀 미만이 되도록 다운 샘플링하여 SAR-SIFT 기법을 적용한다. 적용 후 매칭점들에 RANSAC(RANdom SAmple Consensus) 기법을 적용하여 이상점을 제거하여 정합 오탐지율을 줄인다(Fischler and Bolles, 1981). 개략 정합 단계의 최종 매칭점의 개수가 10개 미만이 될경우 다운샘플링 간격을 1칸 늘려 반복 적용하고 어파인 변환모델을 적용하여 개략정합을 수행한다.
개략정합에서 다운샘플링 통해서 대상영상과 기준 영상에서 영상의 지역적인 패턴이 사라지고, 스펙클 노이즈효과 또한 감소되면서 성공적으로 SAR-SIFT 기법을 적용할 수 있다. 기준 영상과 대상 영상 사이의 큰 이동, 회전, 스케일에 대한 차이들이 제거되면서 정밀 정합 단계의 적용이 가능하다.
2) Fine registration
정밀 정합단계에서 정합 기법별 성능 비교 분석을 위한 강성 정합 모델로는 NCC, PC-NCC, MI3 가지 기법, Non-rigid모델은 Gefolki기법을 적용한다. 강성정합 모델은 윈도우 기반의 패치별 MI나NCC의 유사도값이 최대가 되는 지점을 오프셋으로 추정한다. 영상을 격자 단위로 나눠 해당 패치의 중앙의 윈도우와 주변 픽셀윈도우를 비교하여 오프셋을 추정하며 이때 해당 패치의 픽셀값에서 일정임계값 이상 벗어나는 경우 이상값으로 간주하여 변환모델 생성시에 제거한다(Fig. 3). 개략정합 단계에서 생성된 영상에 다시 변환모델을 적용 할 경우 변환과정에서 생기는 영상왜곡이 중첩되기 때문에 정밀정합 단계의 오프셋 추정이 완료되면 개략정합 단계에서 생성된 변환모델을 통해 추정된정밀정합 단계의 매칭점을 역변환하여 본래 기준영상과 대상영상 좌표에서의 매칭점을 생성한다. Gefolki모델은 조밀비 강성 정합모델로 모든 픽셀에서 오프셋을 계산하기 때문에, 개략정합에서 생성된 영상에서 알고리즘을 적용한다.
Fig. 3. Concept of template matching: Template matching to estimate the offset (upper) at every patch (lower).
(1) NCC
NCC기법을 이용한 템플릿 윈도우 매칭기법은 영상정합에서 기본적인 통계적 접근방법으로 템플릿을통해 오프셋을 추정한다(YooandHan, 2009). NCC 식은 아래와 같다(식4).
\(\gamma(u, v)=\frac{\sum_{x, y}\left[f(x, y)-\bar{f}_{u, v}\right][t(x-u, y-v)-\bar{t}]}{\left.\sum_{x, y}\left[f(x, y)-\bar{f}_{u, v}\right]^{2} \sum_{x, y}[t(x-u, y-v)-t]^{-2}\right\}^{0.5}}\) (4)
여기서 f는 영상의 intensity값, t는 윈도우에서의 intensity 값을 의미한다. \(\bar{f}_{u, v}\)는 윈도우 영역에서의 영상평균값, \(\bar t\)는 윈도우값의 평균을 의미한다. 처리시간을 감안하여 NCC기법의 윈도우사이즈는 64×64, 패치 사이즈는 32×32로 수행했다.
(2) PC-NCC
Intensity영상에 PC기법을 사용하면, 엣지가 강조된 영상생성이 가능하다. PC기법은 기존의gradient를 이용한 엣지강조 기법들은 영상밝기값 변화에 민감하고 올바르게 엣지를 추출하지 못하는 문제점을 주파수 도메인에서의 위상값을 통해 해결할 수 있다. PC-NCC 기법은 기준영상과 대상영상에 PC기법을 사용하여 엣지강조 영상을 생성한 후(Fig. 4), NCC 기법을 적용하는 기법이다. PC기법은 주파수 도메인에서 대푯값을 찾는 접근법으로 위상의 합동정도를 나타내는 0에서1 사이의값을 갖으며 식 5와 같이 표현된다(Kovesi, 2000).
\(P C(x)=\frac{\sum_{n, o} W_{o}(x)\left|A_{n, o}(x) \Delta \emptyset_{n, o}-T\right|}{\sum_{n, o} \mathrm{~A}_{n, o}(x)+\varepsilon}\) (5)
여기서 o는 orientation, n은 스케일을 의미한다. W(x) 는주파수성분에대한가중함수, Δøn, o는위상각(Phase angle), An, o은 주파수 영역에서의 진폭값, ε는 An,o가 작은값이 되었을때 식이 불안정해지는 것을 대비한 상수, T는 threshold를 의미한다. 엣지 성분을 두드러지게 강조하는 Maximummoment PCmap영상(Fig. 4)에 NCC 기법을 적용하여 정밀정합 단계를 수행한다. NCC와 마찬가지로 기법의 윈도우 사이즈는 64×64, 패치 사이즈는 32×32로 수행했다
Fig. 4. Example of amplitude of TSX image and PC (Phase Congruency) map.
(3) MI
MI는 식6과 같이 두 개의 무작위변수의 관계를 엔트로피 형식으로 나타낸다(Studoholme et al., 1999). 기존의 NCC기법이 두 변수간 선형관계여야만 한다는 제약조건에 비해 두 변수가 비선형관계임에도 활용이 가능하다는 장점이있다.
\(M I(A, B)=\frac{H(A)+H(B)}{H(A, B)}\) (6)
H(A)와 H(B)는 Shannon 엔트로피값을 나타내며 H(A, B)는 두 엔트로피의 교집합을 의미한다. H(A)와 H(B)는 한계질량함수(marginal probability mass function), H(A, B)는 결합확률 질량함수(joint probability mass function) 에 의해 계산된다(Chen et al., 2003). MI는 윈도우 사이즈를 128×128로 늘려 신뢰성을 높이는 대신 패치사이즈를 8×8로 줄여 처리시간을 줄여 수행했다.
(4) Gefolki
Gefolki는 optical flow 알고리즘의 일종으로 컴퓨터 비전분야에서 객체의 움직임을 추적하기 위하여 사용되어 진다. 주로 첫번째 프레임의 밝기 값과 다음 프레임에서의 픽셀 밝기값은 식 6과 같이 표현된다.
\(I_{1}=I_{2}(x)+\nabla I_{2} \cdot u\) (6)
여기서 I1는 첫번째 프레임의 영상, I2는 두번째 프레임의 영상을 의미하며 u는 두 영상 사이의 변위를 의미한다.
식 6은 두 영상 프레임간 밝기값 변화량이 아주 작다는 가정과 관측객체의 주변 픽셀분포가 유사하다는 2개의 가정에 의해서 성립하는데 원격탐사 분야에서는 이러한 가정이 모두 위배된다(Lucasand Kanade, 1981). 두 영상간의 시간간격이 크고 관측기하에 따라 영상 밝기값이 변화하기때문인데, 이는 스무딩을 통하여 해당 식을 일반화할 수 있다(Brigot et al., 2016). Gefolki에서는 local window-based 접근법을 기반으로한 Folki 알고리즘을 활용하였으며, 광학흐름(Opticalflow) 식을 1차 테일러전개를 기반으로한 Gauss-Newton기법을 반복하는 것으로 계산했다. Gefolki는 목적함수 인식 7을 최소화 하도록 구성되어있다(Champagnat et al., 2011). Imagepyramid기법을 적용하여 기준영상과 대상영상 사이의 큰 변위부터 미세한 변위를 찾아가는 방식을 추가로 결합하였으며, Gauss-Newton 반복법을 이용하여 J(u; x)를 최소화한다.
\(J(u ; \mathrm{x})=\sum_{x^{\prime} \in S} r\left(x^{\prime}-x\right)\left(R\left(I_{1}\left(x^{\prime}\right)\right)-R\left(I_{2}\left(x^{\prime}+u(x)\right)\right)^{2}\right.\) (7)
Fig. 5. Processing flow chart of Gefolki registration.
이때 S는 이미지 스케일 단계를 의미하고, r은 local window smoothing 필터, R(x) 는 rank transform 함수를 의미한다. Coarse 단계부터 fine 단계로의 피라미드를 구성하고 각 피라미드별J(u;x)값을 최소화하는 식으로 구성된다(Fig. 5).
Rank transform은 각 픽셀값보다 픽셀값이 작은 이웃한 픽셀의 개수로 변환하는 기법으로 식 8과 같이 표현된다.
\(\begin{aligned} R(I)(x)=& \#\left\{\mathrm{x}^{\prime}: \mathrm{x}^{\prime} \in \text { Neighborhood }(x)\right. \text { with }\\ &\left.|I(x)|>\left|I\left(x^{\prime}\right)\right|\right\} \end{aligned}\) (8)
Rank transform은 비선형필터(nonlinear filter)로써 픽셀의 상대적 밝기값에 따라 영상밝기값을압축한다 (Zabihand Woodfill, 1994). Rank transform을 통해 영상을 압축함으로써 더욱 강건한기울기 성분의 추출이 가능하다(Fig. 6).
Fig. 6. Example of (a) amplitude image and (b) rank transform image.
4. 연구 결과
2018년 4월3일, 2018년 4월 8일 2018년 4월 13일에 획득된 총 3개의 관측각 차이가 나는 영상을통하여 정합을 수행했다(Table1). 영상 쌍은 총 3개로 구성하였으며 구성은 Table2와 같다. 관측 각 차이가 가장 많이 나는 쌍은 2018년 4월 3일 영상을 기준영상으로 하고, 2018 년4월13일을대상영상으로하는 1번째 영상 쌍이며, 28.316°의 관측각 차이가 존재한다. 가장 관측각 차이가 적게 나는쌍은 2번째 쌍으로 약 11.963°의 관측각 차이가 있다. 3번째 쌍은 2018년 4월 8일과 4월 3일을 기준영상, 대상영상으로 하는 쌍으로 약 16.353°의 관측각 차이가 있다.
Table 2. List of registration pairs with the difference of incidence angle
1) Coarse registration 단계 정합 결과
개략정합 단계에서는 적응형 샘플링기법과 SAR- SIFT를 이용한 개략정합을 수행하였으며, 정합결과는 Fig. 7과같다. 적응형 샘플링 방식을 적용한 결과 1번 영상 쌍과 2번 영상 쌍에는 다운샘플링 간격을 21칸으로 적용하였으며 3번 영상 쌍은 19칸 간격의 다운 샘플링이 적용되었다. 대상 영상은 기준 영상의 사이즈에 맞춰 리샘플링 되었다. 1, 2, 3영상 쌍은 각 41개, 221개, 89 개의 매칭점을 찾았다. 다운 샘플링을 적용했음에도 불구하고 관측각 차이에 따른 영상 특성이 상이하게 달라지기 때문에 매칭점의 개수 또한 그에 비례하여 나타났다. 모든 영상 쌍의 매칭점의 개수가 10개 이상 추출되었기 때문에 추가 다운 샘플링은 수행하지 않았다.
Fig. 7. Coarse registration results on three pairs in Table 2.
2) Fine registration 단계의 파라미터 세팅
Fine 단계에서의 정합 기법 별 파라미터는 Table3과 같다.
Table 3. Parameter settings and optimized values for each method
NCC기법은 충분한 수의 정합점을 취득하기 위하여 32×32의 패치사이즈를 사용하였고 윈도우사이즈는 64×64로 지정했다. 그러나 MI기법의 경우 처리시간이 지나치게 오래 걸리기 때문에 패치사이즈를 8×8로 줄이는 대신 윈도우 사이즈를 256×256로 늘려 패치별 정합점의 신뢰성을 확보했다.
Gefolki 기법의 Rank filter size와 Gauss Newton iteration 횟수 파라미터는 경험적으로 결과에 큰 영향을주지 않았기 때문에, Plyer et al. (2015)의 기본 파라미터 세팅을 따랐다. Radius 사이즈와 pyramid level은 개략 정합단계의 정확도와 영상 사이즈 크기에 따라 다르게 사용해야 한다. Radius는 큰 값을 가질수록 스펙클 노이즈의 영향을 줄일 수 있고, 두 영상간 상대적으로 큰 오프셋을 추정할 수 있다. Pyramid level 또한 근본적으로 radius 사이즈와 연관되어져 있는데, pyramid level 값을 크게 할 수록 큰 오프셋을 추정할 수 있지만, pyramid level값을 지나치게 높게 주면 영상의 고유한 패턴이 사라질 위험이 있기 때문에, 이를 고려하여 사이즈를 정해야 한다. 본 연구에서 사용한 영상의 사이즈가 충분히 크고, Coarse 단계에서 21번의 다운 샘플링 후 개략 정합을 수행했기 때문에, 이를 고려하여 radius 사이즈와 pyramid level을 파라미터를 세팅했다. Gefolki 적용결과로 구한 u(range), v(azimuth) 성분은 Fig. 8와 같다. Range 방향의 오프셋 성분 u는 연구 지역의 DEM(Fig. 8)과 매우 유사한 형태로 나타난 것으로 보아, 관측 각도 차이에 의한 지형 기복이 잘 반영된 것으로 판단된다.
Fig. 8. Result of Gefolki offset (a) u-component and (b) v-component, and (c) Digital elevation model of study area.
3) 정합 성능 분석
정합기법별 성능분석을 위하여 모자이크 이미지를 통한 정성적 분석, FSIM(Feature similarity) 지수와 RMSE (Root Mean Square Error)를 통한 정량적 분석을 수행 했다. RMSE는 영상 내의 두드러지는 10개 지점을 선정하여 계산하였으며 식 9와 같이 표현된다.FSIM 지수는 영상 품질 평가 방법의 일종으로 대표적인 방법인 SSIM (Structural similarity) 방법을 개선하여 PC(Phase Congruency) 맵과 GM(Gradient magnitude) 맵 두 가지를 사용하여 최종유사도를 측정하는 방식이다(Zhang et al., 2011).
\(R M S E=\sqrt{\sum_{i=1}^{n} \frac{\left(\bar{x}_{i}-x_{i}\right)^{2}+\left(\bar{y}_{i}-y\right)^{2}}{n}}\) (9)
Table 4. Performance comparison between fine registration methods (NCC, PC-NCC, MI, Gefolki)
여기서 n은 총 검사점의 수이고 \(\bar x\), \(\bar y\)는 정합 완료된 대상 영상의 좌표, x, y는 기준 영상에서의 좌표를 의미한다.
영상정합결과, 전반적으로 Gefolki 알고리즘의 FSIM 지수와 RMSE 정확도가 상대적으로 우수하게 나타났다. 3개의 영상 쌍에서 RMSE와 FSIM 지수는 각 각 3.438, 1.651, 3.246 화소, FSIM 지수는 0.7775, 0.86, 0.8115로 가장 우수한 정합 성능을 보였다. 소모시간은NCC기법이 27.98초, 22.58초, 28.39초로 가장 짧은 소모시간을 보여주었으며 PC-NCC의 경우 PCmap을 생성하는 시간을 제외하고 동일한 처리시간을 보여주었다. Gefolki기법은 GPU(Graphics Processing Units)를 사용하여 수행 되었으며, NCC 기법보다 오랜 시간이 걸렸지만 상대적으로 빠른 처리속도를 보여주었다. MI기법은 다른 기법에 비해 소모시간이 길고 RMSE 정확도와 FSIM 지수 또한 다른 강성 정합기법에 비해 저조한 결과를 보였다. NCC 기법과 PC-NCC 기법의차이는 단지 NCC 적용 대상이 PC 맵과 intensity 맵의 차이이기 때문에, 엣지 구조물들이 있는 지역에서는 PC-NCC의 성능이 더 좋을 것으로 사료된다. 전체적인 RMSE와 FSIM 인덱스의 추세는 관측각 차이가 가장 적은 2018년 4월 8일과 4월 13일을 대상으로한 2번 영상 쌍의 정합에서 가장 우수하게 나타났으며, 관측각 차이가 가장 큰 2018년 4월 3일과 4월 13일을 대상으로한 1번 영상 쌍의 정합에서 가장 저조한 결과를 보였다(Table4).
다중 관측 각 영상정합에서 고도가 있는 지역에서는 기하학적 왜곡이 발생하기 때문에 취약하다. 본 논문의 연구 영상 내에서 지형 고도가 위치하는 곳은 우측 산악 지형과 가운데 인공구조물이 존재하는 지역이다(Fig. 8). 정성적 평가를 위하여 고도가 있는 부분을 확대하여 정합 성능을 비교했다.
1 번 영상 쌍은 28.316°의 관측각 차이가 나는 영상 쌍으로, rigid 기법을 적용한 영상의 확대 이미지를 보면 산악지형 주변에서 크게 이격이 존재하는 것에 비해 고도가 없는 지역에서는 상대적으로 정합이 잘 수행되었음을 확인할 수 있다(Fig. 9). 산악지형에서 나타나는 강성 정합 알고리즘 확대 영상에서 어긋나게 나타난 부분들이 Gefolki 알고리즘 적용 영상에서는 상대적으로 정렬된 모습을 확인할 수 있다. 여전히 1-4 pixel 차이가 존재하지만, 고도 지역에서 수십 픽셀이 상이 격이 존재하는 강성 정합기법에 비해 우수한 결과를 보였다.
Fig. 9. First pair registration result and enlarged images. (a) reference image; (b) coarse registered; (c) mosaic (a) and (b) with Gefolki; (d)-(g) the enlarged images of NCC, PC-NCC, MI, Gefolki, (h)-(k) the enlarged images of NCC, PC-NCC, MI, Gefolki.
2 번 영상 쌍은 11.963° 관측각 차이가 나는 영상 쌍으로 정합 결과는 Fig. 10과 같다. 가장 작은 관측각 차이가 존재하는 영상 쌍으로 가장 좋은 정합 결과를 보였으나, 강성 정합 모델은 산악 지형 부근에서 정합 오차가 발생한다.
Fig. 10. Second pair registration result and enlarged images. (a) reference image; (b) coarse registered; (c) mosaic (a) and (b) with Gefolki; (d)-(g) the enlarged images of NCC, PC-NCC, MI, Gefolki.
관측각 차이가 16.353° 나는 Imagepair3의 정합 결과는 Fig. 11과 같다. 마찬가지로 강성정합모델은 산악지형에서 큰 오차를 보였으며, Gefolki 알고리즘이 적정 수준의 정합 결과를 보였다.
Fig. 11. Third pair registration result and enlarged images. (a) reference image; (b) coarse registered; (c) mosaic (a) and (b) with Gefolki; (d)-(g) the enlarged images of NCC, PC-NCC, MI, Gefolki.
5. 결론
본 논문에서는 고해상도 SAR 영상에 대해 다중 관측각 환경에서 정밀 정합 기법 비교 분석을 수행하였다. 실험자료는 촬영각도가 상이한 3개의 TerraSAR-Xstaring modeL1A 영상으로부터, 생성된 3개의 영상 쌍을 사용하였다. 영상정합은 개략 정합 단계와 정밀 정합 단계를 통한 2단계의 정합 알고리즘으로 구성되었다. 개략 정합 단계에서는 적응형 샘플링 기법과 SAR-SIFT 알고리즘을 결합하여 기하학적 왜곡과 스펙클 노이즈에 대한 영향을 줄임으로써 정합을 수행했다. 개략정합 결과 영상 쌍의 관측각 차이에 따라 매칭점의 개수가 반비례하여 나타났으며, 이는 관측각 차이가 커짐에 따라 나타나는 것으로 사료된다. 그럼에도 본 연구에서 사용한 모든 영상 쌍에서, 최소 매칭 점 개수인 10개 이상을 매칭하여 성공적으로 개략정합을 수행하였다.
정밀정합 단계에서는 강성모델인 NCC, PC-NCC, MI와 비강성 모델인 Gefolki를 수행 하였으며, RMSE 와 FSIM 지수를 통해 정량적으로 결과를 비교 분석 하였다. 강성모델들은 4~10 RMSE 화소 결과로 저조한 성적을 보였으나. 비강성모델인 Gefolki 알고리즘은 1~3 RMSE 화소로 우수한 정확도를 보여주었다. NCC 알고리즘은 약 30초 내외의 우수한 자료처리 시간을 보였으며 Gefolki 알고리즘의 자료처리 시간 또한 50초 내외로 준수한 결과를 보였다. 전반적으로 모든 정합기법에서 영상쌍 관측각 차이가 작을수록 우수한 결과를 보여주었으며, 관측각 차이가 가장 큰 2018년 4월 3일과 4월 13일 영상 쌍에서는 가장 저조한 결과를 보였다. 정성적 분석결과 강성모델은 지형 기복이 큰 지역에서 오차가 발생함을 확인할 수 있었으며, 비강성 모델인 Gefolki 알고리즘 또한 관측각 차이에 따른 고도지역에서의 영향을 받긴 하지만, 다른 기법대비 가장 우수한 결과를 보였다.
본 논문에서는 CoarsetoFine 단계를 활용하여 다중 관측각 영상에서 신속하게 자동 정밀 정합을수행하였다. Gefolki 알고리즘은 다중 관측각 SAR 영상에서 정밀 DEM, DSM 없이 1~3 화소 이내의 유의미한 정확도를 취득했다는 점에서 군사, 보안 분야 등 다양한 적용 가능 분야가 있을 것으로 사료된다. 그러나 여전히 관측각 차이가 커질수록 정합성능이 떨어지는 문제가 있고 처리시간 또한 알고리즘 최적화를 통해 개선해야 할 필요성이 있다. 향후 본 연구를 기반으로 다른 관측트랙, 다른 관측 방향과 같은 복잡한 환경에서의 정합 기법 개발 또한 가능할 것으로 기대된다.
사사
이 논문은 국방 과학 기술 연구소의 지원(UD190024 FD)을 받아 수행된 연구이며, 이에 감사드립니다.
References
- Brigot, G., E. Colin-Koeniguer, A. Plyer, and J. Fabrice, 2016. Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9(7): 2923-2939. https://doi.org/10.1109/JSTARS.2016.2578362
- Champagnat, F., A. Plyer, G. Le Besnerais, B. Leclaire, S. Davoust, and Y. Le Sant, 2011. Fast and accurate PIV computation using highly parallel iterative correlation maximization, Experiments in Fluids, 50(4): 1169-1182. https://doi.org/10.1007/s00348-011-1054-x
- Chen, H.M., P.K. Varshney, and M.K. Arora, 2003. Performance of mutual information similarity measure for registration of multitemporal remote sensing images, IEEE Transactions on Geoscience and Remote Sensing, 41(11): 2445-2454. https://doi.org/10.1109/TGRS.2003.817664
- Chen, H.M., M.K. Arora, and P.K. Varshney 2003. Mutual information-based image registration for remote sensing data, International Journal of Remote Sensing, 24(18): 3701-3706. https://doi.org/10.1080/0143116031000117047
- Fischler, M.A. and R.C. Bolles, 1981. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography, Communications of the ACM, 24(6): 381-395. https://doi.org/10.1145/358669.358692
- Goncalves, H., L. Corte-Real, and J.A. Goncalves, 2011. Automatic image registration through image segmentation and SIFT, IEEE Transactions on Geoscience and Remote Sensing, 49(7): 2589-2600. https://doi.org/10.1109/TGRS.2011.2109389
- Jeon, H. and Y. Kim, 2018. Registration Method between High Resolution Optical and SAR Images, Korean Journal of Remote Sensing, 34(5): 739-747 (in Korean with English abstract) https://doi.org/10.7780/kjrs.2018.34.5.3
- Kovesi, P., 2000. Phase congruency: A low-level image invariant, Psychological Research, 64(2): 136-148. https://doi.org/10.1007/s004260000024
- Lowe, D.G., 2004. Distinctive image features from scale-invariant keypoints, International Journal of Computer Vision, 60(2): 91-110. https://doi.org/10.1023/B:VISI.0000029664.99615.94
- Lucas, B.D. and T. Kanade, 1981. An iterative image registration technique with an application to stereo vision, Proc. of International Joint Conference On Artificial Intelligence, Vancouver, BC, CA, Aug. 24-28, vol. 81, pp. 674-694.
- Plyer, A., E. Colin-Koeniguer, and F. Weissgerber, 2015. A new coregistration algorithm for recent applications on urban SAR images, IEEE Geoscience and Remote Sensing Letters, 12(11): 2198-2202. https://doi.org/10.1109/LGRS.2015.2455071
- Sansosti, E., P. Berardino, M. Manunta, F. Serafino, and G. Fornaro, 2006. Geometrical SAR image registration, IEEE Transactions on Geoscience and Remote Sensing, 44(10): 2861-2870. https://doi.org/10.1109/TGRS.2006.875787
- Schwind, P., S. Suri, P. Reinartz, and A. Siebert, 2010. Applicability of the SIFT operator to geometric SAR image registration, International Journal of Remote Sensing, 31(8): 1959-1980. https://doi.org/10.1080/01431160902927622
- Song, Z.L. and J. Zhang, 2010. Remote sensing image registration based on retrofitted SURF algorithm and trajectories generated from Lissajous figures, IEEE Geoscience and Remote Sensing Letters, 7(3): 491-495. https://doi.org/10.1109/LGRS.2009.2039917
- Studholme, C., D.L. Hill, and D.J. Hawkes, 1999. An overlap invariant entropy measure of 3D medical image alignment, Pattern Recognition, 32(1): 71-86. https://doi.org/10.1016/S0031-3203(98)00091-0
- Suri, S. and P. Reinartz, 2010. Mutual-information-based registration of TerraSAR-X and Ikonos imagery in urban areas, IEEE Transactions on Geoscience and Remote Sensing, 48(2): 939-949. https://doi.org/10.1109/TGRS.2009.2034842
- Thiele, A., E. Cadario, K. Schulz, U. Thonnessen, and U. Soergel, 2007. Building recognition from multi-aspect high-resolution InSAR data in urban areas, IEEE Transactions on Geoscience and Remote Sensing, 45(11): 3583-3593. https://doi.org/10.1109/TGRS.2007.898440
- Wang, F., H. You, and X. Fu, 2014. Adapted anisotropic Gaussian SIFT matching strategy for SAR registration, IEEE Geoscience and Remote Sensing Letters, 12(1): 160-164. https://doi.org/10.1109/lgrs.2014.2330593
- Wang, S, H. You, and K. Fu, 2012. BFSIFT: A novel method to find feature matches for SAR image registration, IEEE Geoscience and Remote Sensing Letters, 9(4): 649-653. https://doi.org/10.1109/LGRS.2011.2177437
- Yoo, J.C. and T.H. Han, 2009. Fast normalized cross-correlation, Circuits, Systems and Signal Processing, 28(6): 819-843. https://doi.org/10.1007/s00034-009-9130-7
- Zabih, R. and J. Woodfill, 1994. Non-parametric local transforms for computing visual correspondence, Proc. of Third European conference on computer vision, Stockholm, SE, May. 2-6, vol. 801, pp. 151-158.
- Zhang, L., L. Zhang, X. Mou, and D. Zhang, 2011. FSIM: A feature similarity index for image quality assessment, IEEE transactions on Image Processing, 20(8): 2378-2386. https://doi.org/10.1109/TIP.2011.2109730