基于極化方向角的地形坡度的估計*

袁萬里，張夢杰

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083)

基于極化方向角的地形坡度的估計*

袁萬里，張夢杰

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083)

SAR數(shù)據(jù)處理中，地形坡度信息是重要的地形參數(shù)。利用SAR影像提取地形坡度的主要方法分為雷達(dá)干涉測量(InSAR)和雷達(dá)極化(PolSAR)測量，其中雷達(dá)極化測量包含正交雙軌道全極化SAR測量和單軌道全極化SAR測量。文章結(jié)合極化方向角及Lambertian后向散射模型得到的距離向坡度與方位向坡度之間的關(guān)系，從而反演方位向與距離向坡度。實驗選取了覆蓋青海省治多縣的兩景全極化ALOS-PALSAR數(shù)據(jù)進(jìn)行實驗。結(jié)果表明：反演得到的地形坡度與利用InSAR技術(shù)得到的地形坡度能較好地吻合。

極化方向角；Lambertian后向散射模型；距離向坡度；方位向坡度

0 引言

20世紀(jì)50年代就出現(xiàn)的合成孔徑雷達(dá)(SAR)是一種主動式的微波傳感器，即傳感器本身向地面發(fā)射能量波束并接收其散射回波信號進(jìn)行探測[1]。它不依賴于光照條件，具有全天時、全天候的對地觀測能力，同時它分辨率較高，已經(jīng)成為當(dāng)前遙感研究的一個重要領(lǐng)域。地形坡度信息在生產(chǎn)科研中有很廣泛的應(yīng)用。向茂生等[2]發(fā)現(xiàn)干涉相位中殘余點分布受地形坡度影響較大，經(jīng)過理論分析與模擬實驗，得出地形變化越大，殘余點分布越密集的結(jié)論。Schuler等[3]發(fā)現(xiàn)地形坡度和散射回波極化方向角存在幾何關(guān)系，提出了求解極化方向角的相應(yīng)算法，并對德國Villingen Schwenningen地區(qū)進(jìn)行了成像實驗。王敏等[4]根據(jù)已有的DEM數(shù)據(jù)仿真了正交雙航過全極化SAR圖像，驗證了正交雙航過全極化SAR測量坡度的可行性。Jin等[5]利用傾斜地表水平方位排列產(chǎn)生的圖像紋理進(jìn)行形態(tài)學(xué)細(xì)化算法來確定方位向坡度角，這樣只需要單次飛行的SAR圖像數(shù)據(jù)，并利用中國廣東省惠州地區(qū)的SIR-C數(shù)據(jù)進(jìn)行了實例驗證。Li等[6]基于極化方向角補償前后散射能量比值得到方位向的坡度值，同時利用極化方向角和方位向坡度及距離向坡度之間的幾何關(guān)系得到距離向的坡度值，利用多重網(wǎng)格算法得到地形信息，通過和干涉SAR獲取的地形信息相互比較，驗證了此方法的有效性。綜上所述，國內(nèi)外許多學(xué)者已經(jīng)利用星載或者機(jī)載SAR進(jìn)行了大量極化雷達(dá)提取地形信息的研究。本文結(jié)合極化方向角及Lambertian后向散射模型，根據(jù)地形坡度與極化方向角之間的關(guān)系聯(lián)立得到距離向坡度與方位向坡度，并以2007年5月30日獲取的覆蓋青海省治多縣的ALOS-PALSAR影像為實驗數(shù)據(jù)，進(jìn)行了極化方向角反演地形坡度的研究，并將計算結(jié)果與實驗區(qū)的InSAR獲得的結(jié)果進(jìn)行比較。

1 極化SAR反演地形坡度的方法

全極化SAR圖像，可以得到觀測區(qū)域2×2散射矩陣S2，其表達(dá)式為：

(1)

式中：下標(biāo)h，v為極化通道。

當(dāng)散射矩陣滿足互易性原理時(此時滿足shv=svh)，為了更好地從散射矩陣中提取信息，會將散射矩陣轉(zhuǎn)變?yōu)橄喔删仃嚮騾f(xié)方差矩陣，其表達(dá)式為：

(2)

(3)

式中：*表示復(fù)共軛，<·> 表示整體平均(集平均)。

為了抑制噪聲對實驗結(jié)果的影響，對相干矩陣T3與協(xié)方差矩陣C3進(jìn)行多視濾波處理(方位向視數(shù)為6，距離向視數(shù)為2，濾波采用窗口大小位3×3的增強Lee濾波)。

極化SAR數(shù)據(jù)提取地形坡度中，最重要的參數(shù)為極化方向角(POA)。極化SAR圖像的同極化和交叉極化后向散射是電磁波極化橢圓率角ε和極化方向角φ的函數(shù)。當(dāng)?shù)孛嫫教箷r，后向散射信號一般在極化方向角為零時達(dá)到極大值；當(dāng)?shù)孛鎯A斜時，同極化和交叉極化后向散射信號的極大值會偏離零方向角位置。因此，極化方向角φ的變化可用來推算地表坡度。文獻(xiàn)[3]中給出了極化方向角與地形坡度間的關(guān)系：

(4)

式中：φ為極化方向角；ω為方位向(航向)的坡度角；γ為距離向的(斜距方向)的坡度角；η為雷達(dá)入射角，為已知值。

式(4)中，極化方向角的值可以通過數(shù)據(jù)本身獲得，獲取地形坡度則需要額外的條件方程。

式(4)中極化方向角可以通過極化特征法[3]和圓極化法[7]獲得。圓極化法在實際應(yīng)用中得到了較多的運用，然而其受數(shù)據(jù)本身質(zhì)量的影響較大，得到的極化方向角往往與真實的極化方向角的值相差較大。因此，本文利用已有的DEM(如SRTM DEM)數(shù)據(jù)估算極化方向角作為參考值并結(jié)合Lambertian后向散射模型反演地形坡度。

式(4)中一個方程包含兩個未知數(shù)ω與γ，因此求解兩個未知數(shù)需要額外的方程，本文運用Lambertian后向散射模型得到地形坡度間的幾何關(guān)系，原理如下：

已知，SAR是沿已知路徑運動，以一定間隔發(fā)送脈沖并接收目標(biāo)的散射回波[7]。其散射強度取決于地形坡度和地物類型，Lambertian后向散射模型影像與地面的散射強度表達(dá)式[8]為：

(5)

式中：K為SAR系統(tǒng)配準(zhǔn)常數(shù)；σ0為后向散射系數(shù)；Rω、Rγ為影像方位向和距離向的分辨率；散射強度I(ω，γ)的值為相干矩陣對角矩陣元素的和。當(dāng)?shù)孛嫱耆教箷r(此時ω=0，γ=0)，地面的散射強度表達(dá)式為：

I0=I(0，0)=Kσ0RωRγsin2η

(6)

當(dāng)?shù)孛嫱耆教箷r，散射矩陣的交叉項元素接近于零，因此在沒有先驗信息的情況下，選擇哪些散射矩陣交叉項元素的模接近為零的散射單元作為一個集合，將它們散射強度的均值作為I0的值。

比值R(ω，γ)定義為式(5)與式(6)的比值，其表達(dá)式為：

(7)

式(7)可以寫為：

(8)

其中

(9)

為了便于運算，式(4)可以寫為：

ω=arctan[tanφ(sinη-cosηtanγ)]

(10)

將式(7)和式(10)聯(lián)立，對ω賦初值，并將初值代入式(9)中得到R1的值，將R1代入式(8)得到γ，再代入式(10)得到新的ω，然后以此ω為初值代入式(9)重新進(jìn)行運算，直到滿足abs(ωn+1-ωn)

2 實驗區(qū)概況及數(shù)據(jù)處理

實驗區(qū)位于青海省治多縣的西部，地形地貌以高原為主，平均海拔在4 500 m以上，其境內(nèi)雪山連綿，冰川廣布，河流縱橫交織，湖泊星羅棋布，處于“中華水塔”的核心，是中國長江流域乃至南亞及東南亞國家的“生態(tài)源”，生態(tài)地位極為突出。本文選取的研究范圍為35°21′N～36°14′N，91°48′E～92°23′E。

實驗數(shù)據(jù)為2007年4月14日和2007年5月30日獲取的2景ALOS-PALSAR L波段的全極化單視復(fù)數(shù)影像，其距離向分辨率為4.684 m，方位向分辨率為3.550 m，入射角為23.490°。以2007年5月30日影像為極化方向角地形坡度估計所用的數(shù)據(jù)，2007年4月14日的影像與其形成干涉對，利用InSAR技術(shù)獲取地形坡度信息以驗證基于極化方向角反演的地形坡度的準(zhǔn)確性?；跇O化方向角的地形坡度的估計主要是在Matlab軟件下進(jìn)行操作的，利用InSAR技術(shù)獲取地形坡度信息則是在Gamma軟件下運行的，其中，后者處理流程，如圖1所示。

圖1 利用InSAR技術(shù)獲取地形坡度信息的流程圖Fig.1 Flow of obtaining terrain slope by using InSAR technology

圖1中在獲得SAR影像對為兩景影像HH通道數(shù)據(jù)，首先要將兩個影像精確的配準(zhǔn)，要求配準(zhǔn)精度優(yōu)于0.2像元，保證兩景影像空間重疊，光譜相似，本文中所選取的影像對的配準(zhǔn)精度為方位向0.08，距離向為0.12。影像配準(zhǔn)后將兩幅復(fù)數(shù)影像共軛相乘產(chǎn)生干涉相位圖，此時干涉相位中含有必要的振幅和相位信息，同時也有各種噪聲影響，因此生成的干涉相位需要進(jìn)行相位噪聲濾除處理。接著是去平地相位，又可以稱作去平地效應(yīng)。所謂的平地效應(yīng)是指生成干涉相位時，該干涉相位既包含平坦地面產(chǎn)生的干涉條紋又包含地形起伏所引起的干涉條紋，增加了隨后相位解纏的難度，所以需要計算消除平地效應(yīng)所對應(yīng)的相位分量。然后進(jìn)行InSAR技術(shù)的關(guān)鍵環(huán)節(jié)，即相位解纏。因為某些高程不同的地物目標(biāo)到SAR衛(wèi)星天線的斜距相差為2π的整數(shù)倍，要獲取各地物目標(biāo)的真實高程信息，就要恢復(fù)各等相位差的地物目標(biāo)的真實相位值。之后進(jìn)行高程計算，就是結(jié)合衛(wèi)星軌道參數(shù)和地面控制點等已有信息，將解纏后的相位轉(zhuǎn)換為相應(yīng)的高程。最后將生成的高程信息進(jìn)行差分處理，獲取研究區(qū)域的地形坡度信息。

實驗區(qū)的Google衛(wèi)星圖及2007年5月30日獲取的極化SAR數(shù)據(jù)的Pauli RGB合成圖，如圖2所示。

圖2 研究區(qū)域的衛(wèi)星圖和Pauli RGB影像Fig.2 Google image and Pauli RGB image in study area

3 實驗結(jié)果與分析

為了驗證實驗結(jié)果，將2007年5月30日影像與2007年4月14日影像進(jìn)行干涉處理(其中前者為主影像，后者為從影像)，經(jīng)過圖1所示的流程獲取研究區(qū)域的地形坡度信息。其中，基于極化方向角的地形坡度值的估計顯示在圖3(a)、4(a)中，利用InSAR獲取的地形坡度值顯示在圖3(b)、4(b)中。

圖3 極化方向角反演獲取與InSAR獲取的方位向坡度估值的比較Fig.3 The comparison between valuations of azimuth slope based on POA and InSAR

圖4 極化方向角反演獲取與InSAR獲取的距離向坡度估值的比較Fig.4 The comparison between valuations of range slope based on POA and InSAR

在圖3(b)、4(b)中選取圖中下方四邊形剖面曲線并將其與圖3(a)、4(a)中同一位置的參考數(shù)據(jù)進(jìn)行對比，結(jié)果如圖5所示。

圖5 方位向坡度剖線與距離向坡度剖線Fig.5 The profile lines of azimuth slope and range slope

從圖3、圖4中獲取的坡度值的像素分布來看，方位向坡度估計與參考數(shù)據(jù)地形趨勢比較接近，而距離向坡度結(jié)果在局部區(qū)域與參考數(shù)據(jù)相差較大，選取比較平坦的區(qū)域(圖中上方四邊形)和地形起伏較大的區(qū)域(圖中下方四邊形)與干涉獲取的坡度(參考坡度)進(jìn)行比較，得到的結(jié)果為：方位向坡度整體區(qū)域的RMSE為3.80°，平坦區(qū)域的RMSE為1.21°，地形起伏較大區(qū)域的RMSE為6.56°；距離向坡度整體區(qū)域的RMSE為8.38°，平坦區(qū)域的RMSE為5.24°，地形起伏較大區(qū)域的RMSE為12.3°。

究其原因，由式(4)知，方位向坡度主要取決于極化方向角，而極化方向角的主要誘因就是地形在方位向的變化，因此方位向坡度與參考數(shù)據(jù)相差較?。欢嚯x向坡度則主要憑借雷達(dá)圖像中的灰度信息來計算，其除了受雷達(dá)入射角度和距離向坡度影響以外還受地物自身屬性的影響，因此距離向坡度較參考數(shù)據(jù)有一定的差距。

4 結(jié)論

本文對基于極化方向角反演地形坡度進(jìn)行了研究與分析，通過Lambertian后向散射模型得到距離向坡度與方位向坡度之間的關(guān)系，然后再與極化方向角和地形坡度的關(guān)系相結(jié)合得到地形坡度。經(jīng)過驗證，方位向坡度與參考數(shù)據(jù)較為一致，而距離向坡度與參考數(shù)據(jù)有一定的差距。

[1] 廖明生，林 琿.雷達(dá)干涉測量——原理與信號處理基礎(chǔ)[M].北京：測繪出版社，2003.

[2] 向茂生，李樹楷.地形坡度對干涉相位殘余點發(fā)布的影響[J].測繪學(xué)報，1999，28(2)：139-143.

[3] Schuler，Lee，etal.Measurement of topography using polarimetric SAR images[J].IEEE Trans on Geosciences and Remote sensing，1996，34(5)：1 266-1 277.

[4] 王敏，董臻，梁甸農(nóng).正交雙航過全極化SAR圖像測高研究[J].信號處理，2009，25(1)：48-51.

[5] Jin Yaqiu，Luo Lin.Terrain topographic inversion using single-pass polarimetric SAR image data[J].Science in China Ser.F Information Sciences，2004，47(4)：490-500.

[6] Li Yang，Wen Hong，etal.Topography Retrieval From Single-Pass POLSAR Data Based on the Polarization-Dependent Intensity Ratio [J].IEEE Trans on Geosciences and Remote sensing，2015，53(6)：3 160-3 177.

[7] J.S.Lee，D.L.Schuler and T.L.Ainsworth.Polarimetric SAR Data compensation for Terrain Azimuth Slope Variation[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38(5)：2 153-2 163.

[8] Xi Chen，Chao Wang，etal.DEM Generation Combining SAR Polarimetry and Shape-From-Shading Techniques[J].IEEE Geoscience and Remote Sensing Letters，2008，1(6)：28-32.

[9] Lee J S，Schuler D L，Ainsworth T L，etal.On the estimation of radar polarization orientation shifts induced by terrain slopes[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40(1)：30-41.

[10] J.S.Lee，Dale L.Schuler，T.L.Answorth，etal.A Review of Polarization Orientation Angle Estimation and Applications[C].EUSAR 2006，Dresden，Germany，May 2006.

[11] 金旭晨，韓震，劉瑜，等.基于全極化SAR的浙江衢山島DEM 信息提取[J].海洋學(xué)研究，2015，33(1)：33-38.

[12] 梁淮寧.全極化SAR地形高度測量技術(shù)研究[D].四川：電子科技大學(xué)，2001.

[13] Dale L，SCHULER DL，etal.Compensation of Terrain Azimuthal Slope Effects in Geophysical Parameter Studies Using Polarimetric SAR Data [J].REMOTE SENS.ENVIRON，1999，69(4)：139-155.

[14] D.L.Schuler，J.S.Lee，etal.Measurement of ocean surface slopes and wave spectra using polarimetric SAR image data[J].Remote Sensing of Environment，2004，91(2)：198- 211.

Estimation of Terrain Slope Based on Polarization Orientation Angle

YUAN Wan-li，ZHANG Meng-jie

(SchoolofGeosciencesandInfo-Physics，CentralSouthUniversity，ChangshaHunan410083，China)

Terrain slope information is an important topographic parameter.The main methods for extracting terrain slope by using SAR image include Interferometric SAR(InSAR)and Polarimetric SAR(PolSAR).Where PolSAR method includes two orthogonal track fully polarimetric SAR and single-track fully polarimetric SAR.In this paper，we use Polarimetric orientation angle combining with Lambertian backscattering model to inverse the azimuth slope and the range slope.We select two full polarimetric ALOS-PALSAR images covered Zhiduo County，Qinghai，China for experiments.The results show that the inversed terrain slopes are consistent with the slopes abtained by using InSAR method.

Polarization orientation angle；Lambertian backscattering model；azimuth slope；range slope

2017-03-08

P 236； TN 957.52

A

1007-9394(2017)02-0013-04

袁萬里(1990～)，男，河南睢縣人，碩士，主要研究方向：雷達(dá)遙感對地觀測。