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

孫 翔 宋紅軍 王 宇② 李 寧

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

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

③(河南大學(xué) 開封 475004)

1 引言

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

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

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

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

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

2 極化取向角對極化分解結(jié)果的影響

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

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

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

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

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

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

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

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

下面將對算法的兩個部分進(jìn)行詳細(xì)地介紹：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟6 定義POAnew。

4 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

表1 A, B, C區(qū)域中交叉極化T33項(xiàng)平均值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)過本算法進(jìn)行取向角校正之后，T33在建筑物區(qū)域(A,B區(qū)域)有明顯的下降，在平原地區(qū)(C區(qū)域)基本保持不變。驗(yàn)證了本文提出的算法能夠有針對性地對城市區(qū)域的POA進(jìn)行校正，并獲取最小交叉極化分量的結(jié)果。

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

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

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

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

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

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

5 結(jié)論

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