基于高分辨率全極化SAR圖像的取向角校正方法

孫 翔 宋紅軍 王 宇② 李 寧

①(中國科學院電子學研究所 北京 100190)

②(中國科學院大學 北京 100049)

③(河南大學 開封 475004)

1 引言

極化合成孔徑雷達(Polarimetric Synthetic Aperture Radar, PolSAR)是一種獲取豐富地物散射信息的手段，在軍事和民用領域有著廣泛的應用和研究價值[1]。極化目標分解是極化SAR圖像解譯的一個重要分支，是目標特征參數(shù)反演、目標識別和圖像分類的理論基礎[2]。極化散射矩陣將目標散射的能量特性、相位特性以及極化特性統(tǒng)一起來，相對完整地描述了雷達目標的電磁特性[3]。目標分解理論最早由Huynen提出，它利用極化散射矩陣揭示散射體的物理機理，促進極化信息的充分利用[4]。而后Krogager[5]、vanZly[6]、Cloude[7]和Pottier[8]等人做了大量的研究，為極化分解理論奠定了基礎。

然而，地形目標(散射體)的散射依賴于散射取向、形狀、介電特性、散射機制等。復雜地形表面的散射目標往往是隨機取向的，會引起隨機起伏的回波。隨機取向和隨機分布的散射目標難以分類。不同散射機制的不同取向的散射粒子可能會產(chǎn)生類似的散射；反之，相同散射機制的散射體則在隨機的取向角下會造成不同的散射，導致混亂的分解結果與分類結果[9]。去取向概念的引入是為了減少隨機波動取向的影響，將目標自身的物理特點突出地表現(xiàn)出來[10]。去取向的概念類似于斜坡補償?shù)母拍睿煌氖侨ト∠蜻m用于各種目標，包括目標和目標所在的背景。

2011年，為了將大片植被區(qū)域與城市建筑區(qū)域區(qū)分開來，Yamaguchi提出了使用根據(jù)取向角大小旋轉相干矩陣后，再進行分解的方法[11]。該方法首次將去取向角引入了分解算法中。然而，城區(qū)取向角早在2000年就受到了科學家們的關注，那時就已經(jīng)對方位向斜坡變化的補償方法進行了研究[12]。他們提出了兩種補償方法，并利用反射對稱的概念對估計算法進行了統(tǒng)一的分析[13]。尤其針對建筑物地區(qū)的取向角變化利用了后向散射模型進行了重點的研究[14]。偶次散射模型與奇次散射模型被推廣使用于交叉極化項與非對角項中[15]。

為了研究不同特點的建筑物區(qū)域，2013年陳思偉等人提出了主導極化方向角(Dominant Polarization Orientation Angle,DPOA)的概念[16]。DPOA定義為在某一區(qū)域內(nèi)所有像素點POA分布直方圖的頂點所對應的POA值。在研究了多個不同DPOA值的建筑物區(qū)域去取向結果后，他們發(fā)現(xiàn)DPOA的大小與去取向角后建筑物地物分解的結果有較大的關系。經(jīng)過去取向操作后基于模型的目標分解方法在|DPOA|＜22.5°時較為有效，體散射分量過估的問題得到有效解決，二面角散射分量得到了大幅額度提升，由11%提升到44%。然而，當|DPOA|＞22.5°時，分解結果中體散射分量只有小幅度的減少[16]。即使經(jīng)過了常規(guī)去取向處理，在純粹的建筑區(qū)域中體散射分量仍占主導地位[17]。

分解算法針對取向角問題改進至今，雖然誤分解問題得到了一定的改善，但是由于基于模型的分解算法過于強調(diào)體散射分量，代表建筑物的偶次散射機制依然無法在城區(qū)占主導地位。為了解決這一問題，本文提出了基于高分辨率全極化SAR圖像的取向角校正方法，使用了四川都江堰地區(qū)機載全極化圖像進行算法驗證。

2 極化取向角對極化分解結果的影響

極化方向角(Polarization Orientation Angle,POA)即極化取向角，定義為極化橢圓長軸與水平軸之間的夾角，如圖1所示。

在經(jīng)歷極化基變換后，目標的散射相關信息不會發(fā)生變化。假設將目標散射矩陣變換極化基，變換后的極化基是原本的極化基沿視線旋轉Ψ后得到的，那么新的極化基下的散射矩陣與舊極化基下的散射矩陣之間的關系為：

本文使用的數(shù)據(jù)是中科院電子所使用機載X波段全極化合成孔徑雷達于2009年獲得的四川省都江堰地區(qū)高分辨率數(shù)據(jù)。圖像的斜距分辨率和方位分辨率均為0.5 m。入射角范圍在20°到70°之間變化。圖2(a)展示了該地區(qū)數(shù)據(jù)的Pauli分解結果。

圖中用白色的邊框標出了A, B, C 3個區(qū)域以便進行定量分析，其中A, C都是城市建筑物區(qū)域，B區(qū)域是平原低矮植被區(qū)域，存在少量獨立的房屋建筑。從Pauli分解圖可以看出，有大片的建筑物區(qū)域呈現(xiàn)綠色，反映的是體散射機制的性質(zhì)，只有在C區(qū)域周邊存在一部分較為整塊的紫紅色區(qū)域，反映了建筑物應顯示的二面角散射特性。為了解決這一問題，常用的解決方式是去除使得T33最小的取向角后再重新分解，去除的取向角大小如式(7)所示。

由于計算取向角是使用了相干矩陣，所以去除取向角后進行分解時直接采用了基于相干矩陣的非相干目標分解方法，Yamaguchi分解方法。去取向后的Yamaguchi分解結果如圖2(b)所示。

可以看出經(jīng)過取向角旋轉后的Yamaguchi分解圖像依舊呈現(xiàn)大片綠色區(qū)域，即呈現(xiàn)了體散射特性，這與該區(qū)域的真實地貌不相符。出現(xiàn)該問題的原因在上節(jié)已經(jīng)介紹。在這片區(qū)域中，大多數(shù)的建筑物與飛行方向之間的夾角較大，導致了極化旋轉角度較大。大部分建筑物區(qū)域|DPOA|＞22.5°，使得傳統(tǒng)的去取向角方法不再適用。為了解決這一問題，必須利用圖像的其他特性將建筑物的取向角進行校正。

3 高分辨率圖像POA校正算法

本文針對傳統(tǒng)的去取向角方法無法校正由于取向角導致的分解體散射過估問題進行了研究，提出了新的適用于高分辨率圖像的POA校正算法。算法分為兩個部分：第1個部分利用高分辨率圖像在城區(qū)取向角跳變的現(xiàn)象對需要進行POA校正的區(qū)域進行大致提取；第2部分針對提取出的區(qū)域進行迭代逼近的方法找出新的使得T33最小的POA。完整的算法流程圖如圖3所示。

下面將對算法的兩個部分進行詳細地介紹：

第1部分：大致提取需要進行取向角重校正區(qū)域，將無法通過常用POA計算方法去除取向角影響的區(qū)域選取出來。這個部分利用了POA的隨機性進行研究。POA隨機性的差異是對城市和其他地區(qū)進行分類的適當指標。基于POA的跳變現(xiàn)象，提出一個區(qū)分城市地區(qū)和平原地區(qū)的算法。算法流程如下：

步驟1 將POA分為5個部分。分界線分別為24°, 15°, 3°, –3°, –15°和–24° (圖4(a)所示)；

步驟2 按照傳統(tǒng)方法計算每個像素點的POA大小。將每個像素點根據(jù)POA的大小參照5類進行編號；

步驟3 假設一個像素為參考像素(Reference pixel)A，將其上下左右4個像素的編號與A的編號進行比較。此時定義一個新的跳變參數(shù)(Outburst Parameter, OP)。如果所有4個相鄰像素的編號都與A相鄰或相同，則參考像素A的跳變參數(shù)記為0；否則，記為1。如此計算圖像中每個像素點的跳變參數(shù)(圖4(b)所示)；

步驟4 在圖像中放入一個9×9的窗，將這個窗中跳變參數(shù)為1的像素的個數(shù)記為窗中心像素的異構參數(shù)(Heterogeneous Parameter, HP)；

步驟5 移動窗的位置，再次執(zhí)行步驟4，直到得到所有像素點的異構參數(shù)為止；

步驟6 異構參數(shù)大于10的像素點被歸入需要進行POA重新計算的區(qū)域。

以圖4(b)中標記的A, B兩個像素點作為示范，進一步解釋算法流程。參考像素A呈藍色，根據(jù)圖4(a)被標注為4。參考像素A相鄰的4個像素點分別呈橙、紅、藍、綠色，表示它們分別被標注為1, 5,4, 3。其中左、右兩個像素點的編號與參考像素A的編號相鄰，下面的像素點與參考點的編號相同，位于參考點A上方的像素的編號與A的編號既不相鄰也不相同，所以像素點A相鄰像素中存在不連續(xù)的點，故A的跳變參數(shù)為1。同理，參考點B也按照同樣的步驟進行分析。參考像素B的編號為1，它相鄰4個像素點的編號分別為2, 2, 1, 5。對比圖4(a)可以看出參考像素B的編號與其4個相鄰像素的編號相鄰或者相同，所以像素點B的跳變參數(shù)為0。

跳變參數(shù)的數(shù)值象征了像素點是否為跳變點，異構參數(shù)的大小表征的是像素點的POA隨機性，異構參數(shù)的范圍為[0, 81]。理論上，異構參數(shù)在城市建筑物密集地區(qū)較大，在平坦的草原等地區(qū)較小。圖5展示了都江堰地區(qū)數(shù)據(jù)的異構參數(shù)值。圖中像素點越亮，異構參數(shù)越大。可以看出在城市區(qū)域像素點的異構參數(shù)較大，在平原地區(qū)，高亮度的像素點明顯減少。這與都江堰地區(qū)光學圖像(圖6)展現(xiàn)的地物相符。

在算法的第1部分中存在兩個關鍵點。首先，POA種類數(shù)量的確定十分關鍵，種類的多少，如何劃分會影響到跳變點個數(shù)的多少與跳變的可信程度。如果劃分的種類數(shù)量過少，城市地區(qū)就無法體現(xiàn)POA頻繁變化的特征。如果劃分種類的數(shù)量較多，平坦區(qū)域將顯示出與城市區(qū)域相同的隨機性。因此，組的數(shù)量和它們之間的邊界是至關重要的。圖2(a)中3個區(qū)域分別代表不同的地物，研究這3個區(qū)域的POA分布對確定POA種類與分類原則十分重要。

圖7展示了A, B, C 3個區(qū)域的POA分布情況。其中圖7(b)代表了平原區(qū)域的取向角分布。可以看出，在平原區(qū)域，大多數(shù)像素點的POA都處于–15°到15°之間。因此，POA不屬于該范圍內(nèi)的像素點基本屬于城市建筑區(qū)。代表區(qū)域B的直方圖的頂點位于0°，而在圖7(a)與圖7(c)中，直方圖的頂點分別位于–3°和2°處。因此，POA處于–3°到3°的像素點，多屬于B區(qū)域所代表的類型，而A區(qū)域所代表地物類型的像素點的POA大多處于–15°到–3°之間，C區(qū)域所代表地物類型的像素點的POA大多處于3°到15°之間。綜合以上所述的情況，POA以24°, 15°, 3°, –3°, –15°和–24°為分界線分為5大類是合理的。

算法第1部分中另一個關鍵點是異構參數(shù)的閾值。根據(jù)上述步驟，在理想平原區(qū)域內(nèi)，異構參數(shù)應為0；在另一極端情況下，異構參數(shù)應達到81的最大值。如果閾值設置得太低，會將平原區(qū)域一些較高大且密集的植被像素歸類為城市像素；相反，若閾值設置太高，建筑邊緣的像素將不被算入城市區(qū)域。所以確定閾值非常重要。因此需要對典型城區(qū)和平原區(qū)進行研究，圖8顯示了A, B和C區(qū)域的異構參數(shù)分布情況。

通過比較圖8中的3個直方圖可以看出，B區(qū)域像素的異構參數(shù)多處于0到10之間，而在圖8(a)與圖8(c)中異構參數(shù)相對平均地分布在0到30之間。為了確定閾值，在整個圖像上嘗試了從7到12的異構參數(shù)閾值，并將結果與光學照片(圖6)進行了比較。圖9中給出了7到12的異構參數(shù)閾值時，都江堰地區(qū)圖像的狀態(tài)。考慮到相對集中的大面積城市建筑區(qū)域以及平原地區(qū)的零散分布的建筑物，最終選擇了10作為異構參數(shù)的閾值。

第2部分：在上文中提到當|DPOA|＞22.5°時，使用傳統(tǒng)取向角計算方法獲得的取向角無法使得交叉極化項T33達到最小值。

在使用傳統(tǒng)去取向方法是，無論是相干分解還是非相干分解都無法在網(wǎng)格狀的城市建筑物地區(qū)得到正確的分解結果。如圖2所示，圖像中的大部分建筑物都顯示出典型體散射機制的特點。所以，為了找到一種使得交叉極化項最小化的方法，需要采用迭代逼近法來獲得最適合的取向角。針對每個像素點的具體操作流程如下：

步驟1 在–24°到24°之間以1°的步長逐次對相干矩陣進行旋轉，獲得不同角度旋轉下的相干矩陣；

步驟2 比較根據(jù)不同角度旋轉后的所有相干矩陣的交叉極化T33項，記錄下使得T33最小的兩個角度值標記為α1,α2；

步驟3 將α1,α2之間的角度等分為3份，將兩個等分點記做β1和β2；

步驟4 將初始相干矩陣根據(jù)α1,α2,β1和β2進行旋轉，比較4個矩陣的T33值。將最小的兩個T33值所對應的角度定義為新的α1,α2；

步驟5 返回步驟3進行計算，直到|α1-α2|＜0.1°跳出循環(huán)；

步驟6 定義POAnew。

4 實驗結果分析

通過本文提出的取向角校正方法計算得到的取向角與傳統(tǒng)方法得到的結果在城區(qū)的差別較大，如圖10所示。

使用POAnew進行去取向操作后對都江堰地區(qū)全極化SAR圖像進行分解，結果如圖11所示。與使用傳統(tǒng)算法計算POA后使用Pauli分解(圖2(a))與Yamaguchi(圖2(b))分解算法后得到的分解結果相比，新算法針對體散射過估的城市建筑物區(qū)域進行了POA校正，優(yōu)化了分解結果，有效減輕了分解后城區(qū)顯示體散射特性的現(xiàn)象。

本文提出的校正算法使得分解結果在保持平原、草坪以及植被地區(qū)的體散射機制處于主導地位的同時，將城市建筑物區(qū)域受到的POA影響降到最小，使得城區(qū)的二面角散射分量處于主導地位。圖12展示了A, B, C 3個具有代表性的區(qū)域使用傳統(tǒng)方法和本文方法計算取向角，旋轉后的交叉極化T33項。

3個區(qū)域的T33平均值如表1所示。

表1 A, B, C區(qū)域中交叉極化T33項平均值Tab.1 AverageT33of districts A, B and C after rotation using POAs calculated by the traditional method and the new method

通過圖12和表1對T33的分析可以看出，在經(jīng)過本算法進行取向角校正之后，T33在建筑物區(qū)域(A,B區(qū)域)有明顯的下降，在平原地區(qū)(C區(qū)域)基本保持不變。驗證了本文提出的算法能夠有針對性地對城市區(qū)域的POA進行校正，并獲取最小交叉極化分量的結果。

使用新算法進行取向角校正后的分解結果與傳統(tǒng)算法的分解結果對比如圖13所示。

首先針對城市區(qū)域進行結果分析，區(qū)域A與區(qū)域C是城市區(qū)域。在使用Pauli分解算法時，區(qū)域A完全體現(xiàn)了體散射特性，呈綠色(圖13(a))；在Pauli分解的區(qū)域C中，右半部分的的建筑物呈紫紅色，體現(xiàn)了二面角散射機制的散射特性，而左下角的建筑物呈現(xiàn)綠色(圖13(c))。通過本文提出的POA校正算法后進行分解，區(qū)域A與區(qū)域C都呈現(xiàn)了紅色，即建筑物區(qū)域都正確地顯示出了二面角散射特性(圖13(d)，圖13(f))。

區(qū)域B是典型的平原地區(qū)，在區(qū)域B中存在少數(shù)獨立的建筑，在圖13(b)中可以看出，在Pauli分解結果中，區(qū)域B呈現(xiàn)了綠色；在使用本文提出的POA校正算法后進行分解，區(qū)域B的大部分像素點仍然呈綠色，體現(xiàn)了體散射特性，但是少量的建筑物呈現(xiàn)紅色，準確地表現(xiàn)了建筑物的二面角散射特性(圖13(e))。

圖11(b)展示了使用本文算法進行取向角校正后的Yamaguchi分解結果，通過與圖2(b)的對比可以看出，本文方法使得城市區(qū)域顯示出了較明顯的二面角散射特性，呈紅色。下面對3個特定區(qū)域進行定量分析。表2給出了A, B, C 3個區(qū)域表面散射(Ps)，體散射(Pv)和二面角散射(Pd)占總功率的百分比。從表中可以看出，通過本文提出的POA校正算法后進行分解，3個區(qū)域的二面角散射功率占總功率百分比得到了提升，在建筑物較多的A, C兩個區(qū)域提升幅度較大，變化尤為明顯。

和傳統(tǒng)的POA計算方法進行比較，對所獲得的分解結果進行分析可以看出，本文提出的對高分辨率區(qū)域進行POA校正的算法，能夠使得分解結果更準確。本算法可以有效緩解|DPOA|＞22.5°區(qū)域分解錯誤的問題。

表2 A, B, C區(qū)域中各散射分量占總功率的百分比Tab.2 Percentage of scattering powers with traditional method and the new method

5 結論

基于高分辨率圖像取向角跳變的特性，本文提出了一種適用于城區(qū)體散射過估問題的取向角校正方法。該方法以取向角隨機性為基礎，考慮高分辨率圖像地物分界明確，散射特性清晰對取向角的影響，分兩步對城區(qū)取向角進行校正。本文使用X波段全極化機載SAR實驗數(shù)據(jù)對算法進行了驗證，獲得了與真實地貌相符的分解結果，有效緩解了城區(qū)體散射分量占主導的問題，驗證了算法的有效性。