朱言江韓 震和思海

(1.上海海洋大學(xué) 海洋科學(xué)學(xué)院,上海201306；2.上海市河口海洋測(cè)繪工程技術(shù)研究中心,上海201306)

利用全極化合成孔徑雷達(dá)(Full-polarization Synthetic Aperture Radar)數(shù)據(jù)可以獲得研究對(duì)象的全散射矩陣,該矩陣包含了研究對(duì)象的重要散射信息、空間結(jié)構(gòu)信息和紋理信息,利用這些信息可以深入了解研究對(duì)象的物理及幾何特性[1-3]。目前,全極化SAR數(shù)據(jù)的極化特征提取工作已經(jīng)取得了一系列的研究成果,經(jīng)過數(shù)十年的努力,國(guó)內(nèi)外學(xué)者在不同理論的基礎(chǔ)上,提出多種方法進(jìn)行極化特征分解。如Cloude提出的H/Alpha分解[4]、Freeman等提出的Freeman-Durden分解[5]、Yamaguchi等提出的Yamaguchi分解[6]。對(duì)于不同的地物類型,有不同的分解方法,例如:Kiran等[7]利用不同分解方法對(duì)城市地區(qū)進(jìn)行分類研究,結(jié)果表明,在城市居民點(diǎn)螺旋散射成分更多,Yamaguchi分解效果良好,其次是Freeman-Durden分解,而在水體等表面散射的區(qū)域Pauli分解效果更好。陳軍等[8]以ALOS衛(wèi)星的PALSAR影像為研究數(shù)據(jù),提出一種綜合Pauli極化特征分解和支持向量機(jī)(SVM)的分類策略,研究結(jié)果表明Pauli-SVM算法可以有效地提高分類的準(zhǔn)確性。馬騰等[9]采用Radarsat-2全極化數(shù)據(jù),利用極化目標(biāo)分解方法提取到了散射熵、平均散射角、反熵、平均特征值、單次反射特征值相對(duì)差異度和二次反射特征值相對(duì)差異度,結(jié)合實(shí)測(cè)數(shù)據(jù),分析了各參數(shù)對(duì)于耕地、裸地、含植被水體和建筑等類別的可分離性。Geng等[10]利用Radarsat-2數(shù)據(jù)的目標(biāo)分解分析了在不同的環(huán)境條件下中國(guó)江蘇沿岸潮間帶的極化特征和散射機(jī)制,定量比較分析了Radarsat-2數(shù)據(jù)分解極化參數(shù)對(duì)潮間帶中不同物體的極化特性,表明潮汐對(duì)灘涂和水產(chǎn)養(yǎng)殖場(chǎng)的主要散射機(jī)制具有不同的影響。

如何有效地選擇極化特征和分類算法是全極化SAR影像分類領(lǐng)域亟待解決的問題。經(jīng)過反復(fù)試驗(yàn),同時(shí)結(jié)合長(zhǎng)江口九段沙研究區(qū)域地物類型,利用Pauli分解的極化散射參數(shù)可以準(zhǔn)確地獲取地物實(shí)際物理特性相對(duì)應(yīng)的極化特征。綜合不同提取算法和分類算法的基礎(chǔ),本文選取2015-04-05 Radarsat-2全極化SAR數(shù)據(jù),分析Pauli極化分解出的散射參數(shù)對(duì)九段沙地物散射特性的影響,并利用區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)提取潮溝信息。

1 研究區(qū)域

九段沙(圖1)位于長(zhǎng)江入?？谔?121°46'～122°15'E,31°03'～31°17'N),是上海市的重要濕地之一。由于受潮流的影響,九段沙一直處于變化中,其植物群落主要有3種,分別為互花米草、蘆葦和海三棱藨草[11]。

圖1 研究區(qū)域Fig.1 Geographic position of the study area

2 研究方法

首先對(duì)Radarsat-2全極化數(shù)據(jù)進(jìn)行預(yù)處理,包括數(shù)據(jù)多視處理和相干斑濾波；然后對(duì)得到的數(shù)據(jù)進(jìn)行Pauli極化相干分解,分別得到3種散射機(jī)制系數(shù)K1,K2和K3,并且將這3個(gè)系數(shù)轉(zhuǎn)化為散射強(qiáng)度；隨后統(tǒng)計(jì)分析不同地物類型散射強(qiáng)度值；最后選取合適的散射機(jī)制,利用區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)中的形態(tài)膨脹腐蝕和去除短枝處理提取九段沙潮溝信息(圖2)。

圖2 九段沙潮溝信息提取流程Fig.2 Flow chart of the tidal channel extraction

2.1 數(shù)據(jù)預(yù)處理

對(duì)獲取到的Radarsat-2全極化數(shù)據(jù)進(jìn)行多視處理,比例為4∶2,之后對(duì)多視處理后的圖像采用精極化Lee濾波進(jìn)行相干斑濾波,對(duì)比其結(jié)果,發(fā)現(xiàn)5×5窗口對(duì)降低相干斑噪聲的效果較好,同時(shí)九段沙地物和潮溝的邊緣信息保持得比較清晰完整。

2.2 全極化SAR數(shù)據(jù)的Pauli分解

全極化SAR數(shù)據(jù)目標(biāo)分解主要是基于極化數(shù)據(jù)的散射矩陣來揭示不同散射體的物理表現(xiàn)形式,這樣可以更好地認(rèn)識(shí)和利用極化信息來提取目標(biāo)信息。其優(yōu)勢(shì)是可以將復(fù)雜地物回波信號(hào)分解為多個(gè)簡(jiǎn)單的標(biāo)準(zhǔn)目標(biāo)回波信號(hào)之和,通過分析不同散射機(jī)制下的各分量獲得原目標(biāo)的信息[12]。由于極化分解技術(shù)與獲取目標(biāo)信號(hào)中的相位有關(guān),Pauli目標(biāo)分解的散射特征參數(shù)能夠準(zhǔn)確獲得與其地物物理特性相對(duì)應(yīng)的極化特征,同時(shí)可在一定程度上抑制SAR成像系統(tǒng)自身產(chǎn)生的相干斑噪聲,即使在有噪聲或去極化效應(yīng)的情況下,仍然可以用它進(jìn)行分解,這樣很好地保持了目標(biāo)信息的完整性。

Pauli分解是基于極化散射矩陣S,通過定義不同的極化基矩陣提取極化特征的方法,不同的極化基矩陣代表不同的地物類型[13]。

首先定義基本的散射矩陣,稱為Pauli基{S a,S b,S c,S d},公式為

式中,a為單次散射強(qiáng)度；b為入射角0°的偶次散射強(qiáng)度；c為入射角45°的偶次散射強(qiáng)度；d為S矩陣中反對(duì)稱組分的強(qiáng)度；表示奇次散射,指表面散射,如地面等；表示偶次散射,指二面角散射,如建筑物和樹干等表示π/4偶次散射,指體散射,如灌木和樹冠等；表示不存在的地物類型,因此d取值為零[8]。

基于已定義的Pauli基,對(duì)于本文將研究的散射矩陣,式(1)可以寫成:

式中,H表示水平方向,V表示垂直方向,SVH表示水平極化發(fā)射、垂直極化接收的回波強(qiáng)度?？梢詫懗上蛄縆的形式:

當(dāng)介質(zhì)滿足互易條件時(shí),S HV=S VH,公式(3)變?yōu)?/p>

Pauli相干極化分解將極化散射矩陣S反映到3個(gè)基本散射類型上,分別為奇次散射、繞軸旋轉(zhuǎn)0°的二面角散射和繞軸旋轉(zhuǎn)45°的二面角散射。

經(jīng)過數(shù)據(jù)預(yù)處理后的Pauli分解結(jié)果表明(圖3):在平坦表面的區(qū)域,如河流、潮溝等表現(xiàn)得相對(duì)平滑,呈現(xiàn)出奇次散射,Pauli分解的結(jié)果呈現(xiàn)藍(lán)色,這是因?yàn)槠娲紊⑸涫撬虻闹饕⑸漕愋?；在陸地上由于受到低矮植被的影?絕大部分呈現(xiàn)綠色,主要為二面角散射,植被主要分布在潮溝附近；九段沙的中沙和下沙的中心區(qū)域呈現(xiàn)偏紅色,這是由于這兩處區(qū)域受到高大蘆葦?shù)挠绊?呈現(xiàn)為45°的二面角散射。

潮溝、水體、低矮植被和蘆葦是九段沙典型的4種地物類型,為了更好地分析它們的物理散射特性,將Pauli分解的3個(gè)系數(shù)轉(zhuǎn)化成散射強(qiáng)度(圖4)。每種地物類型選擇2個(gè)區(qū)域,共8個(gè)區(qū)域(圖4c),每個(gè)區(qū)域內(nèi)選擇40個(gè)采樣點(diǎn)。統(tǒng)計(jì)分析這4種典型地物的3種散射強(qiáng)度值,結(jié)果見表1?？梢钥闯?蘆葦?shù)钠娲紊⑸鋸?qiáng)度平均值遠(yuǎn)大于其他3種地物類型的奇次散射強(qiáng)度平均值,潮溝和水體的奇次散射強(qiáng)度值較為接近,分別為0.040和0.038,所以,其可以被看作是平靜的表面。相對(duì)于潮溝和水體,低矮植被的奇次散射強(qiáng)度值并沒有太大的變化。在九段沙典型類型地物中,各種散射機(jī)制分布較為復(fù)雜,其中植被的二面角散射成分更多些,這主要是由于九段沙低矮植被與地面形成散射的結(jié)果。從表1也可以看出,研究目標(biāo)潮溝的二面角散射強(qiáng)度平均值為0.007,而低矮植被和蘆葦?shù)亩娼巧⑸鋸?qiáng)度平均值分別為0.783和0.104,相較于奇次散射和45°二面角散射強(qiáng)度值,潮溝和植被的強(qiáng)度平均值差異較大。因此,Pauli分解中二面角散射系數(shù)更能體現(xiàn)出潮溝和植被之間的差異性。

圖3 Pauli分解結(jié)果的假彩色圖Fig.3 The false color composite image of Pauli decomposition

圖4 Pauli分解的散射強(qiáng)度Fig.4 The scattering intensity image of Pauli decomposition

表1 不同地物類型的散射強(qiáng)度值Table 1 The statistics of the scattering intensities of different land types

2.3 基于區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)的潮溝信息提取

為了進(jìn)一步分析分解參數(shù)的選擇對(duì)潮溝信息的提取效果,我們?cè)赑auli分解的二面角散射強(qiáng)度圖像上選取具有代表性的3條潮溝(圖4b中A,B和C)作為研究對(duì)象,基于區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)提取潮溝信息。

2.3.1 區(qū)域生長(zhǎng)法

區(qū)域生長(zhǎng)法的基本原理是:在每個(gè)分割的區(qū)域找種子像素作為生長(zhǎng)的起點(diǎn),再將種子像素周圍鄰域中與種子像素有相同或相似性質(zhì)的像素(根據(jù)某種事先確定的生長(zhǎng)或相似準(zhǔn)則來判斷)合并到種子像素所在的區(qū)域中,然后將這些新像素當(dāng)作新的種子像素繼續(xù)進(jìn)行上述過程,直到再?zèng)]有滿足條件的像素可被包括進(jìn)去[14]。

區(qū)域生長(zhǎng)法提取潮溝信息的步驟主要有5步[15]:

1)輸入原始圖像,選擇初始種子點(diǎn)并標(biāo)記為潮溝點(diǎn),并對(duì)原始圖像進(jìn)行掩膜；

2)取出種子點(diǎn),計(jì)算其與鄰域中其他像素的灰度差；

3)將灰度差小于給定閾值的像素點(diǎn)與種子點(diǎn)合并,繼續(xù)生長(zhǎng)并且覆蓋掩膜；

4)閾值分割；

5)輸出潮溝提取結(jié)果。

在運(yùn)用區(qū)域生長(zhǎng)法提取潮溝信息過程中,不同潮溝的復(fù)雜程度不一樣,選取種子點(diǎn)的個(gè)數(shù)也是不同的。潮溝A,由于土壤出露,因此選擇2個(gè)種子點(diǎn)；潮溝B,在2個(gè)分叉分別選擇1個(gè)種子點(diǎn),共計(jì)2個(gè)種子點(diǎn)；潮溝C,由于空間結(jié)構(gòu)較復(fù)雜,則在每級(jí)分叉都選擇的潮溝選擇1個(gè)種子點(diǎn),共選擇3個(gè)種子點(diǎn)。

2.3.2 數(shù)學(xué)形態(tài)學(xué)

數(shù)學(xué)形態(tài)學(xué)主要是以膨脹、腐蝕為基礎(chǔ),使用這2種簡(jiǎn)單基礎(chǔ)的運(yùn)算來推導(dǎo)其他較為復(fù)雜的運(yùn)算以實(shí)現(xiàn)圖像處理,具有計(jì)算快、算法較為簡(jiǎn)單的特點(diǎn)。對(duì)潮溝信息進(jìn)行數(shù)學(xué)形態(tài)學(xué)處理可以改善潮溝信息在提取時(shí)發(fā)生的斷裂情況,并且可以降低潮溝信息中的噪聲。4種基本的算子包括膨脹和腐蝕運(yùn)算、結(jié)合膨脹和腐蝕的開運(yùn)算和閉運(yùn)算。

1)膨脹和腐蝕,記M為要處理的圖像,N為結(jié)構(gòu)元素,x為圖像M中的一個(gè)點(diǎn):

式中,Φ為空集。

2)結(jié)合膨脹和腐蝕的其他運(yùn)算主要有:

選取的3條潮溝的發(fā)育形態(tài)和地理位置都不相同,提取這3條潮溝的信息,可以檢驗(yàn)Pauli分解出的二面角散射機(jī)制對(duì)整體潮溝的提取效果。首先以潮溝B為對(duì)象利用該方法提取信息,提取效果如圖5所示。采用同樣的方法對(duì)潮溝A和潮溝C提取信息,提取效果如圖6所示。

從圖6可以看出,利用區(qū)域生長(zhǎng)法對(duì)經(jīng)過Pauli分解后的二面角散射系數(shù)進(jìn)行潮溝A,B和C提取的總體效果較好,但是在提取潮溝信息的過程中,有些潮溝有斷裂現(xiàn)象,并且出現(xiàn)了一些非潮溝,如圖5b和圖6中圓圈標(biāo)記出來的潮溝顯示,它們并沒有生長(zhǎng)到一起。這主要是因?yàn)樵搮^(qū)域原始圖像潮溝邊緣受植被的遮擋或者由潮溝中的土壤出露影響所致。另外,由于區(qū)域生長(zhǎng)法的閾值的設(shè)置,也會(huì)導(dǎo)致一些非潮溝信息被提取出來。針對(duì)區(qū)域生長(zhǎng)法提取過程中出現(xiàn)的問題,利用數(shù)學(xué)形態(tài)學(xué)進(jìn)行了進(jìn)一步的處理。

圖5 潮溝B原始圖像及其提取效果圖Fig.5 Raw image and the extracted result of tidal channel B

圖6 潮溝A和C信息提取效果Fig.6 Extracted tidal channels(A and C)

對(duì)于潮溝斷裂現(xiàn)象,考慮到潮溝在圖像上呈現(xiàn)出形態(tài)細(xì)長(zhǎng)和曲率變化較大的特點(diǎn),選取了尺寸較小的結(jié)構(gòu)元素進(jìn)行膨脹腐蝕處理。這樣不僅可以使斷裂的潮溝得到連接,而且不會(huì)導(dǎo)致提取的潮溝長(zhǎng)度和寬度發(fā)生變化。本文對(duì)選取的潮溝采用3×3正方形結(jié)構(gòu)元素(圖7)進(jìn)行先膨脹后腐蝕處理,直到斷裂潮溝被連接上為止。

由于區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)分別對(duì)潮溝提取和斷裂連接過程中會(huì)產(chǎn)生非潮溝像素點(diǎn),因此在此基礎(chǔ)上利用spur函數(shù)進(jìn)行去除毛刺處理,函數(shù)值根據(jù)實(shí)驗(yàn)要求進(jìn)行設(shè)定,最后得到較為完整的潮溝信息,利用同樣的方法提取了潮溝B和潮溝C的信息(圖8)。

圖7 正方形結(jié)構(gòu)元素Fig.7 Square structural elements

圖8 潮溝A,B和C信息提取效果Fig.8 Extracted tidal channels(A,B and C)

3 結(jié) 語

通過對(duì)長(zhǎng)江口九段沙地區(qū)SAR影像進(jìn)行Pauli目標(biāo)極化分解,分析了研究區(qū)不同地物的散射特性,最后利用區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)進(jìn)行了潮溝信息提取,得到主要結(jié)果:

1)利用目標(biāo)分解技術(shù)對(duì)長(zhǎng)江口九段沙進(jìn)行Pauli分解,可以將九段沙不同地物類型的散射回波簡(jiǎn)單化。

2)Pauli目標(biāo)極化分解出的奇次散射、二面角散射和45°二面角散射的特征參數(shù)能夠反映出水體、潮溝、低矮植被和蘆葦之間的地物差異,其中二面角散射差異最大。潮溝的二面角散射強(qiáng)度值為0.007,低矮植被和蘆葦?shù)纳⑸鋸?qiáng)度值為0.783和0.104。

3)利用二面角散射參數(shù),結(jié)合區(qū)域生長(zhǎng)法和數(shù)學(xué)形態(tài)學(xué)可以有效地提取潮溝信息。