楊雨晴，王冬梅，徐 佳

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098；2.江蘇省水利科學(xué)研究院，江蘇 南京 210017）

洪澇災(zāi)害具有持續(xù)時間短、危害性大等特征，一直危害著人類社會的發(fā)展。因此，及時防御洪澇災(zāi)害，減少災(zāi)害損失，保障人類安全，顯得至關(guān)重要[1]?；谛l(wèi)星遙感技術(shù)，人們可實現(xiàn)快速、準(zhǔn)確的洪澇淹沒范圍提取，直觀地顯示其空間分布及其動態(tài)變化和發(fā)展規(guī)律，從而使得其在洪澇災(zāi)害監(jiān)測中發(fā)揮著越來越重要的作用[2]。經(jīng)過近幾十年的探索和實踐，基于遙感技術(shù)的洪災(zāi)應(yīng)用，逐步形成了貫穿于災(zāi)前、災(zāi)中和災(zāi)后全過程的洪澇災(zāi)害監(jiān)測。水體信息提取是洪水監(jiān)測的出發(fā)點。傳統(tǒng)的監(jiān)測方法主要通過采樣調(diào)查，耗時耗力，花費巨大，遙感技術(shù)的發(fā)展和應(yīng)用則提供了一種新的選擇。光學(xué)傳感器雖然有著較好的洪水監(jiān)測潛力，但時間分辨率和空間分辨率往往不能兼顧，且易受云層和天氣影響，而合成孔徑雷達(dá)（SAR）受多云與陰雨天氣影響小，可實現(xiàn)全天候、全天時的大面積觀測[3]，已成為洪澇災(zāi)害監(jiān)測的重要技術(shù)手段。一直以來，我國在民用SAR遙感衛(wèi)星方面與其他發(fā)達(dá)國家存在較大差距。2016-08-10成功發(fā)射的高分三號衛(wèi)星是我國首顆C頻段多極化高分辨率SAR衛(wèi)星，其具有高分辨率、大幅寬、多成像模式等特點，緩解了我國SAR數(shù)據(jù)的短缺問題，為洪澇災(zāi)害研究提供了重要的數(shù)據(jù)源。

對于洪澇災(zāi)害遙感監(jiān)測而言，從遙感影像中快速準(zhǔn)確地提取洪澇水體信息至關(guān)重要?；赟AR影像的水體信息提取方法[4]主要有閾值分割法、機(jī)器學(xué)習(xí)法、曲線演化法等。由于洪澇災(zāi)害監(jiān)測屬于應(yīng)急監(jiān)測，更側(cè)重于較高的提取速度，因此閾值分割法在SAR影像洪澇水體提取中應(yīng)用最為廣泛，如李景剛等[5]采用改進(jìn)的OTSU法進(jìn)行水體提取研究；郭欣等[6]采用雙峰法和OTSU法對Sentinel-1A SAR影像進(jìn)行洪澇淹沒監(jiān)測等。然而在大場景SAR影像中水體目標(biāo)相對于整幅影像占比一般較小，全局閾值法往往提取效果不佳，特別是當(dāng)研究區(qū)內(nèi)存在地形起伏的山區(qū)時，山體陰影亦會對水體提取造成干擾。因此，不少學(xué)者提出了改進(jìn)算法，如安成錦等[7]利用多閾值分割方法進(jìn)行了SAR圖像水域分割；Junhua等[8]結(jié)合紋理信息，采用非監(jiān)督K-means聚類獲取初始水體范圍；龐科臣等[9]利用改進(jìn)的OTSU法對DEM數(shù)據(jù)進(jìn)行分割來去除誤提為水體的山體陰影。本文在以上研究的基礎(chǔ)上，針對高分三號SAR影像特點，提出了一種DEM輔助下基于局部自適應(yīng)閾值的SAR影像洪澇水體提取方法，利用GF-3精細(xì)條帶2成像模式（FSII）SAR影像進(jìn)行洪澇水體提取實驗，并實現(xiàn)了對湖南省東北區(qū)域洪水淹沒范圍的快速提取。

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

實驗區(qū)位于湖南省東北部，如圖1b所示。該實驗區(qū)位于益陽市和長沙市交接處，其中，寧鄉(xiāng)市所占面積最大。寧鄉(xiāng)市，原名寧鄉(xiāng)縣，2017年撤縣建市后，成為湖南省長沙市管轄下的縣級市。該實驗區(qū)介于112°07′E～112°44′E， 28°03′N～28°31′N 之間，屬中亞熱帶向北亞熱帶過渡的大陸性季風(fēng)濕潤氣候，四季分明，溫度適宜，雨水充足，寒冷期短，炎熱期長。每日平均氣溫16.8℃，每年平均降水量1 358.3 mm。內(nèi)部有溈水河及其他支流穿過，有較多山地分布，山間有較小的湖泊和水洼。境內(nèi)地貌以山地、丘崗、平原為主，平均海拔高度為153 m。

本文實驗數(shù)據(jù)采用一景GF-3 L1影像、SRTM DEM數(shù)據(jù)和兩景Sentinel-2A數(shù)據(jù)。高分三號數(shù)據(jù)的獲取時間為2017-07-09，成像模式為精細(xì)條帶2模式，極化方式為HH和HV雙極化，距離向分辨率為1.1 m，方位向分辨率是5.8 m，入射角為27.8°，如圖1c所示，該RGB圖像由R：HV極化；G：HH極化；B：HV極化組成。Sentinel-2A的獲取時間為2017-05和2017-07，產(chǎn)品級別為Level-1C，是經(jīng)正射校正和亞像元級幾何精校正后的大氣表觀反射率產(chǎn)品，其中5月數(shù)據(jù)為災(zāi)前數(shù)據(jù)，7月為驗證數(shù)據(jù)。

圖1 研究區(qū)示意圖

2 研究方法

為解決大場景SAR影像中水體目標(biāo)占比較小問題，本文首先將原始影像劃分為K×K個局部區(qū)域，然后進(jìn)一步將每個局部區(qū)域分為N×N個子塊，通過統(tǒng)計子塊中初始水體和非水體的比例來篩選子塊，從而提高局部區(qū)域中水體占比。同時，考慮到山體陰影會造成水體誤提和輻射失真，在閾值分割前，利用DEM和SAR軌道參數(shù)信息生成陰影掩膜文件從而去除SAR圖像的輻射失真區(qū)域。本文提出的適用于高分三號SAR影像的洪澇水體提取算法流程如圖2所示，主要包含5個步驟：①影像預(yù)處理：對GF-3數(shù)據(jù)進(jìn)行輻射定標(biāo)、多視、濾波、地理編碼等處理；②山體陰影掩膜文件制作：將DEM和SAR軌道參數(shù)信息相結(jié)合，生成陰影掩膜文件；③初始水體分布獲?。豪蒙襟w陰影掩膜文件去除SAR圖像的輻射失真區(qū)域，然后將K-means聚類方法應(yīng)用于陰影掩膜后的SAR圖像，得到初始的水體分布圖；④洪澇水體自適應(yīng)閾值提?。簩τ跋穹謪^(qū)分塊處理，利用KI方法對優(yōu)化后的局部區(qū)域進(jìn)行閾值分割，提取最終的洪澇災(zāi)區(qū)水體范圍。⑤洪澇災(zāi)害監(jiān)測：以Sentinel-2A數(shù)據(jù)為參考，提取災(zāi)前水體分布圖，并將本文方法提取的災(zāi)后水體分布圖與其進(jìn)行差分處理，得到洪水分布圖。

圖2 水體信息提取流程

2.1 山體陰影檢測

由于SAR成像系統(tǒng)為側(cè)視成像，SAR圖像存在固有的幾何失真，在地形起伏較大的山區(qū)更容易產(chǎn)生。受山體的影響，背向坡將無法獲得雷達(dá)波束，在影像上就形成陰影，陰影區(qū)的大小取決于俯角和斜坡坡度兩個參數(shù)。雷達(dá)陰影在圖像相應(yīng)位置呈現(xiàn)暗色，即使圖像經(jīng)過正射校正，雷達(dá)陰影區(qū)的信息也無法得到補(bǔ)償，其后向散射系數(shù)和水體的后向散射系數(shù)仍然接近，使得圖像分割無法有效區(qū)別陰影和水體信息。因此，為了消除山體陰影的影響，本文首先采用DEM數(shù)據(jù)模擬SAR圖像，然后根據(jù)陰影區(qū)的低亮度特性，從模擬SAR圖像中檢測山體陰影并將其掩膜。

由雷達(dá)成像原理可知，DEM模擬SAR圖最主要的影響因素是本地入射角的大小，其計算公式[10]為：

式中，θ為雷達(dá)入射角；θA、θC為分辨單元在距離向和方位向的坡度角。

在模擬雷達(dá)圖像時，首先根據(jù)分辨率單元的大小將地面分成格網(wǎng)，對地面不同部分的雷達(dá)波束入射角進(jìn)行計算。由于雷達(dá)圖像為斜距成像，需考慮地距和斜距的關(guān)系，即：

式中，R為斜距；G為地距；θ為雷達(dá)波束入射角。

根據(jù)地物的散射特性和成像原理，按雷達(dá)方程和灰度方程，計算每一塊地面單元的散射回波強(qiáng)度[11]。雷達(dá)方程式如式（3）所示：

式中，pr為接收功率；pt為發(fā)射功率；G為天線增益；λ為波長；R為天線和目標(biāo)直接的距離；σ0為后向散射系數(shù)；ΔA為分辨單元。

2.2 KI閾值分割算法原理

陸地水域以鏡面散射為主，比植被、城鎮(zhèn)等非水體表面的后向散射弱，因此可通過閾值分割方法提取水體信息，當(dāng)圖像的后向散射強(qiáng)度x小于閾值時定義為水體，大于閾值定義為非水體。本文采用KI算法進(jìn)行圖像分割，該方法由Kittler和Illingworth于1986年提出，以圖像的概率分布為特征，基于最小錯分率和貝葉斯理論來進(jìn)行閾值選擇[12]。

在代價函數(shù)的輔助下，當(dāng)滿足總體分類代價最小時即可獲得具有最小錯分率的結(jié)果。

在KI算法中，假設(shè)直方圖為h(x)，則通過以下目標(biāo)函數(shù)尋找最優(yōu)閾值：

(2) 專用耐磨軋制鋼板NM360系列。NM360是在Q345基礎(chǔ)上通過調(diào)整合金元素后生產(chǎn)的新鋼種，主要應(yīng)用于礦山機(jī)械、煤礦機(jī)械、環(huán)保機(jī)械、工程機(jī)械等，也常用作為屈服強(qiáng)度不小于700 MPa高強(qiáng)度結(jié)構(gòu)鋼使用。主要是為需要耐磨的場合或部位提供保護(hù)。該襯板的屈服強(qiáng)度大于800 MPa，抗拉強(qiáng)度大于1 000 MPa，布氏硬度為360。其主要化學(xué)成分見表1。

其中

這種利用控制最小化分類誤差而確定閾值的方法可以通過最小化函數(shù)而獲得[13]，即

式中，T*為最優(yōu)閾值，即分割閾值。

2.3 基于局部自適應(yīng)閾值的水體信息提取

考慮到大場景SAR影像中水體目標(biāo)占比較小，全局閾值法難以獲得滿意的提取效果。本文通過對原始影像分區(qū)分塊，在局部區(qū)域中利用KI算法計算自適應(yīng)的局部閾值，實現(xiàn)影像的水體和非水體精確分割。具體步驟如下：

1）提取局部區(qū)域。將影像分為K×K個局部區(qū)域，K值大小可根據(jù)影像大小適當(dāng)選取，本文K取6。

2）獲取新的局部區(qū)域。當(dāng)水體占影像比例較小或較大時，KI方法計算的閾值都不準(zhǔn)確，因此為了提高閾值計算的準(zhǔn)確性，進(jìn)一步將每個局部區(qū)域分為N×N個子塊，根據(jù)初始水體分布圖對子塊進(jìn)行篩選，當(dāng)該子塊滿足水體占圖像一定比例（0.5±0.15）時，將這些子塊中的所有像元作為新的局部區(qū)域。

3）計算局部區(qū)域閾值。根據(jù)上一步中提取的所有像元值進(jìn)行判斷：當(dāng)局部區(qū)域內(nèi)幾乎完全為水體或非水體時，以全局閾值作為該窗口的閾值；當(dāng)局部區(qū)域有水體與非水體時，將提取的所有窗口的像元值（即新的局部區(qū)域）通過KI算法計算出水體和非水體的分割閾值。

3 實驗結(jié)果分析

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

本文使用PIE軟件進(jìn)行GF-3 SAR數(shù)據(jù)預(yù)處理工作，具體包括影像數(shù)據(jù)導(dǎo)入、輻射定標(biāo)、復(fù)數(shù)據(jù)轉(zhuǎn)換、多視處理、濾波、地理編碼、轉(zhuǎn)dB影像、研究區(qū)裁剪等。其中多視處理選擇距離向和方位向的視數(shù)為5和2；濾波過程中選擇7×7窗口的EnLee濾波器進(jìn)行濾波處理。最終獲得空間分辨率為15 m，WGS-84坐標(biāo)系下的后向散射系數(shù)影像，如圖3所示。

圖3 GF-3 SAR影像后向散射系數(shù)影像

3.2 山體陰影掩膜

為了便于后續(xù)實驗的順利進(jìn)行，本研究結(jié)合DEM和SAR軌道參數(shù)生成了SAR模擬圖，如圖4所示，可以發(fā)現(xiàn)在SAR模擬圖中，陰影區(qū)域表現(xiàn)為黑色，后向散射系數(shù)相對較低，因此通過分割閾值的方法可以提取山體陰影區(qū)域，如圖5所示，并制作山體陰影掩膜文件，去除山體陰影對水體提取精度的影響。

圖4 SAR模擬圖

圖5 山體陰影提取

3.3 自適應(yīng)閾值實驗結(jié)果分析

為了解決全局閾值難以準(zhǔn)確提取水體信息的問題，本研究以HV極化為例進(jìn)行自適應(yīng)閾值實驗研究。首先采用非監(jiān)督分類K-means聚類方法對山體陰影掩膜后的影像進(jìn)行分類（K=15），并選擇聚類1和聚類2組合為水體類，其他的為非水體類，得到初始的水體范圍，如圖6所示。從圖6可以發(fā)現(xiàn)該方法得到的水體分布圖存在大量噪聲，有部分裸地、土壤等地物的影響，但是由于山體陰影掩膜的加入，山體陰影的影響被有效消除。

圖6 初始水體分布圖

接下來將山體陰影掩膜后的影像和初始水體分類圖分為N×N個局部區(qū)域（N的值適影像質(zhì)量和大小而定，此實驗N=6），在初始的水體范圍基礎(chǔ)上，用10×10的窗口對每個局部區(qū)域進(jìn)行計算，統(tǒng)計像元值滿足水體占圖像一定比例（0.5±0.15）的所有窗口，將這些窗口中的所有像元作為新的局部區(qū)域的水體和非水體，并用KI法進(jìn)行閾值分割。若局部區(qū)域內(nèi)完全為水體或非水體時，以KI全局閾值作為該窗口的閾值；若局部區(qū)域有水體與非水體時，將提取的所有窗口的像元值（即新的局部區(qū)域）通過KI算法得到水體和非水體的分割閾值。然后對原始影像進(jìn)行閾值分割，得到最終的水體分布圖，如圖7所示。圖8為KI方法對山體陰影掩膜后的影像進(jìn)行全局閾值的結(jié)果，可以發(fā)現(xiàn)全局統(tǒng)一閾值結(jié)果存在部分噪聲，且有裸地等地物的誤提和少量河流出現(xiàn)斷層現(xiàn)象。

圖8 全局閾值分割結(jié)果

從圖7、8可以看出，與全局閾值分割結(jié)果相比，自適應(yīng)閾值分割結(jié)果噪聲有所降低。區(qū)域1為溈水下游區(qū)域，如圖9a和10a所示。使用全局統(tǒng)一閾值提取水體時，閾值偏小，河流呈現(xiàn)斷斷續(xù)續(xù)狀態(tài)，如圖9a紅色框中所示，出現(xiàn)了漏提現(xiàn)象。而自適應(yīng)閾值的水體提取過程中，該現(xiàn)象有所改善，閾值相對較大；區(qū)域2為溈水中下游河流，如圖9b和10b，全局統(tǒng)一閾值法在水體提取過程中，閾值偏大，由于裸地等的后向散射系數(shù)與水體的后向散射系數(shù)相近，則會被分為同一類，誤提現(xiàn)象明顯，如圖9b紅色框中所示。而自適應(yīng)閾值法可以有效降低其他易混淆地物的干擾；區(qū)域3為七家灣區(qū)域如圖9c和10c，該區(qū)域水體內(nèi)部被兩條大壩所分隔，從全局閾值分割的結(jié)果圖可以看出，左邊的大壩出現(xiàn)斷斷續(xù)續(xù)現(xiàn)象，有部分被錯分為水體，右邊的大壩缺失比較嚴(yán)重，大部分被錯分為水體，而在自適應(yīng)閾值水體提取中，左邊的大壩則被完整的分割為非水體部分，右邊的大壩的分割結(jié)果也有所改善。

圖7 自適應(yīng)閾值分割結(jié)果

圖9 全局閾值水體提取局部放大結(jié)果

圖10 自適應(yīng)閾值水體提取局部放大結(jié)果

3.4 水體提取精度驗證

本文以災(zāi)后7月該地區(qū)的Sentinel-2數(shù)據(jù)為參照數(shù)據(jù)，將GF-3 SAR數(shù)據(jù)水體提取結(jié)果與Sentinel-2數(shù)據(jù)的水體提取結(jié)果進(jìn)行比較，采用查全率、查準(zhǔn)率和虛警率3個指標(biāo)對本文提出的自適應(yīng)閾值水體提取方法和全局單一閾值水體提取方法進(jìn)行優(yōu)劣性評價，查全率越高，提取結(jié)果可靠性越高，查準(zhǔn)率越高，監(jiān)測準(zhǔn)確性越高。隨機(jī)選取局部區(qū)域進(jìn)行精度評定，統(tǒng)計結(jié)果如表1所示。

表1 水體提取精度對比

從表1可以看出，在查全率相差不到4%的情況下，從監(jiān)測準(zhǔn)確性出發(fā)，查準(zhǔn)率從82.49%提升到90.46%，上升了7.97%。綜合考慮查全率和查準(zhǔn)率的指數(shù)F-Measure，局部自適應(yīng)閾值方法要優(yōu)于全局統(tǒng)一閾值方法。

3.5 洪澇災(zāi)害監(jiān)測實驗結(jié)果分析

該區(qū)域在2017-06-10～2017-07-03時間段發(fā)生連續(xù)降雨現(xiàn)象。本研究對2017-05-18的Sentinel-2A影像和2017-07-09的影像分別進(jìn)行了水體提取，獲得災(zāi)前和災(zāi)后水體，再將災(zāi)前和災(zāi)后水體疊加，獲得最終的洪水分布圖。比較該區(qū)域的災(zāi)前災(zāi)后的水體疊加圖，可以發(fā)現(xiàn)團(tuán)頭湖附近有明顯的水體變化情況，如圖11所示。

圖11 團(tuán)頭湖洪水分布圖

4 結(jié) 語

本研究基于高分三號數(shù)據(jù)，提出了一種有效的地形數(shù)據(jù)支持下的自適應(yīng)閾值洪水水體提取方法。本文以湖南省東北部地區(qū)為研究對象展開實驗，主要結(jié)論如下：

1）該研究區(qū)水體面積占整幅影像的比例較小，KI法計算的閾值偏大。利用自適應(yīng)閾值水體提取的方法可以有效獲取水體信息，去除部分裸地、土壤等地物的影像，并保留了水體的細(xì)節(jié)信息。

2）在SAR影像中陰影信息必不可少，且會對水體的提取造成一定的干擾。本文加入了地形數(shù)據(jù)，去除了輻射失真區(qū)域，有助于獲得更準(zhǔn)確的水體-非水體信息。

3）從該實驗區(qū)的洪澇災(zāi)害監(jiān)測的實驗結(jié)果來看，該方法對快速、準(zhǔn)確的洪澇災(zāi)害監(jiān)測具有很好的適用性。

本文方法也有不足之處， 由于KI法計算的全局統(tǒng)一閾值偏大，因此保留了河流邊緣的細(xì)節(jié)信息，導(dǎo)致噪聲和相似地物的干擾也會增強(qiáng)，而本文的方法雖然大大減少了噪聲和相似地物的影響，但是有些地方的水體細(xì)節(jié)保留情況不如KI法計算的全局統(tǒng)一閾值方法，使得查全率低于KI法計算的全局統(tǒng)一閾值的查全率。