孟令奎，毛旭東，魏祖帥，張 文

1. 武漢大學遙感信息工程學院，湖北 武漢 430079； 2. 地球空間信息技術(shù)協(xié)同創(chuàng)新中心，湖北 武漢 430079

用遙感影像提取水體信息是水利行業(yè)解決洪澇災(zāi)害監(jiān)測、水資源管理等需求的重要方式。SAR數(shù)據(jù)以其全天候、全天時的成像特點，適用于洪澇災(zāi)害信息的提取和監(jiān)測[1]，在水利遙感領(lǐng)域發(fā)揮著重要作用。作為我國自主研發(fā)的首顆分辨率高達1 m的C波段SAR衛(wèi)星[2]，GF-3的投入運行為水利遙感帶來了更為廣闊的發(fā)展前景。

經(jīng)典的遙感影像水體提取方法包括閾值分割[3-4]、結(jié)合多種特征的分類器分割[5]、基于水平集理論的分割[1,6]、基于能量函數(shù)的圖像分割[7-9]等。當研究區(qū)域環(huán)境復(fù)雜時，水體提取過程常常出現(xiàn)虛警、漏警現(xiàn)象，導(dǎo)致水體提取結(jié)果精度不高。這種環(huán)境復(fù)雜性帶來的精度問題，一定程度上源自像元分類的不確定性[10]。根據(jù)某些屬性對像元分類時，可能存在將某些與水體像元散射特征或是光譜特征相似，但實際上不是水體的像元錯分為水體的現(xiàn)象，使得像元分類結(jié)果存在不確定性；當數(shù)據(jù)的分辨率不高時，影像中某些像元反映的實際區(qū)域可能同時包含了多種地物，在分類上存在著不確定性。此外，在灰度值分割中閾值的選取同樣會造成分類結(jié)果的不確定性。顯然，二值水體分布圖式的水體提取結(jié)果無法反映上述不確定性。

不同于傳統(tǒng)的圖像分割思路，水體分布概率估計[11]基于貝葉斯推斷，根據(jù)SAR影像后向散射系數(shù)值計算水體的概率分布圖。其最大的特點在于結(jié)果圖上每一點取值范圍不再是離散的{0,1}集合，而是連續(xù)的[0,1]區(qū)間，該值表示該點屬于水體的概率。與常規(guī)水體提取方法得到二值水體分布圖形式的結(jié)果不同，水體分布概率估計方法得到的是信息量更豐富的概率圖，能夠定量反映SAR影像中水體像元分類存在的不確定性。水體分布概率圖充分描述了影像中水體像元的分布信息，可以用索引圖的形式將水體概率值映射為不同顏色、進而直觀表示水體像元的分布，也可以通過閾值分割退化為二值水體提取結(jié)果，還可以應(yīng)用于水位估計[12-13]、水動力學模型校準[14-16]等場景。

水體分布先驗概率是水體概率估計方法中的一個重要參數(shù)。在先驗概率估計方面，現(xiàn)有的概率估計方法[11,17]將先驗概率默認置為0.5，這會使推導(dǎo)出的地物后向散射系數(shù)統(tǒng)計模型不夠準確，降低概率估計結(jié)果的精度。針對先驗概率估計不充分的問題，本文提出了改進的水體概率估計方法，借助k-means聚類分析分割影像的結(jié)果估計先驗概率，并利用這一先驗概率估計值優(yōu)化了后向散射系數(shù)分布參數(shù)估計方法。本文試驗以GF-3影像為數(shù)據(jù)源，結(jié)果表明改進的水體概率估計方法能夠有效提升水體概率估計精度。

1 水體概率估計理論

(1)

(2)

對于影像中的某個像元，希望根據(jù)其后向散射系數(shù)σ0判斷其類別，即計算p(W|σ0)。根據(jù)條件概率公式，有

pW|σ0p(σ0)=p(W)pσ0|W

(3)

則p(W|σ0)可以寫作

pW|σ0=Ppriorpσ0|W/p(σ0)

(4)

(5)

(6)

相應(yīng)的，結(jié)合式(2)，邊際分布p(σ0)表現(xiàn)為含有兩個簇的混合Gauss分布形式。理論上其概率密度函數(shù)恰好具有雙峰曲線的形式，這與實際分析得到的直方圖形態(tài)特征一致。

2 方法設(shè)計

2.1 先驗概率估計

關(guān)于水體分布先驗概率估計，在文獻[11]的方法中，先驗概率值默認設(shè)置為0.5，稱為默認先驗方法。其假設(shè)的情況是研究區(qū)域中水體、背景像元數(shù)相同的情況，然而該條件在大多數(shù)情況下遠不成立。注意到含有水體的影像后向散射系數(shù)邊際分布p(σ0)的混合Gauss分布模型形式，考慮采用k-means算法對各像元進行聚類分析，并用聚類的結(jié)果估計水體分布先驗概率。

k-means算法以平方誤差為標準，通過迭代對給定樣本集尋找簇內(nèi)樣本相似度最高的簇劃分，具體算法步驟不再贅述。根據(jù)p(σ0)的混合Gauss分布模型形式，結(jié)合統(tǒng)計分布直方圖h(σ0)呈現(xiàn)出的雙峰形態(tài)，研究區(qū)域所有像元的σ0組成的樣本集合按取值大小聚集為兩簇。并且根據(jù)水體像元后向散射系數(shù)普遍低于背景像元這一特點，可以直接判定聚類得到的低σ0簇為水體像元簇，高σ0簇為背景像元簇。得到聚類結(jié)果后，計算水體像元比例作為水體分布先驗概率Pprior的估計值。此外，統(tǒng)計兩簇樣本的均值和標準差用于初始化后續(xù)參數(shù)估計中的迭代步驟。

2.2 σ0分布參數(shù)估計

(1) 計算研究區(qū)域影像后向散射系數(shù)統(tǒng)計分布直方圖。統(tǒng)計分布直方圖h(σ0)計算需要的參數(shù)為其每個條帶的寬度，稱為帶寬。本次試驗中帶寬的計算公式為

(7)

式中，n為統(tǒng)計樣本個數(shù)；IQR為樣本的四分位差。研究表明[21]，由上式給出的帶寬適用于Gauss分布樣本數(shù)據(jù)。

(2) 疊加理論概率密度函數(shù)p(σ0)。直方圖h(σ0)的取值是像元個數(shù)，概率密度函數(shù)p(σ0)取值是概率值，二者在曲線擬合之前需經(jīng)過疊加[21]操作。具體方法為先對p(σ0)乘以系數(shù)A，再用Ap(σ0)對h(σ0)作曲線擬合，其中A表示直方圖的面積(帶寬×高度)。

(3) 非線性擬合。使用Levenberg-Marquardt算法進行非線性最小二乘擬合[11,22]，迭代計算分布參數(shù)。迭代初始化時，用聚類步驟中得到的水體、背景樣本均值和標準差作為算法初值。

2.3 方法流程

綜上所述，本文方法的主要流程如圖1所示。由于引入了先驗概率估計的步驟，理論上本文方法有兩點優(yōu)勢：

(1) 精度更高。先用聚類估計先驗概率、再用非線性擬合估計分布參數(shù)的流程會估計得到的使p(σ0)更貼合h(σ0)，即更能貼合σ0分布的實際統(tǒng)計規(guī)律，理論上會獲得更高的概率估計精度。

(2) 效率更高。在估計模型參數(shù)初始化值設(shè)置時，默認先驗的方法采用人工設(shè)置經(jīng)驗初值，當初值選取不合適時可能會出現(xiàn)求解過程中迭代次數(shù)過多、甚至找不到局部最優(yōu)解的情況。本文方法用聚類結(jié)果數(shù)據(jù)初始化模型參數(shù)，效率更高，同時也避免了因初值設(shè)置不當而可能出現(xiàn)的找不到最優(yōu)解的情況。

圖1 本文方法流程Fig.1 Workflow of the algorithm

3 試驗結(jié)果

3.1 試驗數(shù)據(jù)

試驗使用了GF-3 L2級產(chǎn)品，包含F(xiàn)SI、FSII、QPSI等3種成像模式，分別對應(yīng)3、10與8 m分辨率。驗證數(shù)據(jù)由研究區(qū)域的GF-1影像人工標注得來，影像的拍攝時間與試驗用GF-3影像相距在15 d以內(nèi)。在開始數(shù)據(jù)分析前對影像應(yīng)用了窗口大小為5×5的Gamma濾波，研究表明[23-24]這一選擇可以在盡可能地保持影像空間細節(jié)的情況下提高信噪比。

3.2 結(jié)果可靠性評價

選取研究區(qū)域的GF-1影像，采用人工標注的方法，得到二值水體區(qū)域圖作為真值，采用可靠性圖(reliability diagram)[11,25]的方式驗證結(jié)果。將概率取值的[0,1]區(qū)間均勻分割為N=10個小區(qū)間Ik，即

(8)

(9)

3.3 結(jié)合先驗概率估計的水體概率估計方法試驗結(jié)果

以河北省鹿泉市黃壁莊水庫周邊區(qū)域和湖北省黃岡市龍感湖周邊區(qū)域為例開展試驗。這些區(qū)域內(nèi)水體集中分布于水庫庫區(qū)或湖區(qū)及周邊河道內(nèi)，其他部分則均為典型的背景地物，成分簡單且分布互不混淆，適用于初步測試本文方法的效果，試驗結(jié)果如圖2、圖3和表1。圖2(a)為河北鹿泉地區(qū)經(jīng)輻射定標得到的GF-3后向散射系數(shù)影像?？梢钥闯雠c非水地物相比，水體因其低σ0值而呈現(xiàn)為明顯的暗色區(qū)域。圖2(b)為用作參考的GF-1影像，分辨率為8 m。圖2(c)為人工標注的水體區(qū)域真值圖。圖2(d)為水體概率估計的結(jié)果，由顏色的深淺變化對應(yīng)條件概率p(W|σ0)不同大小。結(jié)合圖2(c)和圖2(d)來看，概率估計有效提取出了水體區(qū)域。水庫庫區(qū)及其周邊河道的水體條件概率取值很高(pW|σ0>0.9)，與周圍地物區(qū)分明顯。湖北黃岡地區(qū)試驗結(jié)果與鹿泉地區(qū)情況相似，如圖3所示。

表1 參數(shù)估計結(jié)果

圖2 鹿泉地區(qū)試驗結(jié)果Fig.2 The experimental result in Luquan area

圖3 黃岡地區(qū)試驗結(jié)果Fig.3 The experimental result in Huanggang area

圖4 水體概率分析圖Fig.4 Analysis of probabilistic water body map

3.4 與默認先驗方法的結(jié)果精度對比

對比試驗研究區(qū)域地處湖北咸寧地區(qū)長江沿岸。該區(qū)域包含水田、溝渠、河道等零散細小水體以及大范圍的灘涂、濕地等難以與水體區(qū)分的地物，可用于檢驗本文方法在復(fù)雜地物情況下的效果，數(shù)據(jù)及結(jié)果見圖5。同時，使用默認先驗概率估計方法用作對照，對比分析的結(jié)果見圖6。

圖5 咸寧地區(qū)試驗結(jié)果Fig.5 The experimental result in Xianning area

圖6 水體概率圖分析Fig.6 Analysis of probabilistic water body map

3.5 先驗概率對概率估計精度的影響

前文的試驗表明結(jié)合先驗概率估計的水體概率估計較默認先驗的估計方法在精度上的提升，并且推測這一結(jié)果差異與先驗概率估計準確與否有關(guān)。對于先驗概率對結(jié)果帶來的影響，可以作進一步的解釋。

從之前的研究區(qū)域中截取7個大小為200×200像元的子區(qū)域。這些子區(qū)域中水體面積所占區(qū)域總面積的比例各不相同，即各區(qū)域的水體分布先驗概率各不相同。而在這些區(qū)域?qū)嵤┧w概率估計時，默認先驗的方法統(tǒng)一將先驗概率估計值置為0.5，這就造成了估計值與實際值的偏差，且各子區(qū)域偏差值各不相同。根據(jù)7個子區(qū)域的試驗結(jié)果，嘗試分析先驗概率對結(jié)果帶來的影響，試驗結(jié)果見表2及圖7。

表2 各子區(qū)域的試驗結(jié)果

圖7 可靠性度量比較分析Fig.7 Comparison of reliability

默認先驗方法的Re值按1到7的順序整體呈現(xiàn)先減小再增大的“U”形走勢。這與默認先驗概率值與各區(qū)域先驗概率實際值偏差大小的變化規(guī)律一致。說明Re與先驗概率偏差呈正相關(guān)，意味著準確的先驗概率可以提升水體概率估計精度。

根據(jù)可靠性度量的減小比例，除去第5個區(qū)域的結(jié)果，在其他6個區(qū)域上本文方法精度均優(yōu)于默認先驗的方法，且精度提升比例按1到7的順序也呈兩端大中間小的“U”形走勢。這一現(xiàn)象是合理的，根據(jù)表第2、4行，本文方法對先驗概率的估計較為準確，相應(yīng)的先驗概率偏差值保持在較小的范圍內(nèi)，而默認先驗方法的先驗概率偏差值則會出現(xiàn)“U”形波動。這一結(jié)果也再次驗證了水體概率估計精度與先驗概率偏差呈正相關(guān)關(guān)系，結(jié)合先驗概率估計的水體概率估計方法精度更高。

3.6 從概率圖到二值圖

即97.51%的像元被正確分類。剩余2.49%分類錯誤像元在各概率區(qū)間的分布情況如表3所示。結(jié)合圖9，分類錯誤像元在[0,0.1]與[0.9,1]區(qū)間占比較高，但這兩個區(qū)間內(nèi)像元所占像元總數(shù)比例也高，因此分類錯誤像元占區(qū)間像元總數(shù)的比例反而較低。與之相反的是，8個中間概率區(qū)間只含有影像10.3%的像元，分類錯誤部分卻占了錯誤分類像元總數(shù)的57.5%，占各區(qū)間像元數(shù)目的比例也遠遠高于[0,0.1]和[0.9,1]兩個區(qū)間。這說明具有中間概率值的像元更容易被錯誤分類。特別取出這部分像元(水體概率取值于0.1到0.9)，如圖8(b)，發(fā)現(xiàn)普遍位于水體邊界地帶。該類區(qū)域多為濕地、沼澤，常覆蓋有淺層水面，呈現(xiàn)出與水體類似的低后向散射特征，易于出現(xiàn)虛警、漏警現(xiàn)象，僅依據(jù)后向散射系數(shù)無法準確分類此種地物。這是閾值劃分方法不可避免的缺陷，取0.5作為水體概率的分割閾值恰好是一種折中的思路。

圖8 南昌地區(qū)試驗結(jié)果Fig.8 The experimental result in Nanchang area

對水體概率圖做閾值劃分本質(zhì)上是對后向散射系數(shù)影像做閾值劃分。本節(jié)所述的由水體概率圖退化得到水體二值圖的方法屬于閾值法提取水體的范疇。選擇水體概率取0.5時的后向散射系數(shù)值作為分割閾值，對應(yīng)理論分布概率密度曲線上水體與背景兩類像元兩組密度函數(shù)曲線的交點處，也即直方圖兩峰值之間的“波谷”處，與一些閾值劃分方法中通過尋找直方圖“波谷”確定閾值的思路一致。不過與之不同的是，本節(jié)所述水體提取方法根據(jù)兩類像元理論分布推導(dǎo)而來，在像元分布估計準確的前提下，閾值選取的效果必然是好的。與直接在直方圖上尋找“波谷”的方法相比，該方法具有更強的理論基礎(chǔ)，同時會減少直方圖中噪聲帶來的誤差。

表3 各概率區(qū)間像元分布

注：n為區(qū)間像元數(shù)；N為像元總數(shù)；nerror為區(qū)間分類錯誤像元數(shù)；Nerror為分類錯誤像元總數(shù)。

圖9 各概率區(qū)間像元分布Fig.9 The distribution of pixels in each interval

4 結(jié) 語

本文針對水體概率估計中水體分布先驗概率估計不充分的問題，使用k-means聚類估計先驗概率并優(yōu)化了后續(xù)的參數(shù)估計步驟，改進了水體概率估計方法。在水體概率估計方法中，先驗概率估計的偏差會帶來模型對地物后向散射系數(shù)分布規(guī)律描述的偏差，進而影響水體概率估計精度。本文方法由于對先驗概率估計更充分，進而結(jié)果精度更高。此外，本文方法使用聚類結(jié)果數(shù)據(jù)初始化模型參數(shù)讓參數(shù)迭代求解過程能夠更快地收斂至最優(yōu)值，同時也簡化了人工迭代初值設(shè)置的步驟。此外，由于該方法使用的是單極化數(shù)據(jù)，因此同樣適用于多時相或是多種極化方式的影像，具有很好的推廣性。

需要指出的是，本文方法的推導(dǎo)建立在研究區(qū)域后向散射系數(shù)統(tǒng)計分布直方圖呈雙峰形這一前提條件上，根據(jù)水域與背景間后向散射系數(shù)的差異提出模型假設(shè)。但是這一條件未必總是成立，例如當影像中水體與背景像元中的一類數(shù)目遠遠大于另一類時，其后向散射系數(shù)分布會被一類像元主導(dǎo)；或是風浪導(dǎo)致的水體后向散射系數(shù)變高并與背景像元混淆，直方圖也無法呈現(xiàn)雙峰形走勢，影響概率估計效果。對于此類問題及本文方法適用條件的定量描述，是今后研究工作的重點。在實際應(yīng)用中，可以通過設(shè)置河流、湖泊等水體目標緩沖區(qū)等方式避免研究區(qū)域水體比例出現(xiàn)極端值的情況，保證概率估計的效果。