基于Sentinel-1/2數據的洪水淹沒范圍提取模型

摘要：遙感是監測洪水淹沒范圍、掌握洪澇災情演變的重要手段，而光學影像在洪水發生時往往有較多缺失，全天候的SAR影像在提取水體時精度略低。為快速、精準提取洪水淹沒范圍，構建了一種綜合利用Sentinel-2光學影像和Sentinel-1雷達影像數據的洪水淹沒范圍提取模型，采用一種自適應閾值分割算法即大津算法（OTSU）分別對兩種數據以及該模型進行了水體范圍提取試驗，并以河北省保定市為例進行了應用分析。結果顯示：云量較少的Sentinel-2影像水體提取效果最好，總體精度（OA）達到95.6%；所構建的模型在引入部分可用Sentinel-2數據后，OA達到95%，相比單獨使用Sentinel-1數據OA和Kappa系數分別提升1.2%和2.4%。該模型搭載于Google Earth Engine平臺，能實現快速、準確、低成本的地表水體空間范圍連續輸出，不受限于云霧且比單獨使用Sentinel-1影像的提取精度更高，在云覆蓋嚴重導致Sentinel-2數據缺少的情況下，該模型可作為洪水淹沒范圍提取方法的一種選擇。

關 鍵 詞：洪水淹沒范圍； Sentinel-1； Sentinel-2； 自適應閾值分割算法； Google Earth Engine； 保定市

中圖法分類號： P954

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.09.010

0 引 言

由于地理位置和東亞季風的原因，中國一直飽受洪澇災害的侵擾。據《中國水旱災害防御公報2022》概要中統計［1］，2022年因洪澇災害造成171人死亡失蹤，直接經濟損失達1 288.99億元（占當年自然災害總直接經濟損失的54.01%）。及時、大范圍地掌握洪澇災情的演變情況，對災情評估、災后重建有重要的指導意義。

基于水的光譜特性，方便、快速的水體指數模型法可以在不同復雜情景的光學影像中取得較高的水體識別精度［2-4］，可運用于水體的動態監測［5］。但受限于云霧，洪災發生時段實際可用的光學影像較少，因此在監測洪水面積時很難將光學影像作為唯一數據來源。SAR影像具有大范圍、全天候、短周期、長時序的優點，不僅被廣泛應用于城市礦山地面沉降［6］，在洪水動態連續監測方面也具備獨特的優勢。基于上述背景，國內外學者對洪水提取模型進行了廣泛研究。例如賈詩超等［7］分析微波遙感水體信息的影像特點后，提出SDWI指數模型增強水體特征，提取效果顯著；Amitrano等［8］提出了一種基于Sentinel-1數據的無監督洪水提取模型，能夠快速提取洪水發生范圍；Binh等［9］結合SAR數據（Sentinel-1）和光學數據（Landsat-8）的訓練數據開發了一種神經網絡分類模型，并運用于湄公河三角洲，證明該模型有較強魯棒性和較高準確率。然而上述模型也有一定缺陷：Amitrano等所建模型雖然可以快速提取洪水發生范圍，但數據源單一，精度較低；Binh等所建模型雖取得了較高的精度，但訓練模型需要花費較長時間和較大人力，很難運用于強調時效性的抗洪救災場景中。與此同時，數據的獲取和必要的預處理也要盡可能方便，以提高計算效率，而Google Earth Engine（GEE）平臺正好滿足這一需求。GEE是向全球用戶開放，可在線處理全球尺度衛星遙感數據并進行計算和可視化分析的平臺，通過編碼即可遠程調用Landsat系列、Sentinel系列等免費數據進行處理而不占用計算機資源，極大提高了用戶的效率，在水體監測方面得到了廣泛應用［10-12］。

綜上所述，本研究考慮綜合利用光學（Sentinel-2）和SAR（Sentinel-1）兩類數據，理由如下：一方面，覆蓋范圍廣、全天候的Sentinel-1數據可以補足Sentinel-2數據的缺失；另一方面，部分存有的Sentinel-2數據可以提升模型的精度。在方法上選擇不依賴前期訓練、能夠快速得到水體范圍的自適應閾值分割法，搭載于GEE平臺來滿足洪水監測的時效性要求。最后，通過實驗驗證單獨使用Sentinel-1/2數據以及兩者綜合利用構建的本文模型在不同場景的適用性和準確度并進行對比分析。研究成果以期實現洪水范圍的長時序、全自動、高精度提取。

1 研究區概況及數據預處理

1.1 研究區概況

研究區位于河北省保定市東部，包含6個縣級行政區，地理位置約在115°30′E～116°15′E，38°42′N～39°36′N之間。區域內地表永久性水體主要為位于涿州市的北拒馬河與位于安新縣的白洋淀，除了建成區以外，境內大部分土地類型為耕地，通過實地走訪調查得知，種植的主要作物為玉米、小麥。2023年受到臺風“杜蘇芮”等多重因素影響，河北省多地發生特大暴雨從而造成洪水，研究區作為受洪災影響最嚴重的區域，引起了社會廣泛關注。

1.2 數據預處理

本研究SAR影像數據采用Sentinel-1A數據，數據分辨率為10 m，按工作模式為干涉寬幅模式（IW），數據類型為地距多視（GRD），極化方式為VV和VH雙極化，成像軌道為升軌的條件篩選影像。隨后按研究區范圍對影像進行裁剪并用Refined Lee濾波［13］消除一定的SAR影像相干斑噪聲。光學影像數據使用Sentinel-2數據，由多光譜成像儀（MSI）成像，包含可見光、近紅外、短波紅外等13個波段，分辨率包含10，20，60 m三種，本文選用10 m分辨率以獲得更高的提取精度。光學影像去云處理采用GEE平臺自帶的云量分數S2_CLOUD_PROBABILITY數據集，該數據集通過s2cloudless算法計算每一景Sentinel-2影像像元的云量分數，像元云量分數越大，含云的概率越大。本研究所有Sentinel-2影像數據均按照云量分數小于50的準則進行去云處理。經過預處理，得到覆蓋整個研究區的多幅Sentinel-1/2影像，每幅影像的數據類型及其對應的時間信息如表1所列。最后，SAR側視成像還存在山體陰影和水體相互混淆的問題，因此利用美國航天飛機雷達地形測繪計劃SRTM（Shuttle Radar Topography Mission）DEM數據計算坡度柵格［14］，設置坡度閾值對Sentinel-1數據水體范圍進行掩膜以去除山體陰影的影響。

2 原理與方法

閾值法分割圖像的關鍵在于閾值的準確度。在地表水體提取的情景中，當時間跨度較大時，同一地理位置的大氣條件、地面作物的生長階段等不同會引起最佳分割閾值的變化，此時利用固定的閾值分割會導致較為嚴重的誤分。大津算法（OTSU）作為閾值算法的一種，具備計算快速、準確率較高的優點，相比于直方圖雙峰法受限于明顯的雙峰特性，其在閾值的選取上脫離人工分析，能夠滿足大批量遙感影像洪水淹沒范圍自動提取的需求。

2.1 自適應閾值分割原理

大津算法，又被稱為最大類間方差法［15］，通過遍歷圖像的灰度值，尋找使得圖像的前景和背景的類間方差最大的閾值，從而實現二值化閾值分割。其具體計算方式如下：

設圖像的灰度空間為［0，T］，某個灰度級ti的像素個數為f（ti），N為圖像像素的個數，則該灰度級出現的概率Pi為

Pi=ftiN（1）

設灰度級t為當前閾值，圖像被分為A、B兩類，兩類像素的概率PA、PB及均值灰度MA、MB分別為

PA=ti=0Pi（2）

PB=Ti=tPi=1-PA（3）

MA=ti=0i（Pi/PA）（4）

MB=Ti=ti（Pi/PB）（5）

整幅影像的均值灰度M為

M=Ti=0tiPi=MAPA+MBPB（6）

兩類的類間方差δ2為

δ2=PAMA-M2+PBMB-M2

=PAPBMA-MB2（7）

遍歷圖像的灰度空間使得δ2最大的灰度值即為需要求得的最佳閾值。

2.2 Sentinel-1數據的水體提取方法

SAR衛星傳感器發射的微波對于地表目標的散射特性為：在粗糙表面（建筑、裸土等）的反射為漫反射，回波信號強，后向散射系數大，在影像中呈現亮白色；在光滑表面（水體等）的反射為鏡面反射，回波信號弱，后向散射系數較小，在影像中呈現為暗黑色。因此可以選擇合適的閾值來提取SAR影像中的水體。

Sentinel-1A影像由VV和VH兩種極化方式成像，在雙重反彈效應（電磁波在水面散射到部分被淹沒的樹干等處，然后直接反射回SAR傳感器）下，VV極化相比VH極化對表面粗糙度更為敏感，并顯示出更強烈的信號響應，有利于提取洪水中的半淹沒區域，具有更好的洪水監測適用性［16］。因此，本研究采用OTSU算法自動提取VV極化SAR圖像的分割閾值完成水體提取，即VV-OTSU的Sentinel-1A影像水體自動提取方法。

2.3 Sentinel-2數據的水體提取方法

光學影像的水體指數法是水體提取的一種常用方法。其原理是基于水體對綠光波段的高反射率以及近紅外光波段的高吸收率特性，通過比值運算構建水體指數來增強水體在影像中的特征，而后使用簡單的分割算法完成水體范圍的提取。

水體指數中NDWI最為常見和流行［17］，其計算公式如下：

NDWI=（ρ3-ρ8）/（ρ3+ρ8）（8）

式中：ρ3、ρ8分別為Sentinel-2數據中綠光波段B3和近紅外波段B8的反射率值。

Xu［18］通過研究改進NDWI指數提出了MNDWI，可以在水、建筑物以及植被混合的區域取得更高的水體提取精度。由于短波紅外波段中水的吸收比近紅外波段強，因此使用短波紅外波段代替近紅外波段。MNDWI的計算公式如下：

MNDWI=（ρ3-ρ11）/（ρ3+ρ11）（9）

式中：ρ3、ρ11分別為Sentinel-2數據中綠光波段B3和短波紅外波段B11的反射率值。

需要注意的是，Sentinel-2數據中B3、B8波段的分辨率為10 m，而B11波段的分辨率為20 m，在計算MNDWI之前先將B11波段使用最鄰近重采樣至10 m。為了比較NDWI和MNDWI在研究區內的適用性，在同一時期的兩個不同地點的Sentinel-2 NDWI/MNDWI圖像上分別使用OTSU閾值對水體提取范圍進行比較。圖1顯示了兩者在典型區域的視覺結果，可以看出NDWI在建筑和水體混合的拒馬河區域識別出了多余的噪聲斑點，并且將白洋淀外圍的沙丘誤提為水體，整體表達效果不如MNDWI精細、簡潔。其余部位也反映出類似的情況，因此選擇MNDWI-OTSU作為Sentinel-2數據水體提取方法。

3 基于Sentinel-1/2數據的洪水淹沒范圍快速提取模型構建

為了得到準確的水體提取范圍，本文首先開展了不同方法提取效果對比研究。根據對比結果，構建了綜合利用Sentinel-1/2數據的洪水淹沒范圍快速提取模型，并對模型精度進行了評價。

3.1 MNDWI-OTSU與VV-OTSU水體提取效果對比

同時期的地面水體真實范圍難以實際測量［18］，本研究采用目視解譯方式對研究區的水體進行解譯，作為精度驗證的基準。Sentinel-2 FRB（Full Resolution Browse）數據是由短波紅外、近紅外、紅波段進行合成，并經過壓縮和拉伸（gamma = 2.0）生成的優化高分辨率圖像。該圖像在不剪切極值的情況下能夠突出植被和水體的差異［19］，目視解譯效果較好，適用于研究區耕地、濕地、水體混合的復雜情景，如圖2所示為研究區2023年7月16日的Sentinel-2 FRB影像。因此本研究選擇Sentinel-2 FRB圖像進行目視解譯并與水體提取結果進行直觀對比，從而評價本文方法的水體提取精度。Sentinel-2 FRB的時間信息與本研究使用的Sentinel-2影像一一對應，并且同樣經過去云處理，具體見表1。

由于研究區內永久性水體占比較小，為了減少多余計算量，提高計算效率，在研究區內圈定7個水體占比在5%～30%的樣地（A，B，C，D，E，F，G）來代替整個研究區進行精度對比研究。

選擇2023年7月12日的Sentinel-1影像和云量情況較好的7月16日Sentinel-2影像分別使用VV-OTSU和MNDWI-OTSU方法進行水體范圍提取（7月12日至7月16日研究區內無降雨發生），并與相同時段的FRB圖像進行一一比對，由于篇幅原因，僅展示具有代表性的樣地A的比對效果，如圖3所示，紅圈揭示出兩種方法提取水體的差異。總體來講，兩種方法都取得了較好的效果，但是VV-OTSU方法相對于MNDWI-OTSU方法在細小水體較多的場景下仍有不足，主要表現在：

①部分零碎水體漏提、誤提;② 由于SAR圖像濾波處理，弱化了部分水體邊緣，將水體之間的一些空隙連接了起來，因此在此場景下使用VV極化提取的水體面積比光學提取的面積更大。

隨后，為了排除時間上的偶然性，選擇時間相近且時間段內無降雨發生的3組Sentinel-1/2影像對，每個影像對各自影像的成像時間如表2所列。

分別使用VV-OTSU和MNDWI-OTSU對7個樣地的水體范圍進行提取，并統計水體面積占比，共得到21組數據，對其進行擬合分析得到圖4中的結果。由圖4，可以看出兩組數據的水體提取結果具有很強的線性相關性（R2=0.9，RMSE=2.02%），證明在Sentinel-2數據缺失的多云天氣，利用VV-OTSU對Sentinel-1進行提取得到的水體范圍也足夠穩定、可靠，可在該方法的基礎上，拓展出后續模型。

3.2 基于Sentinel-1/2數據的洪水淹沒范圍快速提取模型

從3.1節的實驗結果來看，MNDWI-OTSU比VV-OTSU方法更優，因此本文構建了一種綜合利用Sentinel-1/2數據的洪水淹沒范圍快速提取模型來優化單獨使用Sentinel-1數據的提取精度，其流程如圖5所示。

該模型首先篩選出區域內用戶所需時間段的Sentinel-1/2影像；然后對Sentinel-1影像進行預處理（Lee濾波等），并采用VV-OTSU方法提取水體得到區域水體粗分布圖；接著根據時間間隔不超過3 d、間隔期間不發生降雨的時間配對原則篩選Sentinel-2影像與對應Sentinel-1影像配對，然后對配對的Sentinel-2影像進行預處理，對s2cloudless去云處理（云量分數閾值設為50）和MNDWI-OTSU提取水體，得到部分區域的水體精分布，再對水體粗分布圖進行掩膜得到初步結果。為了防止Sentinel-2影像缺失時，Sentinel-1影像山體陰影與水體相互混淆，采用SRTM數據計算坡度柵格，將坡度大于5°的水體像元去除。模型搭載于GEE平臺，實現數據快速獲取、影像批量調用，極大提高了模型運行效率。

量化驗證MNDWI-OTSU、VV-OTSU、本文模型的精度如下：在7個樣地內隨機采集500個樣本點，包含250個水體點和250個非水體點，然后根據FRB圖像目視解譯結果，得到3種方法的水體提取精度，如表3所列。云量較少的Sentinel-2影像配合MNDWI-OTSU方法取得了更好的水體提取精度（總體精度、Kappa系數均大于90%），而VV-OTSU方法提取精度略低一點，與3.1節中目視對比的結果保持一致；而本文模型相比VV-OTSU精度有一定提升（總體精度、Kappa分別提升1.2%和2.4%），在缺乏可用Sentinel-2影像的情況下，達到了優化VV-OTSU方法的效果。

4 洪水淹沒范圍時間序列提取應用

利用本文提出的基于Sentinel-1/2數據的洪水淹沒范圍快速提取模型，得到研究區2023年5～9月的地表水體分布時間序列，將每幅影像的水體面積進行統計，繪制水體面積變化如圖6所示。

2023年7月24日至8月5日水體面積陡然上升，由107.9 km2擴大到308.0 km2，水體大幅擴張了近2倍。根據保定日報報道：7月29日至8月1日保定市出現大暴雨和特大暴雨，總降水時長超過了80 h。全市平均降水量350 mm，最大雨強82.2 mm/h，250 mm以上（特大暴雨）降水面積占比轄區面積的90%。暴雨引起了嚴重的洪水災情，研究區境內白溝河、小清河等河道水位一路走高。水利部信息中心（http：∥xxzx.mwr.gov.cn/）數據顯示：8月4日08：00，涿州白溝河東茨村水文站水位達27.49 m，超警戒水位（25.50 m）1.99 m。圖7為涿州市洪水現場圖。

以2023年7月24日影像為災前影像，連續繪制后續時間點相比災前影像新增水體范圍圖。其中，離洪水發生時間最近、受災面積最廣的8月5日的研究區水體范圍如圖8所示。

淹沒面積較廣的行政區為高碑店市、雄縣以及安新縣，洪水淹沒的面積分別約為70.1，52.3，40.6 km2。將新增水體與Google地圖進行對比可以發現，洪水淹沒首先發生在漲水的河道周圍，如涿州市東北部的碼頭鎮，該鎮境內琉璃河、拒馬河蜿蜒流淌，受暴雨影響出現嚴重積水；其次是天然的封閉洼地（圖8中高碑店市中部以西、定興東部）；最后是行洪河道（圖8中雄縣中部、安新縣西部），由于地勢差等因素，洪水通過這些河道被引至附近的湖泊和河流中。

2023年8月5日至9月22日，各縣相比7月24日新增水體面積統計如表4所列。

退洪方面，至2023年8月17日，定興縣、涿州市、蓉城縣、高碑店市新增水體面積僅剩約2.2，2.8，0.8，5.6 km2，相比上期分別縮減87.7%，82.1%，77.8%，92.0%；至9月22日，上述四縣水體面積基本恢復到災前水平。雄縣境內主要的積水區域為地勢較低的行洪道，而安新縣內則有白洋淀蓄滯洪水，因此這兩縣的水體消退相比上述四縣更為緩慢。

5 結 論

本研究以Sentinel-1/2數據為數據源，使用NDWI-OTSU、MNDWI-OTSU、VV-OTSU方法分別提取兩類數據的地表水體范圍。實驗結果表明：提取云量較少（<10%）的Sentinel-2影像水體時，MNDWI-OTSU比NDWI-OTSU閾值模型表達更為精細、簡潔，精度達到95.6%；VV-OTSU閾值模型提取Sentinel-1影像水體相比MNDWI-OTSU閾值模型存在零碎水體誤提、邊緣模糊的問題，但總體精度仍達到93.8%。兩種閾值模型在時間上提取的水體范圍擬合度較高（R2=0.9），不隨時間改變發生較大變化，具有足夠的可信度。本文在VV-OTSU方法的基礎上引入Sentinel-2數據和MNDWI-OTSU方法，構建了一種洪水淹沒范圍快速提取模型，并以保定市東部6縣為例，對該區域2023年5～9月時序地表水體進行了提取與分析。結果表明，該模型精度達到95%，不受限于云霧且比VV-OTSU方法的精度更高。該模型可通過編碼搭載于GEE平臺，從而能快速、準確、低成本地實現對一定時間范圍內某一區域的地表水體空間演變圖連續輸出。

參考文獻：

［1］ 《中國水旱災害防御公報》編寫組.《中國水旱災害防御公報2022》概要［J］.中國防汛抗旱，2023，33（10）：78-82.

［2］ 胡榮明，姚曉宙，李朋飛，等.Landast-8影像的水體指數構建及穩定性研究［J］.測繪科學，2022，47（5）：150-155.

［3］ 吳慶雙，汪明秀，申茜，等.Sentinel-2遙感圖像的細小水體提取［J］.遙感學報，2022，26（4）：781-794.

［4］ 趙文舉，李聰聰，馬宏，等.基于Sentinel-2超分辨率影像的干旱區水體提取方法［J］.農業機械學報，2023，54（10）：316-328.

［5］ 劉浩，周萬蓬，張宇佳，等.基于Landsat影像的1999～2019年鄱陽湖面積動態監測［J］.東華理工大學學報（自然科學版），2023，46（1）：68-76.

［6］ 趙蓓蓓，黃海峰，鄧永煌，等.基于Sentinel-1A的三峽庫區范家坪滑坡InSAR監測分析［J］.人民長江，2022，53（10）：103-107.

［7］ 賈詩超，薛東劍，李成繞，等.基于Sentinel-1數據的水體信息提取方法研究［J］.人民長江，2019，50（2）：213-217.

［8］ AMITRANO D，DI MARTINO G，IODICE A，et al.Unsupervised rapid flood mapping using Sentinel-1 GRD SAR images［J］.IEEE Transactions on Geoscience and Remote Sensing，2018，56（6）：3290-3299.

［9］ BINH P D，CATHERINE P，FILIPE A.Surface water monitoring within Cambodia and the Vietnamese Mekong Delta over a year，with Sentinel-1 SAR observations［J］.Water，2017，9（6）：366.

［10］劉小燕，崔耀平，史志方，等.GEE平臺下多源遙感影像對洪災的監測［J］.遙感學報，2023，27（9）：2179-2190.

［11］李健鋒.基于多源遙感影像和Google Earth Engine的斯里蘭卡地表水體及洪水信息快速提取研究［D］.西安：長安大學，2020.

［12］郝金虎，韋瑋，李勝男，等.基于GEE平臺的京津冀長時序水體時空格局及其影響因素［J］.生態環境學報，2023，32（3）：556-566.

［13］LEE J S.Speckle filtering of synthetic aperture radar images ：a review［J］.Remote Sensing Reviews，1994，8（4）：313-340.

［14］孫亞勇，黃詩峰，李紀人，等.Sentinel-1A SAR數據在緬甸伊洛瓦底江下游區洪水監測中的應用［J］.遙感技術與應用，2017，32（2）：282-288.

［15］OTSU N.A threshold selection method from Gray-Level Histograms［J］.IEEE Transactions on Systems Man & Cybernetics，2007，9（1）：62-66.

［16］陳賽楠，蔣彌.Sentinel-1 SAR在洪水范圍提取與極化分析中的應用研究［J］.地球信息科學學報，2021，23（6）：1063-1070.

［17］MCFEETERS S K.The use of the normalized difference water index（NDWI） in the delineation of open water features［J］.International Journal of Remote Sensing，1996，17（7）：1425-1432.

［18］XU H.Modification of normalised difference water index（NDWI） to enhance open water features in remotely sensed imagery［J］.International Journal of Remote Sensing，2006，27（14）：3025-3033.

［19］KHUONG H T，MASSIMO M，LI J.Surface water mapping and flood monitoring in the Mekong Delta using Sentinel-1 SAR time series and Otsu threshold［J］.Remote Sensing，2022，14（22）：5721.

（編輯：謝玲嫻）

Research on extraction model of flood inundation range based on Sentinel-1/2 data

DENG Qirui1，2，ZHANG Ying1，2，LIU Jia1，2，QIAO Qinghua1，2，ZHAI Liang1，2

（1.Chinese Academy of Surveying and Mapping，Beijing 100036，China； 2.Key Laboratory of Surveying and Mapping Science and Geospatial Information Technology of Ministry of Natural Resources，Beijing 100036，China）

Abstract：

Remote sensing is an important means of monitoring the extent of flood inundation and understanding the evolution of flood disasters.However，optical images often have many deficiencies during floods，and all-weather SAR images have slightly lower accuracy in extracting water bodies.A flood inundation range extraction model based on Sentinel-2 optical images and Sentinel-1 radar image data was constructed to extract the flood inundation range quickly and accurately.An adaptive threshold segmentation algorithm，the Otsu algorithm，was used to extract the water body range of two types of data and the proposed model，and an application analysis was conducted using Baoding City，Hebei Province as an example.The results showed that Sentinel-2 images with less cloud cover had the best water extraction effect，with an overall accuracy （OA） of 95.6%.After introducing some available Sentinel-2 data，the OA of the constructed model reached 95.0%，OA and Kappa coefficient were increased by 1.2% and 2.4% respectively compared to using Sentinel-1 data alone.This model is installed on the Google Earth Engine platform and can achieve fast，accurate，and low-cost continuous output of the spatial range of surface water bodies.Clouds and mist do not limit it and have higher extraction accuracy than Sentinel-1 images alone.In the case of severe cloud coverage leading to a lack of Sentinel-2 data，this model can be used as an alternative method for extracting flood inundation areas.

Key words：

flood inundation range； Sentinel-1； Sentinel-2； adaptive threshold segmentation algorithm； Google Earth Engine； Baoding City