劉向銅， 何佳君， 曹秋香， 徐 看， 陳 揚(yáng)

(1.東華理工大學(xué) 測繪工程學(xué)院，江西 南昌 330013；2.東華理工大學(xué) 地球科學(xué)學(xué)院，江西 南昌 330013)

洪災(zāi)是暴雨、風(fēng)暴潮流等自然因素引起的江河湖泊水量迅速增加，或者水位迅猛上漲的一種自然現(xiàn)象。洪水區(qū)域是指暫時性被水體淹沒的陸地區(qū)域(Patel et al.,2013)。洪災(zāi)會給自然環(huán)境和社會經(jīng)濟(jì)帶來巨大破壞，為減輕洪災(zāi)危害，快速準(zhǔn)確地獲取洪災(zāi)淹沒范圍變得尤為重要。水體信息獲取是洪災(zāi)監(jiān)測的基礎(chǔ)，遙感技術(shù)的發(fā)展和應(yīng)用為洪災(zāi)監(jiān)測提供了一種新的選擇(Elkhrachy, 2015)。由于洪災(zāi)發(fā)生時常伴隨著極端天氣，使用可見光和近紅外波段作為傳感器的被動光學(xué)遙感衛(wèi)星易受云層影響，因此在洪災(zāi)發(fā)生期間的作業(yè)效率較低(Kuenzer et al.,2013)。以微波作為傳感器的主動合成孔徑雷達(dá)衛(wèi)星(SAR)，能夠提供不受光照和天氣情況影響的全天時全天候監(jiān)測影像，是理想的洪災(zāi)監(jiān)測數(shù)據(jù)源(O′Grady et al.,2014)。

區(qū)別于光學(xué)衛(wèi)星高光譜特征，雷達(dá)影像主要是根據(jù)地物的后向散射系數(shù)、紋理特征和地物輪廓等進(jìn)行地物識別分類(周晗，2018)。基于SAR數(shù)據(jù)的水體提取主要是依據(jù)水體后向散射系數(shù)較低的特點(diǎn)，通過閾值分割，劃分水體與背景的二分類方法。Martinis等(2015)總結(jié)了一些基于合成孔徑雷達(dá)后向散射系數(shù)及其衍生特征(差分干涉、相干性和極值等)的水體提取方法，證明了SAR數(shù)據(jù)在區(qū)域和全球范圍內(nèi)的水體識別可行性。然而單一極化方式的雷達(dá)后向散射信息所包含地物特征信息較少，難以在洪澇期獲取高精度的水體信息。賈詩超等(2019)組合Sentinel-1雙極化數(shù)據(jù)，提出了一種增強(qiáng)水體特征，同時有效抑制土壤和植被存在的雙極化水體指數(shù)(SDWI)，并以巢湖和鄱陽湖為例驗(yàn)證了該指數(shù)的有效性。郭欣等(2018)基于Sentinel-1A數(shù)據(jù)對比GMM模型雙峰法、Otsu算法和區(qū)域生長法三種閾值分割算法在SAR影像水體提取中的優(yōu)劣性。

SDWI指數(shù)憑借其較強(qiáng)的背景抑制能力，已成為目前SAR影像提取水體信息的主流手段之一，但該指數(shù)在閾值選取上具有較強(qiáng)的主觀性。本研究基于SDWI水體指數(shù)，采用Otsu方法選取最優(yōu)閾值，提取鄱陽湖地區(qū)水體信息，對2020年鄱陽湖特大洪災(zāi)進(jìn)行監(jiān)測，分析洪災(zāi)發(fā)生前后鄱陽湖水面面積變化，繪制洪災(zāi)淹沒區(qū)域圖，并對2016—2020年鄱陽湖地區(qū)歷年水體面積及洪災(zāi)淹沒區(qū)的時空變化進(jìn)行了回顧性分析。

1 研究區(qū)及數(shù)據(jù)概況

1.1 研究區(qū)概況

鄱陽湖是中國第一大淡水湖，位于江西省北部，是受長江、贛江、撫河、信江、饒河、修水等河水水量吞吐平衡而形成的過水性、吞吐型、季節(jié)性的湖泊(米振華等，2021；吳瓊等，2020)。區(qū)域內(nèi)以亞熱帶季風(fēng)氣候?yàn)橹鳎杲邓窟_(dá)1 600～1 900 mm，降水主要集中在每年的4—6月(黃學(xué)輝等，2018)，具有獨(dú)特的豐枯特征，在調(diào)節(jié)長江徑流、維護(hù)區(qū)域生態(tài)平衡、支持經(jīng)濟(jì)發(fā)展發(fā)揮著重要作用(沈泰等，2003)，鄱陽湖地理位置及影像覆蓋范圍如圖1所示。

1.2 數(shù)據(jù)源

本研究采用Sentinel衛(wèi)星數(shù)據(jù)，以6天為時間間隔選取2020年7月2日至8月1日影像以及2016—2020年5月至9月影像(https://search.afs.alaska.edu)，對鄱陽湖區(qū)域洪災(zāi)進(jìn)行監(jiān)測。輔助數(shù)據(jù)包括中國縣級行政區(qū)劃、30 m分辨率SRTM外部DEM數(shù)據(jù)和鄱陽湖地區(qū)國家氣象監(jiān)測站降雨量統(tǒng)計數(shù)據(jù)(http://data.cma.cn/)等。

2 研究方法

2.1 基于SDWI的水體提取模型

在雷達(dá)影像中，像元記錄了目標(biāo)地物的后向散射系數(shù)。后向散射系數(shù)是雷達(dá)系統(tǒng)參數(shù)及地表特征的函數(shù)(吳效勇等，2020)。對水體目標(biāo)而言，雷達(dá)所發(fā)射的微波信號以鏡面反射為主，后向散射系數(shù)較低，在影像中表現(xiàn)為暗色調(diào)。本研究基于SDWI雙極化水體指數(shù)，通過閾值分割法提取水體信息。

SDWI水體指數(shù)由賈詩超等(2019)受DNVI和NDWI方法的啟發(fā)，結(jié)合微波遙感在水體信息影像中的特點(diǎn)，以及Sentinel-1雙極化數(shù)據(jù)(VV和VH)所提出的新型SAR影像水體提取模型。該模型可以有效增強(qiáng)SAR影像中的水體信息并抑制土壤和植被信息。公式如下：

KSDWI=In(10×KVV×KVH)

(1)

式中，KSDWI為SDWI水體指數(shù)，KVV表示雙垂直同向極化雷達(dá)強(qiáng)度影像，KVH表示垂直水平異向極化雷達(dá)強(qiáng)度影像。

在目標(biāo)地物與背景地物有明顯差異的情況下，依靠閾值對影像進(jìn)行分割可以獲得較好的分類精度，即在圖像的灰度統(tǒng)計直方圖中，有兩個“峰值”和一個“谷值”的情況下分割效果較好(陳志國，2017)。計算SDWI指數(shù)，對計算前后影像進(jìn)行灰度統(tǒng)計并作圖(圖2)，發(fā)現(xiàn)SDWI指數(shù)的灰度統(tǒng)計圖出現(xiàn)了更為明顯的峰值和谷值，相較于單極化影像，采用組合了雙極化數(shù)據(jù)的SDWI水體指數(shù)能夠獲取更高精度的水體提取結(jié)果。

2.2 最大類間方差分割

Otsu算法又稱最大類間方差法，是日本學(xué)者Otsu(1979)提出的一種全局最優(yōu)閾值確定方法。Otsu算法的基本思想是利用圖像的灰度特征，將圖像劃分為目標(biāo)與背景兩個部分，以目標(biāo)和背景的最大方差來動態(tài)確定圖像的分割閾值。采用Otsu算法確定閾值可以避免人工經(jīng)驗(yàn)閾值無法顧及全局的限制，是圖像分割中閾值選取的最佳算法，計算簡單，不受圖像亮度和對比度的影響。Otsu最佳閾值獲取公式如下：

Γ(τ)=σ2(τ)τ20(τ)+σ21(τ)

(2)

式中，τ表示臨時閾值的像元值；σ2(τ)表示水體和非水體的兩類類間方差；σ20(τ)表示水體類的平均方差；σ21(τ)表示非水體類的平均方差。Γ(τ)取最大值對應(yīng)的像元值τ即為最優(yōu)閾值T：

T=arg[maxΓ(τ)]

(3)

通過最優(yōu)閾值T，將SDWI指數(shù)值低于T的像元劃分為水體，高于T的像元劃分為非水體。

2.3 水體信息提取與精度驗(yàn)證

對圖像進(jìn)行配準(zhǔn)、多視、地理編碼等預(yù)處理流程，獲取15 m分辨率地圖坐標(biāo)系下的VV、VH強(qiáng)度影像，采用Lee濾波消除平坦地區(qū)的斑點(diǎn)噪聲增強(qiáng)影像邊緣信息(江勇等，2009)。以濾波后的VV和VH兩種極化強(qiáng)度影像為基礎(chǔ)計算SDWI雙極化水體指數(shù)，通過Otsu閾值分割進(jìn)行水體信息提取。

為驗(yàn)證水體提取精度，獲取了云量較少的2020年8月14日Sentinel-2以及8月21日Landsat-8光學(xué)影像數(shù)據(jù)。以光學(xué)影像作為參考影像驗(yàn)證8月13日和8月19日水體提取精度，其結(jié)果如表1所示。

表1 水體提取精度

3 鄱陽湖洪災(zāi)監(jiān)測

3.1 2020年鄱陽湖洪災(zāi)監(jiān)測

利用2020年7月份3景Sentinel1-A以及3景Sentinel1-B影像數(shù)據(jù)，結(jié)合九江站水位進(jìn)行鄱陽湖水體面積及水位統(tǒng)計(圖3)，提取水體信息，并計算其面積。根據(jù)中國地面氣候資料數(shù)據(jù)集(http://data.cma.cn/)獲取鄱陽湖區(qū)域九江站、廬山站、南昌站和鄱陽站4個國家基準(zhǔn)氣象站的降水量數(shù)據(jù)以及長江流域九江站水流量數(shù)據(jù)(圖4)，對鄱陽湖特大洪災(zāi)進(jìn)行分析。

圖4顯示自2020年7月2日以來，鄱陽湖地區(qū)降水量增加，僅7月10日單日降水量高達(dá)664.7 mm，降水量的激增導(dǎo)致鄱陽湖地區(qū)水量迅速累積。同時自2020年7月3日起長江流量也呈現(xiàn)快速增長趨勢，來自長江中上游的大量洪水匯入鄱陽湖，加劇了鄱陽湖水量累積，導(dǎo)致鄱陽湖九江水位站水位迅速攀升，于7月12日達(dá)到22.74 m的最高水位，鄱陽湖水體面積快速增加，于7月14日達(dá)到最大。2020年7月14日之后，鄱陽湖地區(qū)降水量明顯減少，但長江流量仍處于較高水平，致使鄱陽湖地區(qū)退洪速度緩慢，九江站水位以及鄱陽湖地區(qū)水面面積呈現(xiàn)緩慢下降趨勢。大量的降水以及長江中上游匯入的大量洪水是導(dǎo)致鄱陽湖地區(qū)特大洪災(zāi)的直接原因。

為更為直觀地反映洪災(zāi)受災(zāi)情況，以Landsat-8多光譜影像為底圖，2020年7月2日水體作為災(zāi)前參考水體，7月14日水體作為災(zāi)中水體，結(jié)合2005—2006年鄱陽湖1∶25萬矢量邊界進(jìn)行鄱陽湖洪災(zāi)淹沒區(qū)提取，得到2020年鄱陽湖洪災(zāi)淹沒區(qū)域圖(圖5)。識別出A、B兩處大型淹沒區(qū)域，分別位于永修縣、鄱陽縣境內(nèi)，兩處分別于2020年7月12日、7月11日發(fā)生潰堤，為此次洪災(zāi)災(zāi)情中心。

3.2 2016—2020年鄱陽湖洪災(zāi)回顧

鄱陽湖區(qū)域?qū)儆诤闉?zāi)高發(fā)地帶，每年6至7月屬于洪災(zāi)高發(fā)月，收集2016年至2020年的5月至9月的存檔Sentinel-1衛(wèi)星數(shù)據(jù)，對鄱陽湖區(qū)域水體面積變化及洪災(zāi)進(jìn)行回顧分析。

對存檔Sentinel-1影像進(jìn)行水體提取并進(jìn)行面積統(tǒng)計(圖6)，同時收集2016年至2020年的6至7月降水量數(shù)據(jù)(圖7)。在查閱2016年至2020年最大洪災(zāi)發(fā)生時間的基礎(chǔ)上，選擇鄰近最大洪災(zāi)發(fā)生日影像(20160716、20170706、20180713、20190714、20200714)作為洪災(zāi)影像，按照上文矢量疊加方法，繪制2016—2020年各年洪災(zāi)淹沒區(qū)域圖(圖8)，分析淹沒區(qū)域時空變化。

圖6顯示，鄱陽湖地區(qū)水體面積呈現(xiàn)明顯的季節(jié)性規(guī)律，在月份變化上呈現(xiàn)5至6月增加，8至9月減少的總體趨勢，最大水體面積均出現(xiàn)在每年7月份。2016年和2020年水體面積總體上處于較高水平，為豐水年份；2018年水體面積處于較低水平，屬于干旱年份。通過圖7可以得出，鄱陽湖地區(qū)水體在面積上與區(qū)域內(nèi)降水量成正比。2018年6至7月累積降水量僅為1 377.4 mm，為近5年最低水平，同時其6月水體面積與災(zāi)中水體面積均處于較低水平，2020年6至7月累計降水量高達(dá)3 428.8 mm，導(dǎo)致其出現(xiàn)了5 354.91 km2的最大災(zāi)中面積。受災(zāi)面積即水體變化面積，2018年災(zāi)前水體與災(zāi)中水體均處于較低水平，但變化面積高達(dá)627.64 km2，變化面積僅次于2020年的1 174.19 km2。結(jié)合圖8發(fā)現(xiàn)2018年水體面積雖有較大的變化，但由于水體面積整體處于較低水平，其變化區(qū)域絕大部分處于鄱陽湖邊界之內(nèi)，屬于干旱年份的自然變化，因此造成的陸地淹沒面積較小。從淹沒區(qū)域的空間變化上看，B區(qū)域位于鄱陽湖邊界內(nèi)部，最易受鄱陽湖水位影響，屬常駐淹沒區(qū)，A、C兩區(qū)域分別位于永修縣和鄱陽縣，為新增淹沒區(qū)域。

4 結(jié)論

(1)利用多極化的SAR影像后向散射信息，可有效地獲取水體信息，能在洪災(zāi)發(fā)生期間快速獲取洪水淹沒面積及范圍。

(2)2020年鄱陽湖洪災(zāi)最大災(zāi)情發(fā)生于7月14日前后，最大受災(zāi)面積達(dá)1 174.19 km2，受災(zāi)嚴(yán)重區(qū)域位于永修縣與鄱陽縣境內(nèi)，區(qū)域內(nèi)超高的降水量以及長江中上游源源不斷的洪水匯入是導(dǎo)致洪災(zāi)發(fā)生的主要原因。

(3)2016—2020年鄱陽湖地區(qū)水體在面積上與區(qū)域內(nèi)降水量高度相關(guān)，在空間分布上存在常駐淹沒區(qū)，2018年淹沒區(qū)均位于鄱陽湖內(nèi)部，2020年由于潰堤事件新增了兩處淹沒區(qū)域。