基于Q-OTSU算法高分三號水體快速提取方法研究
——以河南鶴壁洪澇災害監測為例

沈蘭芝，諶 華，范鳳云，王澤宇，于 飛

（二十一世紀空間技術應用股份有限公司，北京 100096）

引言

基于雷達遙感數據的水體信息提取方法主要包括灰度閾值分割法、濾波法、機器學習法和結合輔助信息的提取方法等［1-3］。Klemenjak等［4］結合形態學濾波和監督分類方法，應用TerraSAR-X圖像提取河網，結果表明該方法可以應用于不同的數據集（如多時間、多極化等）。楊存建等［5］采用閾值分割法基于DEM對SAR影像洪水進行提取。Huang等［6］基于隨機森林方法對哨兵一號（Sentinel-1）影像中的水體信息進行提取，實現了較高的自動化程度；谷鑫志等［7］基于GF-3影像，將閾值分割法與馬爾可夫隨機場（MRF）相結合，提出了一種自動化的水體信息提取方法，提取精度較高。欒玉潔等［8-9］采用面向對象閾值分割法對壽光洪災進行監測分析。李晟銘等［10］基于Sentinel-1影像對2016年長江中下游重災區洪水進行了監測及災情評估。湯玲英等［11］基于面向對象方法研究Sentinel-1數據在洪水監測中的應用。辜曉青等［12］利用Sentinel-1影像監測了2018年6月25日江西撫順暴雨的空間分布信息。崔倩等［13］利用GF-3影像對全國3個典型湖庫地區進行了水體提取。周帆等［14］對GF-3和Sentinel-1影像在斯里蘭卡2017年5月底發生的洪災洪水信息提取結果進行對比。這些方法中，濾波法可以較好地抑制SAR圖像的斑點噪聲，但得到的水體邊緣線不夠連續光滑，易引起邊緣特征的誤差，檢測精度低［3］；機器學習方法可針對特定情境建立模型，提取精度高，但獲取樣本及特征的代價大，分類器訓練時間花費時間長［15］；而灰度閾值分割法速度快、原理簡單、計算量小。以上研究雖取得了一些成果，但方法相對復雜，在實際應急情況中，無法協同精度和效率，在保證精度的同時，較快地完成對大范圍洪水信息的快速提取與決策分析，難以滿足應急時效性的需求。

高分三號衛星由中國自主研發，所搭載的SAR傳感器具有大成像寬幅、超高分辨率、多成像模式、長壽命運行等特點［16］，能夠全天候和全天時實現全球海洋和陸地信息的監視監測，不僅可以大范圍普查，也能對特定目標進行詳查，是世界上工作模式最多分辨率最高的C頻段多極化合成孔徑雷達衛星［17］。

本文基于GF-3影像，對2021年7月中下旬河南省鶴壁市的洪澇災害進行洪水遙感監測，同時利用同期哨兵二號（Sentinel-2B）影像進行精度驗證，利用Q-OTSU算法快速提取洪澇災情信息，完成專題圖制作，為河南省抗洪救災提供了有力支撐。

1 研究區和數據介紹

1.1 研究區概況

研究區為河南省鶴壁市，地理坐標113°59'-114°45'E，35°26'-36°02'N，位于河南省北部，太行山東麓向華北平原過渡地帶，總面積為2 182平方公里，常住人口為1 565 973萬人。鶴壁市屬于暖溫帶半濕潤型季風氣候，年均氣溫14.2～15.5℃，年際降水變化大，多年平均降水量為571 mm，汛期6～9月份平均降雨量占全年降雨量的70%左右［18］。

境內主要河流共有13條，其中流域面積超過1 000平方公里的河流有3條（淇河、衛河、共產主義渠），3條河流境內總長192.5公里。共有6個蓄滯洪區，總滯洪量5.87億立方米，洪區總面積327平方公里。共有14座水庫，總庫容6.37億立方米，其中最大水庫盤石頭水庫，總庫容達6.08億立方米［19］。

1.2 數據介紹

研究主要采用鶴壁區GF-3 L1A級數據（下載自http://36.112.130.153:7777/DSSPlatform/index.html）及Sentinel-2 L2A級數據（下載自https://search.asf.alaska.edu/），包括18景精細條帶1成像模式（FSII）的HH/HV雙極化GF-3影像及1景Sentinel-2B影像。18景GF-3影像成像時間從7月20日至8月30日，共包含12期。成像信息詳見表1。

表1 鶴壁市影像信息Table 1 Hebi image information

FSII模式GF-3影像分辨率10 m，幅寬100 km，入射角為27.8°。Sentinel-2B影像，包括13個光譜波段，分辨率為10 m、20 m和60 m，幅寬達290 km［20］。另外，采用30 m分辨率ASTER DEM及河南省鶴壁矢量文件作為輔助數據。

為驗證Q-OTSU方法的普適性，采用河南省其他3個區域的GF-3 L1A級數據影像進行水體信息提取，并采用相近日期的Sentinel-2影像做精度驗證。影像詳細信息見表2。

表2 其他區域影像信息Table 2 Image information of other areas

2 研究方法

HV極化的穿透性要比HH極化的強，對水體信息識別效果更好，因此，本研究首先對各期HV極化GF-3影像進行批量定量化預處理，轉換到地理坐標系下（圖1）。之后通過Q-OTSU算法對影像進行閾值分割，初步得到水體信息圖。再對水體信息進行小斑點去噪，得到最終水體信息分布圖，快速制作水體信息成果圖。最后選取與GF-3影像成像時間相近的Sentinel-2B光學影像作為參考，評價水體提取精度。

圖1 水體提取總流程圖Fig.1 The general flow chart of water extraction

2.1 數據預處理

首先對GF-3影像進行定量化預處理，依次進行輻射定標、多視、Lee濾波、地理編碼、拼接、db格式轉換及裁剪操作。GF-3原始數據采用兩個波段將實部與虛部存儲在一個文件中，從元數據中獲取定標參數，分別對實部和虛部進行輻射定標，組合為單視復數（SLC）數據，將原始信號轉化為雷達后向散射系數［7］。結合雷達入射角和距離向、方位向的幾何關系，確定采用距離向3視、方位向3視進行多視化處理。相干斑噪聲，是SAR影像的固有特點，嚴重影響了SAR影像的地物可解譯性，采用Lee濾波進行斑點噪聲濾波處理。采用空間分辨率為30 m的ASTER DEM數據進行地理編碼，將GF-3影像從雷達坐標系轉到地理坐標系。對于單景影像不能全覆蓋研究區的部分數據，利用同期多景影像進行拼接，保證拼接后影像盡可能多地覆蓋全研究區域。之后經過取對數運算轉換到db格式，并按河南省鶴壁市矢量范圍裁剪，最終生成空間分辨率為10 m的SAR正射影像產品，完成自動化預處理過程。

Sentinel-2B L2A級數據已經過輻射定標和大氣校正操作，數據產品中為13個單波段產品。選擇其中4個10 m分辨率的波段B2（Blue）、B3（Green）、B4（Red）及B8（NIR）進行波段合成，并按河南省鶴壁市矢量范圍裁剪，形成10 m分辨率Sentinel 2B預處理產品。

2.2 Q-OTSU算法

大津法（OTSU）由日本學者大津于1979年提出，又稱為最大類間方差法［21-22］，是一種確定圖像二值化分割閾值的算法。該方法所選閾值應滿足類間方差最大，類內方差最小的準則［23］。

Q-OTSU在影像分割前進行8鄰域均值濾波，相比OTSU更多地考慮了圖像空間鄰域信息，可以有效濾除噪聲影響；在分割時基于OTSU選擇的閾值重新自適應選擇最優閾值，使得最終選擇的閾值更接近于人工選取的閾值，提高水體提取精度；在分割后進行小斑點去噪，進一步去除由噪聲導致的誤提水體，提高精度。

本研究采用Q-OTSU算法進行研究區水體信息提取，具體流程見圖2。先對預處理后的GF-3影像進行8鄰域均值濾波操作，之后將影像像元值轉化到［0，255］范圍，轉換公式為：

圖2 Q-OTSU算法水體提取流程圖Fig.2 The flow chart of water extraction based on Q-OTSU algorithm

式中：F（x）為影像中x點轉換后的像元值；f（x）為影像中x點原始的像元值；fmax（x）為影像中原始最大像元值；fmin（x）為影像中原始最小像元值。

然后選擇自適應最優分割閾值。先設定臨時閾值，臨時閾值獲取公式為：

式中：x為臨時閾值對應的像元值，取值范圍為［0，255］，σ2(x)為水體與非水體間的類間方差，σ2(x)的獲取方式為：

等價于：

式中：ω1和ω2分別為水體和非水體像元所占比重；μ1和μ2分別為水體和非水體的像元均值；μ為總體像元均值。

OTSU算法選取的閾值為t，t為M（x）取最大值時對應的像元值：

基于OTSU算法選取的閾值t，Q-OTSU算法最優閾值T的獲取公式為：

式中：s為圖像像元直方圖統計結果呈現雙峰狀情況下（圖3），兩峰之間谷底所對應的像元值；t為OTSU算法選取的閾值。

圖3 最優閾值選取Fig.3 Optimal threshold selection

獲得最優分割閾值T后，根據水體像元后向散射系數特征，取像元值低于閾值T的像元為水體像元，否則，為非水體像元，得到初步水體信息圖。對初步水體信息分布圖進行小斑塊去噪處理，去除面積小于10個像元的小斑塊，得到最終水體信息分布圖。

本研究中Q-OTSU算法通過自主編輯的python代碼實現。

2.3 水體提取精度驗證指標

為定量地評價本研究中所用水體信息提取方法的精度，選取了研究區7月31日成像時間的Sentinel-2B影像，對同期GF-3影像水體信息提取結果做精度驗證。在Sentinel-2B影像上研究區范圍內隨機選擇300個檢驗樣本（圖4），結合GF-3假彩色合成影像、Sentinel-2B NDWI影像等目視解譯識別樣本點的水體-非水體類別屬性，作為后期實驗結果的精度評價參考數據。

圖4 水體提取樣點分布Fig.4 Distribution of water extraction samples

圖4中顯示了隨機樣點的分布情況，底圖為Sentinel-2假彩色合成影像。比較樣點位置在實驗GF-3影像上的水體提取結果和Sentinel-2參考影像上的水體-非水體類別信息，若類別一致則定義為正確的檢驗樣本（Correct Sample，CS），不一致則定義為錯誤的檢驗樣本（Wrong Sample，WS）。計算混淆矩陣，用混淆矩陣中的總體分類精度（Overall Accuracy，OA）代表研究水體信息提取精度。

3 實驗結果分析

3.1 水體提取結果精度驗證

為了驗證本文提出的Q-OTSU算法在水體提取方面的有效性，利用傳統OTSU算法和Q-OTSU算法分別對7月31日河南省鶴壁市預處理GF-3影像進行水體提?。▓D5）。從整體上對比可知，傳統OTSU算法提取出來的水體信息相對于Q-OTSU算法提取出的水體有很多分散的小噪點，尤其是鶴壁市中部和西部。經后期仔細比對原始影像，發現這些小噪點多為其他地物錯分為水體。

圖5 水體分布圖Fig.5 Water distribution map

為進一步對比Q-OTSU算法的水體提取優勢，選取其中一小塊區域做詳細對比，選取的區域為圖中綠色框標出的區域。對比發現，在傳統OTSU算法水體提取結果中，部分柏油道路被錯分為水體，且非水體區域有少量斑點噪聲，這些斑點噪聲多為比較大的十字路口或彩鋼棚被誤分為水體。原因是這些被誤分的地物表面也相對比較光滑，后向散射能力弱，它們的像素值與水體信息接近，導致被誤分。還有部分斑點噪聲為影像本身噪聲引起。而在Q-OTSU水體提取結果中，道路分類基本正確，且非水體區域的斑點噪聲現象得到很大改善。

利用Sentinel-2光學影像上采集的隨機樣點和分類精度指標OA，定量評價本研究中提出的Q-OTSU水體提取算法精度。表3列出了傳統OTSU算法及Q-OTSU算法水體提取精度信息。Q-OTSU算法的水體提取精度比傳統OTSU算法的精度高，達到了96.7%，精度要求滿足應急水體提取所需。

表3 精度評價表Table 3 Accuracy evaluation table

另外，本研究中，從數據獲取、下載，利用服務器對每景GF-3影像預處理，時間約3 h，水體提取時間約為9 min，在接到影像24 h內高效完成影像水體信息提取成果制圖并進行洪水分布分析，編制分析報告。

為進一步驗證Q-OTSU方法的普適性，選取河南省其他3個不同時段不同區域的GF-3影像進行QOTSU方法水體信息提取，并采用時間相近的Sentinel-2影像做精度驗證。分別在3個區域范圍內隨機選擇200個檢驗樣本（圖6（a）、圖7（a）、圖8（a）），結合Sentinel-2假彩色影像、Sentinel-2 NDWI影像等目視解譯識別樣本點的水體-非水體類別屬性，作為對應GF-3影像水體提取結果的精度評價參考數據。

圖6 區域1樣點分布及水體分布圖Fig.6 Sample points map and water distribution map in area 1

圖7 區域2樣點分布及水體分布圖Fig.7 Sample points map and water distribution map in area 2

圖8 區域3樣點分布及水體分布圖Fig.8 Sample points map and water distribution map in area 3

圖6（b）、圖7（b）和圖8（b）中，水體信息已基本完整提取。對3個區域的水體信息提取結果分別做精度驗證（見表4），發現精度都在93%以上，進一步驗證了Q-OTSU方法對水體提取的普適性。

表4 其他區域精度評價表Table 4 Accuracy evaluation table of other areas

3.2 多時域水體信息變化檢測分析

本研究通過Q-OTSU算法分別對河南省鶴壁市7月20日至8月30日共12期GF-3影像進行水體信息提取，各期水體信息分布如圖9。為了量化分析受災期間鶴壁市及其各縣區水體信息變化情況，將鶴壁市及其鶴山區、浚縣、淇濱區、淇縣及山城區的水體信息分布面積做了相關統計（圖10）。8月23日與8月30日的GF-3影像不能全覆蓋研究區域，不參與水體面積統計。根據應急需求，本文將7月20日影像作為災前影像。

結合圖9和圖10可知，整體上，鶴壁市的水體分布范圍呈現先增大后減小的趨勢，在7月31日達到最大（圖10），且受災區域主要集中在浚縣和淇縣。7月20日，鶴壁市水體面積15.45 km2，水體信息基本分布在淇河、衛河、共產主義渠、盤石頭水庫及其他一些小水庫或小河流等河流水庫周圍（圖9（a））。連續大暴雨后，7月25浚縣及淇縣南部已有大量積水（圖9（b））。7月31日，水體分布范圍在7月28日的基礎上進一步擴大（圖9（d）），全市水體分布面積已高達209.44 km2，?？h東南部、共產主義渠周圍及淇縣東南部已基本被積水覆蓋。其中?？h是鶴壁市受災最嚴重的一個區縣，7月31日水體分布面積達到170 km2（圖9），是其他區縣水體分布面積總和的4倍多。至8月28日，暴雨過后將近40天，積水明顯大范圍褪去（圖9（k）），但水體分布面積仍高于災前7月20日，為35.9 km2，洪水還未完全退去。

圖9 水體信息分布圖Fig.9 Water information distribution map

圖10 鶴壁市水體信息變化圖Fig.10 Water information change map of Hebi

為了更好地顯示此次降雨造成的洪水淹沒范圍，選取積水分布范圍最大的7月31日水體分布圖與災前7月20日水體分布圖做對比（圖9）。圖11中，藍色區域為災前水體分布，紅色區域為受災洪水淹沒區域，主要集中在?？h和淇縣，受災淹沒區域總面積為193.99 km2。淇縣東南部和?？h地勢低，被淹沒的大部分地區為農田和村莊，洪水不易退去。鶴山區、山城區和淇濱區大部分地區為山地，地勢高，不易積水，災前災后水體分布面積無明顯變化。

圖11 災前災后水體信息分布對比圖Fig.11 Comparison of water information distribution before and after disaster

4 結論

本文以12期國產高分三號衛星10 m分辨率影像為主要數據源，輔以1期10 m分辨率哨兵二號影像及30 m分辨率DEM數據，利用Q-OTSU算法對2021年7月中下旬因暴雨導致特大洪災的河南省鶴壁市進行連續40天洪澇災害遙感監測，對災情期間鶴壁市及其5個縣區的水體分布變化情況進行分析。研究表明，Q-OTSU算法相比傳統OTSU算法，可改善道路、彩鋼棚等地物被誤分為水體和非水體區域存在斑點噪聲的現象，自動化程度高，水體提取精度高達96.7%。在保證精度的同時提取效率也高，在獲取影像的24 h內可完成影像水體提取結果出圖并進行災情分析，在鶴壁市洪災的水體提取應急數據處理中具有良好的使用性。另外，此次洪災中，鶴壁市的淇縣和浚縣由于地勢相對其他縣區較低，且被淹沒區域多為農田和村莊，受泄洪影響，排水不易，受災嚴重。GF-3影像可高效開展洪水災情快速評估及災后重建，為政府和相關部門及時、精確地掌握洪災信息，科學防災救災提供了決策依據。

通過對水體分布圖的仔細檢查對比，發現當遇到細小水體時水體與背景像元不易區分，由于難以識別會出現河流間斷或小面積水體漏分等現象，影響水體提取精度。在之后的研究中應更多地選擇圖像特征以及關注水體像元與鄰近背景像元的關系，實現水體信息連貫提取，提高水體提取精度。