基于DEM的SAR影像直接與間接正射校正方法比較

李艷艷，唐 娉，胡昌苗，單小軍

1.中國科學(xué)院 遙感與數(shù)字地球研究所，北京 100101

2.中國科學(xué)院大學(xué)，北京 100049

1 引言

合成孔徑雷達(dá)（Synthetic Aperture Radar，SAR）具有全天時(shí)、全天候的觀測特點(diǎn)，日益受到人們的關(guān)注。然而，由于其側(cè)視成像的成像機(jī)理，SAR影像的質(zhì)量受地形影響嚴(yán)重，影響了人們的理解與應(yīng)用，因此，需要對SAR影像進(jìn)行正射校正處理來消除幾何變形[1]。SAR影像正射校正是消除由于側(cè)視成像和復(fù)雜地形引起的變形，使其具有真實(shí)地理位置和地表信息的過程[2]。

基于數(shù)字高程模型（Digital Elevation Model，DEM）數(shù)據(jù)構(gòu)建距離-多普勒（Range Doppler，RD）模型對SAR影像進(jìn)行正射校正[3]是近年來發(fā)展較快的一類方法，根據(jù)出發(fā)點(diǎn)的不同，可以分為直接校正方法和間接校正方法。直接校正[4-7]是從SAR成像幾何空間坐標(biāo)出發(fā)，通過RD模型直接解算出每一個(gè)像元所對應(yīng)的大地參考空間坐標(biāo)，完成校正工作；間接校正[8-11]是從大地參考空間坐標(biāo)出發(fā)，通過RD模型解算出每一個(gè)像元對應(yīng)的模擬SAR影像坐標(biāo)，生成模擬SAR影像，再通過模擬SAR影像與真實(shí)SAR影像的配準(zhǔn)，進(jìn)而建立SAR成像幾何空間坐標(biāo)與大地參考空間坐標(biāo)的映射關(guān)系，間接完成SAR影像的校正工作。

已發(fā)表的關(guān)于SAR影像正射糾正的文獻(xiàn)，大多集中于自身模型原理與應(yīng)用方法的闡述，對于兩種方法校正效果與適用性的比較分析，鮮有文獻(xiàn)進(jìn)行論述。因此，本文將利用同一地區(qū)DEM數(shù)據(jù)對兩種方法的校正步驟、校正效果與精度進(jìn)行分析比較。

2 RD模型簡介

RD定位模型是根據(jù)SAR影像的成像機(jī)制，對SAR影像進(jìn)行定位的方法。對于斜距影像上任意一點(diǎn)，滿足以下條件：

其中，N為方位向視數(shù)，PRF為SAR的脈沖重復(fù)頻率，R0是SAR近地點(diǎn)斜距，Mslant是SAR影像的斜距分辨率，均可在參數(shù)文件中讀??；t為當(dāng)前時(shí)間；R為衛(wèi)星和地面點(diǎn)之間的斜距。

實(shí)驗(yàn)中的衛(wèi)星軌道記錄在地心旋轉(zhuǎn)坐標(biāo)系（ECR）中，地面點(diǎn)相對于坐標(biāo)軸是靜止的。運(yùn)動(dòng)速度為零，因此建立RD模型如下：

其中，(XS,YS,ZS,VX,VY,VZ)為衛(wèi)星的三維位置和三維速度，(X ,Y,Z)為地面點(diǎn)P的三維位置，fD為地面點(diǎn)P的多普勒頻移，λ為雷達(dá)波長。

式（2）的第一個(gè)公式是距離公式，幾何意義是地面點(diǎn)P坐標(biāo)位置分布在距離衛(wèi)星為R的圓上。第二個(gè)公式是多普勒頻移公式，如果 fD已知，其幾何意義表示為一簇雙曲線（多普勒曲線）。由于地球自轉(zhuǎn)的原因，這些曲線相對星下點(diǎn)呈不對稱分布[12]。

3 直接校正與間接校正

3.1 直接校正

對于直接校正方法，是從SAR成像幾何空間坐標(biāo)(x ,y )出發(fā)，根據(jù)公式（1）、（2），直接計(jì)算出每個(gè)像元對應(yīng)的地面點(diǎn)坐標(biāo)(X ,Y,Z )，完成正射校正工作，主要經(jīng)過以下幾個(gè)環(huán)節(jié)[13]：

（1）從SAR影像讀取行列坐標(biāo)(x ,y)，根據(jù)SAR的頭文件，進(jìn)行衛(wèi)星軌道的擬合，解算出衛(wèi)星位置隨時(shí)間的變化關(guān)系。

（2）根據(jù)公式（2）的兩個(gè)公式構(gòu)建模型函數(shù)記為F1(x ,y,X,Y,Z)和F2(x ,y,X,Y,Z )，建立線性化的RD模型，C?ΔG-l=0，其中l(wèi)=[F10,F20]T，ΔX,ΔY 分別為X、Y的坐標(biāo)改正量，F(xiàn)10,F20分別為F1,F2的初值，建立誤差方程式，V=C?ΔG-l，記V=[VF1,VF2]T，對誤差方程式建立并求解法方程式，可求得地面點(diǎn)平面坐標(biāo)改正量 ΔG=(CTC)-1CTl，最終得地面點(diǎn)平面坐標(biāo)為：X=X0+ΔX，Y=Y0+ΔY，其中(X0,Y0)為地面點(diǎn)平面坐標(biāo)初值。在已知DEM中，由平面坐標(biāo)(X ,Y )通過雙線性內(nèi)插求得地面點(diǎn)高程Z，給定閾值T，若| Z-Z0|>T，則將內(nèi)插的高程作為高程初值，重新迭代計(jì)算平面坐標(biāo)，否則終止迭代，求得地面點(diǎn)空間直角坐標(biāo)( )X,Y,Z ，并通過坐標(biāo)系轉(zhuǎn)換為大地坐標(biāo)，其中，B為大地緯度，L為大地經(jīng)度，H為大地高，使用最近鄰法確定影像灰度值。

（3）遍歷SAR影像，生成處理前正射校正影像。

（4）以上生成的處理前的校正影像中存在大量的未被賦值的像元，即“空洞”像元，平坦區(qū)域零星點(diǎn)綴，地形起伏較大的區(qū)域呈現(xiàn)大面積的空洞像元，之所以產(chǎn)生空洞像元是側(cè)視成像迎坡像元疊加壓縮造成的，地形起伏越大的地方，疊加壓縮的效應(yīng)越顯著，空洞區(qū)域越大。首先對處理前校正結(jié)果影像進(jìn)行3×3鄰域的內(nèi)插處理，零星空洞的像元處理后可以得到明顯改善，但是對于大區(qū)域的空洞像元是遠(yuǎn)遠(yuǎn)不夠的，還需要進(jìn)一步處理；處理的原則依然是借附近像元的值進(jìn)行填充，只是內(nèi)插鄰域會(huì)根據(jù)空洞區(qū)域的面積自適應(yīng)變大。具體處理的步驟如下：從SAR影像空間出發(fā)，對于每一個(gè)像點(diǎn)坐標(biāo)，解算出像點(diǎn)坐標(biāo)所對應(yīng)的地面點(diǎn)坐標(biāo)，和以上工作不同的地方在于灰度值的確定，此時(shí)要判定校正影像中B、L所在像元周圍的一片區(qū)域K×K（如10×10像元）灰度值是否為零，同時(shí)判定SAR影像像點(diǎn)坐標(biāo)(x,y)的灰度值是否大于某一限定值T1（如200），若校正影像中B、L所在像元周圍的一片區(qū)域K×K內(nèi)某一像元灰度值為零，同時(shí)SAR影像像點(diǎn)坐標(biāo)(x,y)的灰度值大于限定值T1，則將像點(diǎn)坐標(biāo)(x,y)的灰度值賦給此像元，不斷調(diào)整區(qū)域大小，直至所有空洞像元都能被賦值。這樣即可生成處理后的正射校正影像，并對校正結(jié)果進(jìn)行精度驗(yàn)證與分析。

流程圖如圖1所示。

3.2 間接校正

對于間接校正方法，是從地面點(diǎn)空間直角坐標(biāo)(X ,Y,Z )出發(fā)，根據(jù)公式（1）、（2），計(jì)算每個(gè)像元所對應(yīng)的模擬SAR影像的坐標(biāo)(x ,y)，生成模擬SAR影像，建立地面點(diǎn)與模擬SAR影像的映射關(guān)系；再通過模擬SAR影像與真實(shí)SAR影像的配準(zhǔn)，建立模擬SAR影像與真實(shí)SAR影像的映射關(guān)系，進(jìn)而建立地面點(diǎn)與真實(shí)SAR影像的對應(yīng)關(guān)系，間接完成校正工作，主要經(jīng)過以下幾個(gè)環(huán)節(jié)[14]：

圖 1 直接校正流程

（1）從DEM影像讀取地面點(diǎn)大地坐標(biāo)( )B,L,H，并通過坐標(biāo)系轉(zhuǎn)換為空間直角坐標(biāo)( )X,Y,Z。

（2）建立RD模型，對公式（2）稍作變形，構(gòu)建函數(shù)其中，VS為衛(wèi)星速度矢量，組成元素為(VX,VY,VZ)，rS為衛(wèi)星位置矢量，組成元素為(XS,YS,ZS)，rP為地面點(diǎn)位置矢量，組成元素為(X ,Y,Z )。使用牛頓迭代法求解函數(shù)，從t=0開始迭代，給定閾值T，若g(t)>T ，則 Δt=-g(t)/g′(t) ，t=t+Δt，重新迭代，否則結(jié)束迭代，返回t，再利用t計(jì)算斜距R，根據(jù)公式（1）以及t，R值求得對應(yīng)的模擬影像行列坐標(biāo)(x ,y)，使用雙線性內(nèi)插法確定影像灰度值。

（3）遍歷DEM，生成模擬SAR影像。

（4）模擬SAR影像與真實(shí)SAR影像配準(zhǔn)，建立模擬SAR與真實(shí)SAR的映射關(guān)系其中，is，js為模擬影像坐標(biāo)，ir，jr為真實(shí)影像坐標(biāo)，a,b,c,d,e,f,g,h為擬合得到的系數(shù)。

（5）采用雙線性內(nèi)插法確定每一地面點(diǎn)所對應(yīng)的灰度值，生成最終的正射校正影像，并對校正結(jié)果進(jìn)行精度驗(yàn)證與分析。

流程圖如圖2所示。

4 實(shí)驗(yàn)比較

4.1 實(shí)驗(yàn)數(shù)據(jù)

圖2 間接校正流程

為了對兩種方法的效果與精度進(jìn)行分析，采用日本富士山地區(qū)的ALOS-2 PALSAR-2影像和SRTM 30 m分辨率的DEM（來源于USGS網(wǎng)站，http：//earthexplorer.usgs.gov/）進(jìn)行了實(shí)驗(yàn)。

由于單視復(fù)數(shù)據(jù)（Single-Look Complex，SLC）圖像數(shù)據(jù)中包含大量的斑點(diǎn)噪聲，為了降噪，提高輻射分辨率，同時(shí)兼顧DEM的分辨率，對SLC影像進(jìn)行了方位向5視，距離向2視處理，多視處理后的影像見圖3，影像大小為4 642×3 748，分辨率為14.2 m（方位向）×5.72 m（斜距向），入射角為31.1°，因此對應(yīng)的地距向分辨率為11.1 m。

圖4為該地區(qū)的DEM，為了兼顧SAR影像的分辨率，對DEM進(jìn)行了兩倍過采樣處理，有4 746×4 684個(gè)格網(wǎng)點(diǎn)，格網(wǎng)大小為15 m，高程范圍是0～3 759 m，DEM與SAR影像區(qū)域并不完全重合。

圖3 PALSAR-2影像

圖4 DEM影像

4.2 正射校正結(jié)果對比分析

間接校正過程相對來說比較復(fù)雜，需要先基于DEM進(jìn)行影像模擬，建立DEM格網(wǎng)與模擬SAR影像間的映射關(guān)系，圖5為生成的模擬影像，從影像上可以看出，平坦區(qū)域紋理信息較弱，地形起伏較大的區(qū)域，紋理信息很明顯，和真實(shí)SAR影像極為相似。模擬SAR影像與真實(shí)SAR影像進(jìn)行自動(dòng)特征點(diǎn)的檢測與匹配，采用最小二乘的方法建立模擬SAR影像與真實(shí)SAR影像之間的仿射變換關(guān)系：

其中，is為模擬影像橫坐標(biāo)，js為模擬影像縱坐標(biāo)，ir為真實(shí)影像橫坐標(biāo)，jr為真實(shí)影像縱坐標(biāo)。

根據(jù)兩個(gè)變換關(guān)系就可以建立真實(shí)SAR影像與DEM格網(wǎng)間的變換關(guān)系，從而將SAR影像校正到DEM所在的地理空間。

而直接校正過程相對來說比較直觀，從SAR影像空間出發(fā)，建立RD模型。依據(jù)DEM輔助數(shù)據(jù)，直接計(jì)算出所對應(yīng)的地理空間坐標(biāo)即可，但處理前校正影像如圖6所示，存在大量的未被賦值的像元（稱為空洞像元），平坦區(qū)域零星點(diǎn)綴，圖7所示，地形起伏較大的區(qū)域呈現(xiàn)大面積的空洞像元，圖8所示，這就需要對初步校正影像進(jìn)行進(jìn)一步處理，生成最終校正影像。

圖9給出了間接校正后的校正影像，圖10給出了直接校正處理后影像，圖11為崎嶇區(qū)域局部放大圖，圖12為崎嶇區(qū)域局部放大圖，圖13為平坦區(qū)域局部放大圖，圖14為平坦區(qū)域局部放大圖。

對比圖9和圖6，發(fā)現(xiàn)兩者紋理基本一致，但輻射值相差很多，尤其是地形變化較大的區(qū)域會(huì)更加明顯，如圖11、圖8所示，由于SAR影像側(cè)視成像的成像原理，SAR影像中存在疊掩、陰影、透視收縮等現(xiàn)象，直觀的表現(xiàn)就是迎坡被壓縮，且輻射值是多個(gè)值的疊加值，會(huì)很大。間接校正的過程中，是從DEM空間出發(fā)，在SAR影像中尋找其對應(yīng)點(diǎn)，并將其輻射值賦給校正影像。因此，校正影像中不會(huì)出現(xiàn)“空洞”像元，且迎坡會(huì)被賦予SAR影像中疊加后的輻射值，普遍偏亮，高于真正的地表輻射值。而在直接校正過程中，由于是從SAR影像空間出發(fā)，尋找正射影像中所對應(yīng)的坐標(biāo)，又由于SAR采用的是側(cè)視方式按距離方向進(jìn)行數(shù)據(jù)采樣的，即按地面目標(biāo)到飛機(jī)平臺的距離進(jìn)行等間隔采樣，而正射校正是要按高度方向在水平面上進(jìn)行等間隔采樣，因此將SAR的側(cè)視方向轉(zhuǎn)換成正垂直水平面的方向，會(huì)造成采樣不均勻，如圖8所示，迎坡方向會(huì)出現(xiàn)欠采樣[15]，表現(xiàn)為空洞，在地形變化很大的區(qū)域，即使經(jīng)過內(nèi)插處理，也沒能得到很好的改善，需要進(jìn)行進(jìn)一步處理，圖12為本文處理后的結(jié)果，可以看出處理方法可以較為正確地反映地表信息，但難免會(huì)有一些不準(zhǔn)確現(xiàn)象。

圖5 模擬SAR影像

圖6 直接校正處理前影像

圖7 直接校正處理前子影像

圖8 直接校正處理前子影像

圖9 間接校正影像

圖10 直接校正處理后影像

圖11 間接校正子影像

圖12 直接校正處理后子影像

圖13 間接校正子影像

圖14 直接校正處理后子影像

在平坦區(qū)域，如圖13和圖14所示，兩種方法均較好地反映了地表信息。由于本文選取的區(qū)域既有平坦區(qū)域，又有崎嶇區(qū)域，在間接校正過程中，模擬SAR與真實(shí)SAR進(jìn)行特征點(diǎn)檢測與匹配時(shí)，即使平坦區(qū)域沒有太多的紋理信息，但在整體環(huán)境的約束下，也得到了不錯(cuò)的校正效果，但當(dāng)平坦區(qū)域范圍很大時(shí)，由于模擬影像紋理信息不明顯，間接校正的方法就難免顯得有些遜色了。

4.3 正射校正精度比較

采用正射校正后的分辨率為4.78 m的光學(xué)影像作為參考基準(zhǔn)，在研究區(qū)均勻選取18個(gè)明顯地物點(diǎn)作為檢驗(yàn)點(diǎn)，分布情況如圖15所示，在校正的SAR影像上和衛(wèi)星影像圖上都能準(zhǔn)確識別。分別讀取每一個(gè)檢驗(yàn)點(diǎn)在校正SAR影像上的坐標(biāo)和衛(wèi)星影像上的坐標(biāo)，并計(jì)算兩者的坐標(biāo)差，結(jié)果如表1、表2所示，由于校正影像與內(nèi)插后的DEM分辨率相同，均為15 m，將表1、表2中的距離與方位向坐標(biāo)差換算成以像元為單位，并用圖（圖16、圖17）的形式更直觀地顯示兩種方法檢驗(yàn)點(diǎn)的精度結(jié)果。

圖15 檢驗(yàn)點(diǎn)分布情況

圖16 間接校正精度

圖17 直接校正精度

對比圖16和圖17可以直觀地看出，間接校正的精度整體要較好于直接校正，根據(jù)中誤差公式M=計(jì)算得，間接校正中檢驗(yàn)點(diǎn)的中誤差為0.62像素，直接校正中檢驗(yàn)點(diǎn)的中誤差為0.95像素，其中M為中誤差，N為檢驗(yàn)點(diǎn)的個(gè)數(shù)，Δx為距離向坐標(biāo)差，Δy為方位向坐標(biāo)差。由此可以得出：間接校正結(jié)果較優(yōu)于直接校正，但兩者差別不大。這與理論也相符，由于間接校正有配準(zhǔn)的過程，也相當(dāng)于一次精度優(yōu)化的過程，而直接校正沒有經(jīng)過配準(zhǔn)的過程，因此間接校正結(jié)果會(huì)較優(yōu)于直接校正結(jié)果；至于兩者差別不大，是因?yàn)閮煞N方法采用的是相同的數(shù)據(jù)集、相同的定位模型，且ALOS-2衛(wèi)星有偏航導(dǎo)引姿態(tài)控制，使得衛(wèi)星軌道精度很高，這對于影像質(zhì)量和定位模型精度有直接的影響，因此，兩者差別不會(huì)很大。

表1 間接校正精度檢驗(yàn)表

5 結(jié)束語

直接校正方法是從SAR影像空間出發(fā)，直接計(jì)算所對應(yīng)的真實(shí)地理空間坐標(biāo)，相對來說較為直觀，校正結(jié)果中出現(xiàn)很多“空洞”像元，地形復(fù)雜的區(qū)域尤為明顯，需要對影像進(jìn)行進(jìn)一步處理，平坦區(qū)域可以較好地反映地表信息。間接校正是從真實(shí)地理坐標(biāo)空間出發(fā)，經(jīng)過影像模擬、匹配、校正，生成最終的正射影像，相對來說較為復(fù)雜，校正結(jié)果中沒有“空洞”像元，地形復(fù)雜的區(qū)域能較為真實(shí)地反映地表信息，但當(dāng)平坦區(qū)域范圍很大時(shí)，由于模擬影像紋理信息較弱，進(jìn)行特征點(diǎn)檢測與匹配時(shí)就顯得有些困難，校正精度也難免會(huì)受到影響。從實(shí)驗(yàn)影像校正精度來看，間接校正結(jié)果較優(yōu)于直接校正，但差別不大，兩種方法均可滿足常規(guī)SAR正射校正精度的需求。因此，兩種方法各有利弊，在實(shí)際應(yīng)用中需要根據(jù)實(shí)際情況選擇合適的方法。

由于本文只對ALOS-2 PALSAR-2一種傳感器一景影像做了實(shí)驗(yàn)，可能不具有足夠的說服力，后續(xù)的研究將是對多種傳感器多景影像做實(shí)驗(yàn)，來佐證本文的結(jié)論。

表2 直接校正精度檢驗(yàn)表

[1]Zi W，Yun S，Chou X，et al.Ortho-rectification of highresolution SAR image in mountain area by DEM[C]//2010 18th International Conference on Geoinformatics，2010：1-6.

[2]Yang B，Wang C，Zhang H，et al.Application and comparison of resolving methods in SAR image ortho-rectification based on G.Konecny model-art.no.67871P[C]//Proceedings of SPIE，2007，6787.

[3]Toutin T，Wang H，Chomaz P，et al.Ortho-rectification of polarimetric radarsat-2/RCM data with accurate lidar DSM[C]//2012 IEEE International Geoscience and Remote Sensing Symposium，2012：5880-5883.

[4]尤紅建，向茂生.機(jī)載干涉雷達(dá)數(shù)據(jù)的直接地學(xué)編碼處理[J].現(xiàn)代雷達(dá)，2006（1）：38-40.

[5]尤紅建，丁赤彪，向茂生.機(jī)載高分辨率SAR圖像直接對地定位原理及精度分析[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2005（8）：712-715.

[6]魏鉅杰，張繼賢，黃國滿，等.面向西部測圖困難地區(qū)稀少控制下的TerraSAR-X影像快速正射糾正[J].遙感信息，2010（3）：43-49.

[7]陳爾學(xué)，李增元.分析法和數(shù)值解算法相結(jié)合的星載SAR直接定位算法[J].中國圖象圖形學(xué)報(bào)，2006（8）：1105-1109.

[8]Xu L，Yang W，Pu G.Ortho-rectification of satellite sar image in mountain area by dem[J].Computing Techniques for Geophysical and Geochemical Exploration，2004，26（2）：145-148.

[9]Zhang G，Qiang Q，Zhu X，et al.Ortho-rectification of satellite-borne SAR image based on image simulation[J].Acta Geodetica et Cartographica Sinica，2010，39（6）：554-560.

[10]黃志楊，曹永鋒.基于影像模擬的多山地區(qū)SAR影像正射校正[J].現(xiàn)代計(jì)算機(jī)，2015（24）：21-25.

[11]王慶，曾琪明，焦健，等.基于DEM和查找表的高分辨率機(jī)載SAR圖像正射校正[J].測繪通報(bào)，2013（11）：10-13.

[12]張波，張紅，王超，等.一種新的星載SAR圖像定位求解方法[J].電波科學(xué)學(xué)報(bào)，2006（1）：1-5.

[13]張紅敏.SAR圖像高精度定位技術(shù)研究[D].鄭州：解放軍信息工程大學(xué)，2013.

[14]敖祖銳.星載SAR影像模擬及其在正射校正中的應(yīng)用研究[D].北京：首都師范大學(xué)，2014.

[15]尤紅建，丁赤飚，吳一戎.基于DEM的星載SAR圖像模擬以及用于圖像精校正[J].中國空間科學(xué)技術(shù)，2006（1）：65-71.