基于高分三號衛星數據的水體自動提取及應用

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

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

0 引言

我國洪澇和干旱災害和水資源短缺等水問題突出，地面監測系統目前還不能實現全面覆蓋。高分三號衛星具有全天候全天時對地觀測能力，可以不分晝夜、陰雨天氣進行對地觀測。針對洪澇災害監測、水資源管理等水利應用需求，基于可以常態化接收的國產雷達數據，實現水體的自動化提取，可以補充地面監測能力的不足，在提升我國水利監測能力方面發揮了重要作用。

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

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

出于自動化、業務化的考慮，所應用的雷達水體提取方法應具有簡單、穩定、快速的特點，并具有滿足業務需求的提取精度。結合上述常用水體提取方法，文章采用雙峰法進行業務化水體提取工作。直觀上說，雙峰法尋找閾值的過程，是根據雷達影像后向散射系數統計分布直方圖呈現出的雙峰值特征、在兩峰之間尋找波谷的過程。

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

1 雷達水體自動化提取算法

1.1 后向散射系數分布模型

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

對于水體來說，未擾動水面的后向散射系數非常低，在所有雷達波段都接近于鏡面反射，微波能量對于高介電常數幾乎以等于入射角的方向反射出去。就粗糙度看，平靜水面相對于雷達波來說幾乎是水平的，因此水的雷達信號非常弱，區別于其他地物，在影像中呈現為低后向散射系數的均勻區域。在水體提取任務中，通常只對某個像元是否屬于水體感興趣，對于非水體的背景像元的具體地物類型并不在意。考慮到水體像元的低后向散射特性，將像元劃分為水體與背景2類，并進一步討論這2類像元的特征。

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

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

1.2 自動化閾值尋找算法

統計學上認為直方圖是對概率密度函數的估計［12］。經過預處理后的研究對象，逐個統計單一研究區域內像元值的最小值和最大值，在值域內定義相應對象的后向散射系數直方圖。直方圖計算需要的參數為其每個條帶的寬度，稱為帶寬。帶寬的計算公式為［12］：

式（2）中，n為統計樣本個數，IQR為樣本的四分位差。研究表明，由式（2）給出的帶寬適用于Gauss分布樣本數據。

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

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

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

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

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

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

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

2 實驗與分析

2.1 實驗區域與數據

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

表1 GF-3影像數據Table 1 GF-3 images data

2.2 實驗區域與數據

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

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

圖1 陽宗海實驗區：a后向散射系數影像；b水體提取結果；c后向散射系數直方圖Fig.1 Yangzonghai Lake area

圖3 密云水庫實驗區：a后向散射系數影像；b水體提取結果；c后向散射系數直方圖Fig.3 Miyun Reservoir area

2.3 極化方式的選取

在實驗中用于業務化水體提取場景的精細條帶、標準條帶模式數據均是多極化方式。該文將進一步分析極化方式對水體提取結果的影響。

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

圖4 陽宗海4種極化影像提取結果：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展示了各影像所提水體區域的面積。除HH極化外，其他3種極化方式所提面積均與在GF-1影像上人工標注所得的面積相差不大于2%，一定程度上反映出水體提取精度對于極化方式的選取并不敏感。

3 高分三號業務化水體提取流程

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

圖5 業務化水體提取流程圖Fig.5 Workflow of operational water body extraction

（1）數據訂閱接收：根據任務需求向中國資源衛星應用中心申請訂購監測范圍內的高分三號影像數據，通過專用光纖接收由中國資源衛星中心推送的影像數據，數據為L1A級產品，即斜距復數影像。

（2）數據預處理：對獲取的數據進行數據預處理，包括多視處理、輻射定標、地理編碼、濾波和影像裁剪。根據元信息文件和衛星姿軌參數，經過輻射定標和地理編碼，斜距復數影像轉換為地距后向散射系數影像。利用濾波技術，消除影像中的斑點噪聲，提高信噪比。利用水利普查數據，確定水體的大概范圍，以此為基礎，采用一定規則確定緩沖區大小，進行影像裁剪，得到一個主要包含水體信息的矩形影像。

（3）水體自動提取：基于上一步裁剪后的影像，利用雷達水體提取算法進行自動水體提取工作。

（4）水體提取結果后處理：利用形態學模型進行簡單的后處理工作，優化水體矢量數據。

（5）產品入庫：將經過后處理并審核通過的水體提取產品自動歸檔數據庫中，包括水體的空間矢量數據、面積數據和相關的元數據，根據監測時間，如水文數據庫（包括該水體），同步提取水文監測數據，用于后續分析和展示。

（6）數據分析和產品展示：根據業務分析需求，主要完成單個面積變化分析、區域水體面積變化分析、結合水文監測數據的水量變化分析。設計產品展示形式，開發標準服務接口，為各業務應用系統提供接口服務。

4 結論

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