基于局部自適應閾值的GF-3 SAR影像洪澇水體提取

楊雨晴，王冬梅，徐 佳

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

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

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

1 研究區和數據介紹

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

本文實驗數據采用一景GF-3 L1影像、SRTM DEM數據和兩景Sentinel-2A數據。高分三號數據的獲取時間為2017-07-09，成像模式為精細條帶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，產品級別為Level-1C，是經正射校正和亞像元級幾何精校正后的大氣表觀反射率產品，其中5月數據為災前數據，7月為驗證數據。

圖1 研究區示意圖

2 研究方法

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

圖2 水體信息提取流程

2.1 山體陰影檢測

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

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

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

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

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

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

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

2.2 KI閾值分割算法原理

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

在代價函數的輔助下，當滿足總體分類代價最小時即可獲得具有最小錯分率的結果。

在KI算法中，假設直方圖為h(x)，則通過以下目標函數尋找最優閾值：

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

其中

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

式中，T*為最優閾值，即分割閾值。

2.3 基于局部自適應閾值的水體信息提取

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

1）提取局部區域。將影像分為K×K個局部區域，K值大小可根據影像大小適當選取，本文K取6。

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

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

3 實驗結果分析

3.1 數據預處理

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

圖3 GF-3 SAR影像后向散射系數影像

3.2 山體陰影掩膜

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

圖4 SAR模擬圖

圖5 山體陰影提取

3.3 自適應閾值實驗結果分析

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

圖6 初始水體分布圖

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

圖8 全局閾值分割結果

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

圖7 自適應閾值分割結果

圖9 全局閾值水體提取局部放大結果

圖10 自適應閾值水體提取局部放大結果

3.4 水體提取精度驗證

本文以災后7月該地區的Sentinel-2數據為參照數據，將GF-3 SAR數據水體提取結果與Sentinel-2數據的水體提取結果進行比較，采用查全率、查準率和虛警率3個指標對本文提出的自適應閾值水體提取方法和全局單一閾值水體提取方法進行優劣性評價，查全率越高，提取結果可靠性越高，查準率越高，監測準確性越高。隨機選取局部區域進行精度評定，統計結果如表1所示。

表1 水體提取精度對比

從表1可以看出，在查全率相差不到4%的情況下，從監測準確性出發，查準率從82.49%提升到90.46%，上升了7.97%。綜合考慮查全率和查準率的指數F-Measure，局部自適應閾值方法要優于全局統一閾值方法。

3.5 洪澇災害監測實驗結果分析

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

圖11 團頭湖洪水分布圖

4 結 語

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

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

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

3）從該實驗區的洪澇災害監測的實驗結果來看，該方法對快速、準確的洪澇災害監測具有很好的適用性。

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