結合先驗概率估計的GF-3影像水體概率估計方法
孟令奎1,2, 毛旭東1, 魏祖帥1, 張文1
1. 武漢大學遙感信息工程學院, 湖北 武漢 430079;
2. 地球空間信息技術協同創新中心, 湖北 武漢 430079
收稿日期：2017-12-20；修回日期：2018-08-24
基金項目：國家重點研發計劃（2017YFC0405806）
第一作者簡介：孟令奎(1967-), 男, 教授, 研究方向為網絡GIS、水利遙感技術及應用。E-mail:lkmeng@whu.edu.cn
通信作者：張文, E-mail：wen_zhang@whu.edu.cn
摘要：綜合SAR（synthetic aperture radar）影像的統計模型假設與k-means聚類算法，提出了一種結合水體分布先驗概率估計的水體概率估計方法。首先，用貝葉斯推斷對研究區域後向散射係數做統計模型假設。隨後，結合聚類算法對像元作分類，估計水體分布先驗概率，結合統計分布直方圖使用非線性最小二乘擬合完成模型參數估計。試驗選取了高分三號（GF-3）多種工作模式數據，並用高分一號（GF-1）影像進行驗證。結果表明，該方法可有效實現SAR影像的高精度水體概率估計。
關鍵詞：水體概率估計 參數估計 合成孔徑雷達(SAR) 高分三號
Probabilistic water body mapping of GF-3 images based on prior probability estimation
MENG Lingkui1,2, MAO Xudong1, WEI Zushuai1, ZHANG Wen1
1. School of Remote Sensing and Information Engineering, Wuhan University, Wuhan 430079, China;
2. Collaborative Innovation Center of Geospatial Technology, Wuhan 430079, China
Foundation support: The National Key Research and Development Program of China (No. 2017YFC0405806)
First author: MENG Lingkui(1967—), male, professor, majors in web GIS, water resources remote sensing technology and application.E-mail:lkmeng@whu.edu.cn.
Corresponding author: ZHANG Wen, E-mail: wen_zhang@whu.edu.cn.
Abstract: We combine k-means cluster algorithm with the statistical model of SAR(synthetic aperture radar) images and develop the probabilistic water body mapping algorithm based on the priori probability estimation. Firstly, we make the statistical model assumption about backscatter values based on Bayesian theory. Then, we classify the images based on cluster algorithm, calculate the prior probability of the water body mapping and estimate the parameters of the statistical model of water distribution.Thewater body probabilistic maps based on GF-3 images in Luquan and Xianning are calculated and then validated with GF-1 images. The algorithm is effective on high-precision probabilistic water body mapping of SAR images.
Key words: probabilistic water body mapping parameter estimation synthetic aperture radar (SAR) GF-3
用遙感影像提取水體信息是水利行業解決洪澇災害監測、水資源管理等需求的重要方式。SAR數據以其全天候、全天時的成像特點，適用於洪澇災害信息的提取和監測[1]，在水利遙感領域發揮著重要作用。作為我國自主研發的首顆解析度高達1 m的C波段SAR衛星[2]，GF-3的投入運行為水利遙感帶來了更為廣闊的發展前景。
經典的遙感影像水體提取方法包括閾值分割[3-4]、結合多種特徵的分類器分割[5]、基於水平集理論的分割[1, 6]、基於能量函數的圖像分割[7-9]等。當研究區域環境複雜時，水體提取過程常常出現虛警、漏警現象，導致水體提取結果精度不高。這種環境複雜性帶來的精度問題，一定程度上源自像元分類的不確定性[10]。根據某些屬性對像元分類時，可能存在將某些與水體像元散射特徵或是光譜特徵相似，但實際上不是水體的像元錯分為水體的現象，使得像元分類結果存在不確定性；當數據的解析度不高時，影像中某些像元反映的實際區域可能同時包含了多種地物，在分類上存在著不確定性。此外，在灰度值分割中閾值的選取同樣會造成分類結果的不確定性。顯然，二值水體分布圖式的水體提取結果無法反映上述不確定性。
不同於傳統的圖像分割思路，水體分布概率估計[11]基於貝葉斯推斷，根據SAR影像後向散射係數值計算水體的概率分布圖。其最大的特點在於結果圖上每一點取值範圍不再是離散的{0, 1}集合，而是連續的[0, 1]區間，該值表示該點屬於水體的概率。與常規水體提取方法得到二值水體分布圖形式的結果不同，水體分布概率估計方法得到的是信息量更豐富的概率圖，能夠定量反映SAR影像中水體像元分類存在的不確定性。水體分布概率圖充分描述了影像中水體像元的分布信息，可以用索引圖的形式將水體概率值映射為不同顏色、進而直觀表示水體像元的分布，也可以通過閾值分割退化為二值水體提取結果，還可以應用於水位估計[12-13]、水動力學模型校準[14-16]等場景。
水體分布先驗概率是水體概率估計方法中的一個重要參數。在先驗概率估計方面，現有的概率估計方法[11, 17]將先驗概率默認置為0.5，這會使推導出的地物後向散射係數統計模型不夠準確，降低概率估計結果的精度。針對先驗概率估計不充分的問題，本文提出了改進的水體概率估計方法，藉助k-means聚類分析分割影像的結果估計先驗概率，並利用這一先驗概率估計值優化了後向散射係數分布參數估計方法。本文試驗以GF-3影像為數據源，結果表明改進的水體概率估計方法能夠有效提升水體概率估計精度。
1 水體概率估計理論
水體概率估計理論依賴於像元後向散射係數統計分布假設與貝葉斯推斷。將研究區域所有像元組成的樣本集合視作水體(W)和背景(W)兩部分的非交併集[18]，並用後向散射係數(σ0)描述其特徵。記水體像元的後向散射係數分布概率密度函數為p(σ0|W)，背景像元的後向散射係數分布概率密度函數為p(σ0|W)，研究區域後向散射係數邊際分布的概率密度函數為p(σ0)，則
(1)
式中，p(W)與p(W)分別表示水體與背景像元所佔比例，p(W)+p(W)=1。即p(σ0)可視為兩類像元的混合分布。
記Pprior=p(W)，p(W)=1-Pprior，式(1)可寫作
(2)
對於影像中的某個像元，希望根據其後向散射係數σ0判斷其類別，即計算p(W|σ0)。根據條件概率公式，有
(3)
則p(W|σ0)可以寫作
(4)
根據式(2)、式(4)，為了得到p(W|σ0)，需要分別估計Pprior、p(σ0|W)、p(σ0|W)。為此，下面對p(σ0|W)、p(σ0|W)作進一步的分布假設。
水體由於其鏡面反射現象，在SAR影像中呈現為低後向散射係數的均勻區域，含有一定範圍水體的SAR影像後向散射係數統計分布直方圖h(σ0)常聚為大小兩類，整體呈雙峰狀。這種山峰型走勢恰好與Gauss分布的概率密度函數一致。事實上，已有的研究發現經過對數變換的SAR強度數據服從Gauss分布[19-20]。根據文獻[11]在洪水概率估計(probabilistic flood mapping)中對水體、背景兩類像元的Gauss分布假設，本文同樣假設p(σ0|W)與p(σ0|W)均服從Gauss分布，即
(5)
(6)
相應的，結合式(2)，邊際分布p(σ0)表現為含有兩個簇的混合Gauss分布形式。理論上其概率密度函數恰好具有雙峰曲線的形式，這與實際分析得到的直方圖形態特徵一致。
確定後向散射係數分布假設後，估計模型中的未知參數即可確定其統計特徵，給定的σ0即可計算其屬於水體的概率值p(W|σ0)。從上述水體概率估計的理論來看，本文方法的核心是參數估計，包括水體分布先驗概率Pprior和後向散射係數分布參數(μW, sW, μW, sW)。本文從先驗概率估計方法上著手改進。考慮到p(σ0)的混合Gauss分布形式，引入k-means聚類算法估計先驗概率，並將其運用於分布參數估計方法中，優化參數估計方法。
2 方法設計2.1 先驗概率估計
關於水體分布先驗概率估計，在文獻[11]的方法中，先驗概率值默認設置為0.5，稱為默認先驗方法。其假設的情況是研究區域中水體、背景像元數相同的情況，然而該條件在大多數情況下遠不成立。注意到含有水體的影像後向散射係數邊際分布p(σ0)的混合Gauss分布模型形式，考慮採用k-means算法對各像元進行聚類分析，並用聚類的結果估計水體分布先驗概率。
k-means算法以平方誤差為標準，通過迭代對給定樣本集尋找簇內樣本相似度最高的簇劃分，具體算法步驟不再贅述。根據p(σ0)的混合Gauss分布模型形式，結合統計分布直方圖h(σ0)呈現出的雙峰形態，研究區域所有像元的σ0組成的樣本集合按取值大小聚集為兩簇。並且根據水體像元後向散射係數普遍低於背景像元這一特點，可以直接判定聚類得到的低σ0簇為水體像元簇，高σ0簇為背景像元簇。得到聚類結果後，計算水體像元比例作為水體分布先驗概率Pprior的估計值。此外，統計兩簇樣本的均值和標準差用於初始化後續參數估計中的迭代步驟。
2.2 σ0分布參數估計
統計學上認為直方圖是對概率密度函數的估計[21]。之前的推導已經得到了Pprior的估計值，接下來直接將理論概率密度函數p(σ0)疊加到實際分析研究區域影像數據得到的直方圖h(σ0)之上進行曲線擬合，估計未知分布參數(μW, sW, μW, sW)。參數估計過程主要為：
(1) 計算研究區域影像後向散射係數統計分布直方圖。統計分布直方圖h(σ0)計算需要的參數為其每個條帶的寬度，稱為帶寬。本次試驗中帶寬的計算公式為
(7)
式中，n為統計樣本個數；IQR為樣本的四分位差。研究表明[21]，由上式給出的帶寬適用於Gauss分布樣本數據。
(2) 疊加理論概率密度函數p(σ0)。直方圖h(σ0)的取值是像元個數，概率密度函數p(σ0)取值是概率值，二者在曲線擬合之前需經過疊加[21]操作。具體方法為先對p(σ0)乘以係數A，再用Ap(σ0)對h(σ0)作曲線擬合，其中A表示直方圖的面積(帶寬×高度)。
(3) 非線性擬合。使用Levenberg-Marquardt算法進行非線性最小二乘擬合[11, 22]，迭代計算分布參數。迭代初始化時，用聚類步驟中得到的水體、背景樣本均值和標準差作為算法初值。
得到分布參數(μW, sW, μW, sW)的估計值後，代入式(4)即得到水體概率估計的結果式。給定某像元的後向散射係數σ0，由式(4)即可得到該像元屬於水體的概率p(W|σ0)。
2.3 方法流程
綜上所述，本文方法的主要流程如圖 1所示。由於引入了先驗概率估計的步驟，理論上本文方法有兩點優勢：
圖 1 本文方法流程Fig. 1 Workflow of the algorithm(1) 精度更高。先用聚類估計先驗概率、再用非線性擬合估計分布參數的流程會估計得到的使p(σ0)更貼合h(σ0)，即更能貼合σ0分布的實際統計規律，理論上會獲得更高的概率估計精度。
(2) 效率更高。在估計模型參數初始化值設置時，默認先驗的方法採用人工設置經驗初值，當初值選取不合適時可能會出現求解過程中迭代次數過多、甚至找不到局部最優解的情況。本文方法用聚類結果數據初始化模型參數，效率更高，同時也避免了因初值設置不當而可能出現的找不到最優解的情況。
3 試驗結果3.1 試驗數據
試驗使用了GF-3 L2級產品，包含FSI、FSII、QPSI等3種成像模式，分別對應3、10與8 m解析度。驗證數據由研究區域的GF-1影像人工標註得來，影像的拍攝時間與試驗用GF-3影像相距在15 d以內。在開始數據分析前對影像應用了窗口大小為5×5的Gamma濾波，研究表明[23-24]這一選擇可以在儘可能地保持影像空間細節的情況下提高信噪比。
3.2 結果可靠性評價
選取研究區域的GF-1影像，採用人工標註的方法，得到二值水體區域圖作為真值，採用可靠性圖(reliability diagram)[11, 25]的方式驗證結果。將概率取值的[0, 1]區間均勻分割為N=10個小區間Ik，即
(8)
水體概率圖上任意一個像元的值則屬於上述10個小區間中的某一個。記取值在Ik的像元為集合Ωk，則Ik中每個像元屬於水體的概率在Ik區間範圍內。那麼理論上Ωk中水體像元的比例在Ik區間範圍內，取其中值記為Pk。同時，以二值水體區域圖為真值參考，統計Ωk中水體像元的實際比例Pk，則可以用理論比例Pk與實際比例Pk的偏差對概率圖的精度做定量評價，以判斷本文方法所計算的水體概率值是否準確。以Pk和Pk(k=1, 2, …, N)為橫縱坐標所作的圖即是可靠性圖。具體的，對於水體概率圖和給定的二值水體區域圖，有如下的可靠性度量[25-27]
(9)
式中，nk為研究區域內取值在Ik上的像元數目。可靠性度量是以加權均方根誤差的形式給出的，表示理論比例Pk與實際比例Pk的偏差程度。Re值越小，水體概率圖精度越高。
3.3 結合先驗概率估計的水體概率估計方法試驗結果
以河北省鹿泉市黃壁莊水庫周邊區域和湖北省黃岡市龍感湖周邊區域為例開展試驗。這些區域內水體集中分布於水庫庫區或湖區及周邊河道內，其他部分則均為典型的背景地物，成分簡單且分布互不混淆，適用於初步測試本文方法的效果，試驗結果如圖 2、圖 3和表 1。圖 2(a)為河北鹿泉地區經輻射定標得到的GF-3後向散射係數影像。可以看出與非水地物相比，水體因其低σ0值而呈現為明顯的暗色區域。圖 2(b)為用作參考的GF-1影像，解析度為8 m。圖 2(c)為人工標註的水體區域真值圖。圖 2(d)為水體概率估計的結果，由顏色的深淺變化對應條件概率p(W|σ0)不同大小。結合圖 2(c)和圖 2(d)來看，概率估計有效提取出了水體區域。水庫庫區及其周邊河道的水體條件概率取值很高(p(W|σ0>0.9)，與周圍地物區分明顯。湖北黃岡地區試驗結果與鹿泉地區情況相似，如圖 3所示。
表 1 參數估計結果Tab. 1 Parameter estimate
研究區域PpriorμWsWμWsW鹿泉0.277 7-30.912 21.736 8-19.331 61.542 7黃岡0.402 3-28.341 61.329 3-17.640 01.425 0圖 4(a)為兩個區域試驗的p(W|σ0)取值在I1到I10各個區間內像元數目的柱狀圖。從圖上看，兩個邊緣區間I1和I10包含了絕大部分像元，表明影像中絕大部分區域都被可靠地區分為了水體或背景，而這一理想結果的產生得益於該地區簡單的地物成分以及互不混淆的地物分布狀態。圖 4(b)為可靠性圖，其中(Pk, Pk點均靠近於直線y=x，表明根據人工標註的驗證數據，本文方法對像元屬於水體的概率值的估計較為準確。
圖 4 水體概率分析圖Fig. 4 Analysis of probabilistic water body map
3.4 與默認先驗方法的結果精度對比
對比試驗研究區域地處湖北鹹寧地區長江沿岸。該區域包含水田、溝渠、河道等零散細小水體以及大範圍的灘涂、溼地等難以與水體區分的地物，可用於檢驗本文方法在複雜地物情況下的效果，數據及結果見圖 5。同時，使用默認先驗概率估計方法用作對照，對比分析的結果見圖 6。
圖 5 鹹寧地區試驗結果Fig. 5 The experimental result in Xianning area圖 6(a)在σ0灰度直方圖上同時疊加了本文方法估計出的σ0理論概率密度函數p1(σ0)與默認先驗方法估計出的σ0理論概率密度函數p2(σ0)。從曲線形狀上看，p1(σ0)與直方圖高度相同、走勢一致，契合度高，說明本文方法對σ0分布參數的估計基本反映出了真實規律；p2(σ0)曲線雖然也呈現出與直方圖一致的雙峰走勢，但是對一個波峰峰值估計過高、對第二個波峰峰值估計過低，造成曲線整體與直方圖的契合度偏低，未能刻畫出σ0的分布規律。從參數估計結果上看，兩種方法的估計值幾乎是相同的，唯一的區別在於先驗概率。根據式(2)，當兩個分布的參數(μW, μW, sW, sW)一致時，先驗概率Pprior作為比例係數決定了理論σ0分布概率密度函數曲線p(σ0)的形狀。因此當先驗概率偏差較大時，其對概率估計結果的影響直觀上就體現為直方圖h(σ0)與疊加其上的p(σ0)曲線不匹配，本質上則是對研究區域地物後向散射係數σ0分布規律描述不夠準確的表現。圖 6(b)可靠性圖從定量分析的角度描述了兩種概率估計方法精度上的差異。從圖上看，本文方法的Pk值在k=1, 2, …, 10處均高於默認先驗方法得到的結果、對應的點(Pk, Pk)均更靠近直線y=x，這反映出不論σ0取值高低，即在各種地物狀態條件下，本文方法對像元屬於水體的概率值的估計均取得了更好的效果。採用可靠性度量Re的評價方式，本文方法較默認先驗方法結果誤差降低了17.6%。
3.5 先驗概率對概率估計精度的影響
前文的試驗表明結合先驗概率估計的水體概率估計較默認先驗的估計方法在精度上的提升，並且推測這一結果差異與先驗概率估計準確與否有關。對於先驗概率對結果帶來的影響，可以作進一步的解釋。
從之前的研究區域中截取7個大小為200×200像元的子區域。這些子區域中水體面積所佔區域總面積的比例各不相同，即各區域的水體分布先驗概率各不相同。而在這些區域實施水體概率估計時，默認先驗的方法統一將先驗概率估計值置為0.5，這就造成了估計值與實際值的偏差，且各子區域偏差值各不相同。根據7個子區域的試驗結果，嘗試分析先驗概率對結果帶來的影響，試驗結果見表 2及圖 7。
表 2 各子區域的試驗結果Tab. 2 Results of sub-regions
研究區域1234567先驗概率實際值0.1530.2540.3420.5050.6020.6960.832先驗概率估計值(默認先驗)0.5000.5000.5000.5000.5000.5000.500先驗概率估計值(本文方法)0.1630.3230.3910.5570.6810.7110.831默認先驗Re10.201 30.102 10.105 90.075 50.070 00.071 60.096 3本文方法Re20.142 10.089 50.089 90.075 50.076 00.064 10.065 8(Re1-Re2)/Re129.4%12.4%15.1%0-8.6%10.6%31.7%默認先驗方法的Re值按1到7的順序整體呈現先減小再增大的「U」形走勢。這與默認先驗概率值與各區域先驗概率實際值偏差大小的變化規律一致。說明Re與先驗概率偏差呈正相關，意味著準確的先驗概率可以提升水體概率估計精度。
根據可靠性度量的減小比例，除去第5個區域的結果，在其他6個區域上本文方法精度均優於默認先驗的方法，且精度提升比例按1到7的順序也呈兩端大中間小的「U」形走勢。這一現象是合理的，根據表第2、4行，本文方法對先驗概率的估計較為準確，相應的先驗概率偏差值保持在較小的範圍內，而默認先驗方法的先驗概率偏差值則會出現「U」形波動。這一結果也再次驗證了水體概率估計精度與先驗概率偏差呈正相關關係，結合先驗概率估計的水體概率估計方法精度更高。
3.6 從概率圖到二值圖
水體概率圖可以通過閾值分割的方法退化為二值化水體分布圖。考慮如下的劃分準則：若像元屬於水體的概率不小於屬於背景的概率，即p(W|σ0)≥p(W|σ0)時，認為該像元為水體像元。因為p(W|σ0)+p(W|σ0)=1，上述準則即選取p(W|σ0)=0.5為決策面對水體概率圖進行閾值劃分。以一景江西南昌附近的鄱陽湖區域GF-3影像為例，其水體提取結果如圖 8(a)。從圖上看，水體提取結果與水體範圍真值基本一致。其像元分類結果的混淆矩陣是
圖 8 南昌地區試驗結果Fig. 8 The experimental result in Nanchang area即97.51%的像元被正確分類。剩餘2.49%分類錯誤像元在各概率區間的分布情況如表 3所示。結合圖 9，分類錯誤像元在[0, 0.1]與[0.9, 1]區間佔比較高，但這兩個區間內像元所佔像元總數比例也高，因此分類錯誤像元佔區間像元總數的比例反而較低。與之相反的是，8個中間概率區間只含有影像10.3%的像元，分類錯誤部分卻佔了錯誤分類像元總數的57.5%，佔各區間像元數目的比例也遠遠高於[0, 0.1]和[0.9, 1]兩個區間。這說明具有中間概率值的像元更容易被錯誤分類。特別取出這部分像元(水體概率取值於0.1到0.9)，如圖 8(b)，發現普遍位於水體邊界地帶。該類區域多為溼地、沼澤，常覆蓋有淺層水面，呈現出與水體類似的低後向散射特徵，易於出現虛警、漏警現象，僅依據後向散射係數無法準確分類此種地物。這是閾值劃分方法不可避免的缺陷，取0.5作為水體概率的分割閾值恰好是一種折中的思路。
表 3 各概率區間像元分布Tab. 3 The distribution of pixels in each interval
概率區間n/Nnerror/Nerrornerror/N0.0~0.10.5590.4060.0130.1~0.20.0190.1010.0980.2~0.30.0120.0810.1250.3~0.40.0090.0760.1470.4~0.50.0080.0770.1630.5~0.60.0090.0140.0300.6~0.70.0100.0180.0320.7~0.80.0130.0220.0310.8~0.90.0230.0360.0290.9~1.00.3380.1690.009註：n為區間像元數；N為像元總數；nerror為區間分類錯誤像元數；Nerror為分類錯誤像元總數。對水體概率圖做閾值劃分本質上是對後向散射係數影像做閾值劃分。本節所述的由水體概率圖退化得到水體二值圖的方法屬於閾值法提取水體的範疇。選擇水體概率取0.5時的後向散射係數值作為分割閾值，對應理論分布概率密度曲線上水體與背景兩類像元兩組密度函數曲線的交點處，也即直方圖兩峰值之間的「波谷」處，與一些閾值劃分方法中通過尋找直方圖「波谷」確定閾值的思路一致。不過與之不同的是，本節所述水體提取方法根據兩類像元理論分布推導而來，在像元分布估計準確的前提下，閾值選取的效果必然是好的。與直接在直方圖上尋找「波谷」的方法相比，該方法具有更強的理論基礎，同時會減少直方圖中噪聲帶來的誤差。
4 結語
本文針對水體概率估計中水體分布先驗概率估計不充分的問題，使用k-means聚類估計先驗概率並優化了後續的參數估計步驟，改進了水體概率估計方法。在水體概率估計方法中，先驗概率估計的偏差會帶來模型對地物後向散射係數分布規律描述的偏差，進而影響水體概率估計精度。本文方法由於對先驗概率估計更充分，進而結果精度更高。此外，本文方法使用聚類結果數據初始化模型參數讓參數迭代求解過程能夠更快地收斂至最優值，同時也簡化了人工迭代初值設置的步驟。此外，由於該方法使用的是單極化數據，因此同樣適用於多時相或是多種極化方式的影像，具有很好的推廣性。
需要指出的是，本文方法的推導建立在研究區域後向散射係數統計分布直方圖呈雙峰形這一前提條件上，根據水域與背景間後向散射係數的差異提出模型假設。但是這一條件未必總是成立，例如當影像中水體與背景像元中的一類數目遠遠大於另一類時，其後向散射係數分布會被一類像元主導；或是風浪導致的水體後向散射係數變高並與背景像元混淆，直方圖也無法呈現雙峰形走勢，影響概率估計效果。對於此類問題及本文方法適用條件的定量描述，是今後研究工作的重點。在實際應用中，可以通過設置河流、湖泊等水體目標緩衝區等方式避免研究區域水體比例出現極端值的情況，保證概率估計的效果。
【引文格式】孟令奎, 毛旭東, 魏祖帥, 等. 結合先驗概率估計的GF-3影像水體概率估計方法. 測繪學報，2019，48(4)：439-447. DOI: 10.11947/j.AGCS.2019.20170732
權威 | 專業 | 學術 | 前沿
微信投稿郵箱 | song_qi_fan@163.com
微信公眾號中搜索「測繪學報」，關注我們，長按上圖二維碼，關注學術前沿動態。
歡迎加入《測繪學報》作者QQ群： 297834524
進群請備註：姓名+單位+稿件編號