崔 倩，毛旭東，陳德清，瞿 孟，白玨瑩

（1. 水利部信息中心，北京100053；2. 武漢大學(xué)，湖北武漢430072）

0 引言

我國洪澇和干旱災(zāi)害和水資源短缺等水問題突出，地面監(jiān)測系統(tǒng)目前還不能實(shí)現(xiàn)全面覆蓋。高分三號衛(wèi)星具有全天候全天時(shí)對地觀測能力，可以不分晝夜、陰雨天氣進(jìn)行對地觀測。針對洪澇災(zāi)害監(jiān)測、水資源管理等水利應(yīng)用需求，基于可以常態(tài)化接收的國產(chǎn)雷達(dá)數(shù)據(jù)，實(shí)現(xiàn)水體的自動(dòng)化提取，可以補(bǔ)充地面監(jiān)測能力的不足，在提升我國水利監(jiān)測能力方面發(fā)揮了重要作用。

常用的雷達(dá)水體提取技術(shù)包括目視解譯［1］、監(jiān)督分類［2-4］、紋理特征分析［4-5］、動(dòng)態(tài)輪廓模型［6-7］、閾值分割等。監(jiān)督分類法通過構(gòu)造灰度、紋理等特征，標(biāo)注樣本訓(xùn)練支持向量機(jī)等監(jiān)督分類器提取目標(biāo)地物，算法穩(wěn)定有效，但需要人工標(biāo)注樣本，不適用于大批量水體提取工作。紋理特征分析方法從地物粗糙度的差異入手，利用灰度共生矩陣等紋理特征提取地物，提取精度高，但計(jì)算量大，并且效果對紋理窗口尺寸依賴較大。動(dòng)態(tài)輪廓模型將分割圖像提取水體的問題轉(zhuǎn)換為求解能量泛函最小值問題，對噪聲有很強(qiáng)的魯棒性，但對參數(shù)設(shè)置敏感且效率欠佳。閾值分割法是應(yīng)用最為廣泛的水體提取方法，其原理設(shè)定閾值分割雷達(dá)影像，將圖像中小于閾值部分標(biāo)記為水體，大于閾值部分標(biāo)記為背景，獲得水體范圍二值圖。該方法易于實(shí)現(xiàn)，速度快，但難點(diǎn)在于閾值的確定。最直接的方式通過人機(jī)交互的方式選取閾值，但當(dāng)處理大批量數(shù)據(jù)時(shí)效率低下，并且獲取的結(jié)果主觀性強(qiáng)，可追蹤性差。

針對人機(jī)交互閾值選取方法上的不便，自動(dòng)化閾值選取方法被應(yīng)用到水體提取算法中。常用的自動(dòng)閾值分割選取方法有雙峰法、OTSU方法［8-9］、最大熵閾值法［10］和最小誤差閾值法。

出于自動(dòng)化、業(yè)務(wù)化的考慮，所應(yīng)用的雷達(dá)水體提取方法應(yīng)具有簡單、穩(wěn)定、快速的特點(diǎn)，并具有滿足業(yè)務(wù)需求的提取精度。結(jié)合上述常用水體提取方法，文章采用雙峰法進(jìn)行業(yè)務(wù)化水體提取工作。直觀上說，雙峰法尋找閾值的過程，是根據(jù)雷達(dá)影像后向散射系數(shù)統(tǒng)計(jì)分布直方圖呈現(xiàn)出的雙峰值特征、在兩峰之間尋找波谷的過程。

值得注意的是，閾值自動(dòng)選取方法理論上的可行性依賴水體與背景地物的反射特征差異，即雷達(dá)影像中水體像元與背景像元后向散射系數(shù)大小的顯著差異，反應(yīng)在統(tǒng)計(jì)分布直方圖上即圖形的雙峰形走勢。然而，當(dāng)影像中水體地物所占比例很小時(shí)，在直方圖中水體像元的特征被背景像元掩蓋、圖形不再呈現(xiàn)出雙峰形，水體提取效果會受到很大影響。為了避免這一情況的出現(xiàn)，先使用全國水利普查成果湖庫水體緩沖區(qū)矢量對影像做裁剪，再在裁剪結(jié)果上用雙峰法提取水體，保證業(yè)務(wù)化雷達(dá)水體提取方法的精度。

1 雷達(dá)水體自動(dòng)化提取算法

1.1 后向散射系數(shù)分布模型

與可見光和近紅外主要對地物的化學(xué)和熱性質(zhì)敏感不同，微波遙感主要對地物的表面、亞表面及覆蓋物的物理性質(zhì)和介電性質(zhì)敏感。含水量和粗糙度是決定雷達(dá)后向散射的重要因素，一般含水量越大，粗糙度越大，后向散射越強(qiáng)。另外，不同頻率、入射角和不同極化方式的后向散射也有明顯不同。

對于水體來說，未擾動(dòng)水面的后向散射系數(shù)非常低，在所有雷達(dá)波段都接近于鏡面反射，微波能量對于高介電常數(shù)幾乎以等于入射角的方向反射出去。就粗糙度看，平靜水面相對于雷達(dá)波來說幾乎是水平的，因此水的雷達(dá)信號非常弱，區(qū)別于其他地物，在影像中呈現(xiàn)為低后向散射系數(shù)的均勻區(qū)域。在水體提取任務(wù)中，通常只對某個(gè)像元是否屬于水體感興趣，對于非水體的背景像元的具體地物類型并不在意?？紤]到水體像元的低后向散射特性，將像元?jiǎng)澐譃樗w與背景2類，并進(jìn)一步討論這2類像元的特征。

用后向散射系數(shù)（σ0）描述SAR影像的散射特征，則整個(gè)研究區(qū)域所有像元的σ0組成的樣本集合可以視作水體與背景兩部分的非交并集［11］。記水體像元的后向散射系數(shù)分布概率密度函數(shù)為p（σ0| W），背景像元的后向散射系數(shù)分布概率密度函數(shù)為p（σ0），研究區(qū)域后向散射系數(shù)邊際分布的概率密度函數(shù)為p（σ0），則以上表述可以概括為：

式（1）中，p（W）與p（）分別表示水體與背景像元所占比例，p（W）+p（W）=1。即p（σ0）可視為2類像元的混合分布。該分布形式與影像直方圖的雙峰形走勢是一致的。

1.2 自動(dòng)化閾值尋找算法

統(tǒng)計(jì)學(xué)上認(rèn)為直方圖是對概率密度函數(shù)的估計(jì)［12］。經(jīng)過預(yù)處理后的研究對象，逐個(gè)統(tǒng)計(jì)單一研究區(qū)域內(nèi)像元值的最小值和最大值，在值域內(nèi)定義相應(yīng)對象的后向散射系數(shù)直方圖。直方圖計(jì)算需要的參數(shù)為其每個(gè)條帶的寬度，稱為帶寬。帶寬的計(jì)算公式為［12］：

式（2）中，n為統(tǒng)計(jì)樣本個(gè)數(shù)，IQR為樣本的四分位差。研究表明，由式（2）給出的帶寬適用于Gauss分布樣本數(shù)據(jù)。

對具有明顯峰谷形狀的圖像直方圖來說，閾值法的分割效果較好。如果圖像的直方圖呈現(xiàn)出多峰多谷形狀，則存在多個(gè)直方圖分割閾值，此時(shí)圖像可分割成多個(gè)目標(biāo)區(qū)域。在后向散射系數(shù)直方圖上緊鄰目標(biāo)或背景峰值所發(fā)生畸變的地方可以認(rèn)為是最佳分割閾值?？紤]到水體提取是在全國水利普查成果矢量緩沖區(qū)范圍內(nèi)進(jìn)行，緩沖區(qū)范圍內(nèi)地物類型較少，后向散射系數(shù)直方圖分布較為簡單，表現(xiàn)為典型的雙峰狀。水體像元由于其低散射特性主要分布在直方圖的左側(cè)，因此只需要確定一個(gè)合適的閾值分割目標(biāo)水體與背景，而不用考慮背景中不同地物之間的分割閾值。

以直方圖橫坐標(biāo)第1個(gè)灰度級為初始峰值，逐步比較每個(gè)灰度級所對應(yīng)的像元個(gè)數(shù)，通過循環(huán)迭代遍歷整個(gè)直方圖，尋找直方圖上最大峰值。利用公式（3）將直方圖在灰度級HistSize范圍內(nèi)平均劃分為左、中、右3部分，公式（4）定義了遍歷步長。

式（3）～（4）中，dfmax表示研究區(qū)域內(nèi)像元灰度級的最大值，dfmin表示研究區(qū)域內(nèi)像元灰度級的最小值。

在后向散射系數(shù)直方圖上，目標(biāo)或背景峰值位置分布有以下3種情況，針對不同情況，有相應(yīng)的分割閾值判斷方法。

（1）最大的峰值位置位于直方圖的右側(cè)（直方圖曲線最大峰值位置NmaxIndex＞2HistSize/3），說明研究區(qū)域內(nèi)目標(biāo)水體區(qū)域面積較大，此時(shí)以dstep為步長，從峰值處往左側(cè)遍歷，把直方圖曲線上第1次發(fā)生突變點(diǎn)判斷為分割閾值。

（2）最大的峰值位置位于直方圖的左側(cè)（直方圖曲線最大峰值位置NmaxIndex≤HistSize/3），則說明研究區(qū)域內(nèi)背景區(qū)域面積較大，此時(shí)以dstep為步長，從峰值處往右側(cè)遍歷找出分割閾值。考慮到右側(cè)區(qū)域內(nèi)可能還存在其他背景地物的影響，把距該峰值最遠(yuǎn)處的突變點(diǎn)判斷為分割閾值。

（3）最大的峰值位置位于直方圖的中部（直方圖曲線最大峰值位置HistSize/3 ＜NmaxIndex≤2HistSize/3），如果該峰值右側(cè)仍然存在一個(gè)較明顯的峰值，則該最大峰值為背景峰值，此時(shí)采用步驟（2）中的方法尋找分割閾值，否則認(rèn)為該最大峰值為目標(biāo)峰值，此時(shí)采用步驟（1）中的方法尋找分割閾值。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)區(qū)域與數(shù)據(jù)

為了驗(yàn)證雷達(dá)水體提取算法的有效性，在全國范圍內(nèi)選取了3個(gè)典型湖庫進(jìn)行實(shí)驗(yàn)，分別是云南省陽宗海、江蘇省洪澤湖和北京市密云水庫。綜合考慮影像的標(biāo)稱分辨率、成像帶寬、湖庫覆蓋程度以及數(shù)據(jù)處理難度，精細(xì)條帶、標(biāo)準(zhǔn)條帶模式下的中高分辨率、成像帶寬為50～150 km的影像數(shù)據(jù)更適用于業(yè)務(wù)化水體提取需求。該文選取了3座湖庫的GF-3精細(xì)條帶2影像用于實(shí)驗(yàn)，實(shí)驗(yàn)數(shù)據(jù)如表1所示。

表1 GF-3影像數(shù)據(jù)Table 1 GF-3 images data

2.2 實(shí)驗(yàn)區(qū)域與數(shù)據(jù)

圖1～3分別為陽宗海、洪澤湖和密云水庫的算法處理結(jié)果。圖1a、圖2a和圖3a為裁剪后的后向散射系數(shù)圖，使用的是精細(xì)條帶2模式HH單極化影像?？梢钥闯雠c非水地物相比，水體因其低值而呈現(xiàn)為暗色區(qū)域。圖1b、圖2b和圖3b為水體提取結(jié)果圖，從影像直觀來看，水體范圍得到了有效提取。湖泊、水庫地區(qū)多云霧天氣，尤其是云貴川等地，常年有云，光學(xué)影像難以做到常態(tài)化監(jiān)測，對比之下利用雷達(dá)影像提取水體更具優(yōu)勢。此外，水體提取結(jié)果圖中存在若干噪點(diǎn)，可以通過后續(xù)的形態(tài)學(xué)處理消除。圖1c、圖2c和圖3c為后向散射系數(shù)的統(tǒng)計(jì)分布直方圖，經(jīng)過裁剪后的影像與原始影像相比，區(qū)域內(nèi)水體與背景像元的比例更加均衡，因此直方圖呈現(xiàn)典型的雙峰形走勢，為自動(dòng)閾值尋找算法的有效性提供了保障，其中紅色直線表示自動(dòng)閾值尋找算法所尋閾值。從圖上看，該方法在3幅影像中均較為準(zhǔn)確地找到了直方圖波谷所在的位置，有效分割了影像中的水體與背景像元，證明了自動(dòng)化閾值尋找算法的有效性。

通過尋找直方圖波谷確定最佳分割閾值的方法的合理性在于直方圖的類條件概率密度可正確模擬類的分布的假設(shè)［13］。為了從理論上驗(yàn)證雷達(dá)水體提取算法的有效性，對后向散射系數(shù)的統(tǒng)計(jì)分布模型做進(jìn)一步的估計(jì)。假設(shè)影像中水體與背景像元均服從高斯分布，結(jié)合公式（1），使用Levenberg-Marquardt方法集合直方圖估計(jì)模型參數(shù)得到影像后向散射系數(shù)的理論分布［14］，疊加在直方圖上，顯示為圖中綠色的曲線。如果該曲線反映了后向散射系數(shù)的真實(shí)統(tǒng)計(jì)規(guī)律，那么理論上最佳的分割閾值應(yīng)處在曲線中2個(gè)高斯分布眾數(shù)之間的波谷處。從圖上看，該算法求得閾值所在的紅線與理論分布綠色曲線的波谷位置相近，相距均不超過2 dB，說明了閾值選取的合理性。

圖1 陽宗海實(shí)驗(yàn)區(qū)：a后向散射系數(shù)影像；b水體提取結(jié)果；c后向散射系數(shù)直方圖Fig.1 Yangzonghai Lake area

圖3 密云水庫實(shí)驗(yàn)區(qū)：a后向散射系數(shù)影像；b水體提取結(jié)果；c后向散射系數(shù)直方圖Fig.3 Miyun Reservoir area

2.3 極化方式的選取

在實(shí)驗(yàn)中用于業(yè)務(wù)化水體提取場景的精細(xì)條帶、標(biāo)準(zhǔn)條帶模式數(shù)據(jù)均是多極化方式。該文將進(jìn)一步分析極化方式對水體提取結(jié)果的影響。

為了減少實(shí)驗(yàn)的不確定性，選用陽宗海區(qū)域的全極化條帶模式數(shù)據(jù)，分別嘗試用4種極化影像進(jìn)行水體提取，實(shí)驗(yàn)結(jié)果如圖4。從圖4看，除HH極化影像因成像質(zhì)量問題噪聲較大導(dǎo)致提取效果不佳外，其余3種極化影像均有效地提取出了水體區(qū)域，其差異部分主要集中在水體邊緣地區(qū)，有兩方面的原因：（1）該類區(qū)域內(nèi)植被對像元散射特征的影響；（2）不同極化方式對水陸邊界區(qū)域響應(yīng)不同。

圖4 陽宗海4種極化影像提取結(jié)果：a HH；b HV；c VH；d VVFig.4 Extraction results of 4 polarization images in Yangzonghai Lake

表2 各極化方式影像提取水體面積Table 2 Water area of extraction results in each polarized image

表2展示了各影像所提水體區(qū)域的面積。除HH極化外，其他3種極化方式所提面積均與在GF-1影像上人工標(biāo)注所得的面積相差不大于2%，一定程度上反映出水體提取精度對于極化方式的選取并不敏感。

3 高分三號業(yè)務(wù)化水體提取流程

實(shí)驗(yàn)證明基于單景雷達(dá)影像自動(dòng)提取水體的理論和方法具有可行性，該方法速度快，不需人工選取閾值，具有較高精度，能滿足洪澇監(jiān)測、水資源管理等業(yè)務(wù)化需求。在合適的軟、硬件運(yùn)行環(huán)境下，利用該自動(dòng)化算法實(shí)現(xiàn)業(yè)務(wù)化水體提取的流程如圖5所示。

圖5 業(yè)務(wù)化水體提取流程圖Fig.5 Workflow of operational water body extraction

（1）數(shù)據(jù)訂閱接收：根據(jù)任務(wù)需求向中國資源衛(wèi)星應(yīng)用中心申請訂購監(jiān)測范圍內(nèi)的高分三號影像數(shù)據(jù)，通過專用光纖接收由中國資源衛(wèi)星中心推送的影像數(shù)據(jù)，數(shù)據(jù)為L1A級產(chǎn)品，即斜距復(fù)數(shù)影像。

（2）數(shù)據(jù)預(yù)處理：對獲取的數(shù)據(jù)進(jìn)行數(shù)據(jù)預(yù)處理，包括多視處理、輻射定標(biāo)、地理編碼、濾波和影像裁剪。根據(jù)元信息文件和衛(wèi)星姿軌參數(shù)，經(jīng)過輻射定標(biāo)和地理編碼，斜距復(fù)數(shù)影像轉(zhuǎn)換為地距后向散射系數(shù)影像。利用濾波技術(shù)，消除影像中的斑點(diǎn)噪聲，提高信噪比。利用水利普查數(shù)據(jù)，確定水體的大概范圍，以此為基礎(chǔ)，采用一定規(guī)則確定緩沖區(qū)大小，進(jìn)行影像裁剪，得到一個(gè)主要包含水體信息的矩形影像。

（3）水體自動(dòng)提取：基于上一步裁剪后的影像，利用雷達(dá)水體提取算法進(jìn)行自動(dòng)水體提取工作。

（4）水體提取結(jié)果后處理：利用形態(tài)學(xué)模型進(jìn)行簡單的后處理工作，優(yōu)化水體矢量數(shù)據(jù)。

（5）產(chǎn)品入庫：將經(jīng)過后處理并審核通過的水體提取產(chǎn)品自動(dòng)歸檔數(shù)據(jù)庫中，包括水體的空間矢量數(shù)據(jù)、面積數(shù)據(jù)和相關(guān)的元數(shù)據(jù)，根據(jù)監(jiān)測時(shí)間，如水文數(shù)據(jù)庫（包括該水體），同步提取水文監(jiān)測數(shù)據(jù)，用于后續(xù)分析和展示。

（6）數(shù)據(jù)分析和產(chǎn)品展示：根據(jù)業(yè)務(wù)分析需求，主要完成單個(gè)面積變化分析、區(qū)域水體面積變化分析、結(jié)合水文監(jiān)測數(shù)據(jù)的水量變化分析。設(shè)計(jì)產(chǎn)品展示形式，開發(fā)標(biāo)準(zhǔn)服務(wù)接口，為各業(yè)務(wù)應(yīng)用系統(tǒng)提供接口服務(wù)。

4 結(jié)論

GF-3衛(wèi)星為洪澇監(jiān)測、旱情監(jiān)測等水資源管理業(yè)務(wù)提供了有力的數(shù)據(jù)支撐。基于GF-3影像在水體提取中的巨大應(yīng)用潛力，該文設(shè)計(jì)了一種基于影像直方圖自動(dòng)化閾值分割的水體提取方法。該方法簡單、快速且無需人工干預(yù)，適用于業(yè)務(wù)化水體提取任務(wù)。分別從理論和實(shí)驗(yàn)兩方面驗(yàn)證了該方法的有效性，分析了不同極化方式對水體提取進(jìn)度的影響，發(fā)現(xiàn)4種極化方式下的提取精度均能滿足水體提取任務(wù)的需求，在批量生產(chǎn)中可以根據(jù)實(shí)際情況任選一種極化方式的數(shù)據(jù)提取水體。該文根據(jù)實(shí)際業(yè)務(wù)需求整理了高分三號業(yè)務(wù)化水體提取流程，為自動(dòng)化完成高分三號SAR數(shù)據(jù)從數(shù)據(jù)接收到獲得水體提取結(jié)果提供了規(guī)范性操作指導(dǎo)，為湖庫日常監(jiān)測、洪澇災(zāi)害監(jiān)測等涉水事件提供了有力的技術(shù)支撐，高效的水體自動(dòng)提取也為農(nóng)田水利管理提供了技術(shù)支持。后續(xù)將進(jìn)一步開展基于極化SAR的水體提取研究，提高GF-3衛(wèi)星在水利領(lǐng)域的應(yīng)用水平。