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

（二十一世紀(jì)空間技術(shù)應(yīng)用股份有限公司，北京 100096）

引言

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

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

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

1 研究區(qū)和數(shù)據(jù)介紹

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)介紹

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

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

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

為驗(yàn)證Q-OTSU方法的普適性，采用河南省其他3個(gè)區(qū)域的GF-3 L1A級(jí)數(shù)據(jù)影像進(jìn)行水體信息提取，并采用相近日期的Sentinel-2影像做精度驗(yàn)證。影像詳細(xì)信息見表2。

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

2 研究方法

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

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

2.1 數(shù)據(jù)預(yù)處理

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

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

2.2 Q-OTSU算法

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

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

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

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

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

然后選擇自適應(yīng)最優(yōu)分割閾值。先設(shè)定臨時(shí)閾值，臨時(shí)閾值獲取公式為：

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

等價(jià)于：

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

OTSU算法選取的閾值為t，t為M（x）取最大值時(shí)對應(yīng)的像元值：

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

式中：s為圖像像元直方圖統(tǒng)計(jì)結(jié)果呈現(xiàn)雙峰狀情況下（圖3），兩峰之間谷底所對應(yīng)的像元值；t為OTSU算法選取的閾值。

圖3 最優(yōu)閾值選取Fig.3 Optimal threshold selection

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

本研究中Q-OTSU算法通過自主編輯的python代碼實(shí)現(xiàn)。

2.3 水體提取精度驗(yàn)證指標(biāo)

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

圖4 水體提取樣點(diǎn)分布Fig.4 Distribution of water extraction samples

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

3 實(shí)驗(yàn)結(jié)果分析

3.1 水體提取結(jié)果精度驗(yàn)證

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

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

為進(jìn)一步對比Q-OTSU算法的水體提取優(yōu)勢，選取其中一小塊區(qū)域做詳細(xì)對比，選取的區(qū)域?yàn)閳D中綠色框標(biāo)出的區(qū)域。對比發(fā)現(xiàn)，在傳統(tǒng)OTSU算法水體提取結(jié)果中，部分柏油道路被錯(cuò)分為水體，且非水體區(qū)域有少量斑點(diǎn)噪聲，這些斑點(diǎn)噪聲多為比較大的十字路口或彩鋼棚被誤分為水體。原因是這些被誤分的地物表面也相對比較光滑，后向散射能力弱，它們的像素值與水體信息接近，導(dǎo)致被誤分。還有部分斑點(diǎn)噪聲為影像本身噪聲引起。而在Q-OTSU水體提取結(jié)果中，道路分類基本正確，且非水體區(qū)域的斑點(diǎn)噪聲現(xiàn)象得到很大改善。

利用Sentinel-2光學(xué)影像上采集的隨機(jī)樣點(diǎn)和分類精度指標(biāo)OA，定量評價(jià)本研究中提出的Q-OTSU水體提取算法精度。表3列出了傳統(tǒng)OTSU算法及Q-OTSU算法水體提取精度信息。Q-OTSU算法的水體提取精度比傳統(tǒng)OTSU算法的精度高，達(dá)到了96.7%，精度要求滿足應(yīng)急水體提取所需。

表3 精度評價(jià)表Table 3 Accuracy evaluation table

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

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

圖6 區(qū)域1樣點(diǎn)分布及水體分布圖Fig.6 Sample points map and water distribution map in area 1

圖7 區(qū)域2樣點(diǎn)分布及水體分布圖Fig.7 Sample points map and water distribution map in area 2

圖8 區(qū)域3樣點(diǎn)分布及水體分布圖Fig.8 Sample points map and water distribution map in area 3

圖6（b）、圖7（b）和圖8（b）中，水體信息已基本完整提取。對3個(gè)區(qū)域的水體信息提取結(jié)果分別做精度驗(yàn)證（見表4），發(fā)現(xiàn)精度都在93%以上，進(jìn)一步驗(yàn)證了Q-OTSU方法對水體提取的普適性。

表4 其他區(qū)域精度評價(jià)表Table 4 Accuracy evaluation table of other areas

3.2 多時(shí)域水體信息變化檢測分析

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

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

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

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

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

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

4 結(jié)論

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

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