改進Otsu與形態學相結合的水體信息提取

陳景玨，劉瑞，楊鑫，楊梅，楊遠陶

(1.成都理工大學 地球科學學院，成都 610059；2.成都理工大學 地球物理學院，成都 610059)

0 引言

星載雷達衛星具有全天時、全天候、穿透云層的能力，可以利用星載雷達的優勢獲取多云暴雨情況下的合成孔徑雷達(synthetic aperture radar，SAR)影像用于快速洪澇災害評估。隨著SAR數據空間分辨率的提高，出現了種類繁多的超小尺度圖像目標，這給圖像處理和分析帶來了更大的挑戰。一方面，由于相干斑點的存在，影響了SAR影像質量；另一方面，建筑物陰影和道路等混亂的目標可能會影響分割結果，因此SAR影像分割是一個公認的難題[1]。

SAR影像分割廣泛使用的方法有：最小誤差閾值、熵閾值和Otsu等。Otsu也被稱為大津法，于1979年由Otsu[2]提出。Otsu因其原理易理解且穩定有效，被廣泛用于圖像分割[3]。SAR圖像中水體和背景內的像元光譜測度各具有同質性，可視為圖像中只有水體和背景兩類地物[4]，而Otsu可以快速地獲取兩類目標的分割閾值，因此Otsu可以快速分割水陸邊界。基于此特點，Otsu被廣泛地應用于SAR圖像分割領域[5-8]。雖然Otsu穩定高效，但是在SAR影像中目標與背景的分布通常是不平衡的，這造成Otsu法得到的閾值會偏向方差較大的一類，偏離正確的分割閾值。

SAR衛星是側式成像，使得圖像中包含植物、建筑物和山體所產生的陰影，并且呈現與水體相似的暗色調，造成分割結果中存在一些偽水體[9]。SAR圖像相干斑的存在，也會使得分割結果存在雜散點[10]。采用閾值分割法得到的水體分布圖，陰影和雜散點是主要的影響對象，剔除它們的影響顯得尤為重要，因此對閾值法得到的水體分布圖進行優化是準確提取水體的重要技術環節。針對水體在分割結果中大面積連通的特點，利用形態學中的連通區域標定算法除去小面積的雜散點和偽水體，可以有效保留真實水體。但是在海拔較高區域，產生的山體陰影可能比部分真實水體面積要大。在該條件下，使用形態學中連通區標定算法刪除偽水體的同時，面積較小的真實水體也會被刪除。

綜上所述，為解決傳統Otsu閾值分割法會偏向方差較大一類的問題，本文提出目標權重改進的Otsu算法。加入數字高程模型(digital elevation model，DEM)進行形態學連通區標定，提出結合高程連通標定算法的形態學優化方法。利用高程連通標定算法刪除面積較大的山體陰影，并繼續對小面積偽水體進行刪除，從而得到精細的水體分布圖。

1 研究區域與數據處理

1.1 研究區概況

2020年6—8月長江干流先后發生五次編號洪水，我國27省3 789萬人受災，195.8萬人緊急轉移安置，直接造成的經濟損失高達700億元，隸屬于江西省九江市的彭澤縣地區內澇嚴重，多地村民受災。連續暴雨、大暴雨引起的山洪暴發或江水倒灌對該地區會造成嚴重的威脅[11]。

1.2 數據來源

Sentinel-1A衛星于2014年成功發射升空，重訪周期為12 d。本文獲取了Sentinel-1A兩景單視復數(single look complex，SLC)影像，詳細參數如表1所示。使用航天飛機雷達地形測繪使命(shuttle radar topography mission，SRTM)90 m分辨率DEM(https：//srtm.csi.cgiar.org/)。

表1 研究數據

1.3 數據預處理

Sentinel-1A影像預處理包括：條帶拼接生成SLC數據、多視處理、濾波、地理編碼和輻射校正。多視處理犧牲了空間分辨率以提高影像質量[12]。本文使用Frost濾波方法，可以減少斑點噪聲的影響，同時保留了邊緣信息[13]。地理編碼將SAR影像雷達坐標轉為地理坐標，輻射校正根據入射角信息計算得出雷達回波后向散射系數。

1.4 哨兵1雙極化水體指數計算

賈詩超等[14]提出了哨兵1雙極化水體指數(sentinel-1 dual-polarize water index，SDWI)。SDWI將Sentinel-1衛星兩種極化影像后向散射系數相乘并擴大十倍，再以自然對數作為函數式。將研究區Sentinel-1A 進行SDWI運算后得到的圖像(圖1)，經SDWI運算后水體從暗色調轉換成了亮色調。本文把Sentinel-1A影像使用SDWI計算得到的圖像進行水體提取研究，將ND值量化到256個灰度級。如圖2所示，VV和VH極化影像直方圖雙峰形狀基本一致，波谷與第一個波峰起伏較小；SDWI圖像的直方圖比VV和VH極化影像的直方圖雙峰現象更為明顯。 其中，T1是Otsu計算出地分割閾值；T2是改進Otsu計算出地分割閾值；T3是目視解譯地分割閾值。

圖1 SDWI圖像與湖泊目視解譯結果

圖2 直方圖統計與分割閾值結果

2 研究方法

2.1 改進Otsu算法

周迪等[15]提出了改進Otsu閾值分割方法并命名為KOtsu閾值法。KOtsu考慮了目標類在圖像中約占的權重，并結合了類內方差的影響，可以有效地解決Otsu算法會偏向方差大的問題。KOtsu閾值分割方法已成功實踐在Lena圖像上，本文將該方法應用到SDWI圖像，進行閾值分割。經實驗發現，類內方差的加入將更依賴權重K的取值，甚至要設置到極低權重值，才能平衡Otsu算法閾值偏向問題。因此，本文對KOtsu進行了改造，如式(1)所示，改造算法沒有考慮類內方差，引入了權重K的思想。當使用原始SAR影像進行閾值分割時，水體類在閾值的左邊，非水體類在閾值的右邊，而經過SDWI運算后，水體類和非水體類則反之。本文提出的改進Otsu與傳統的Otsu相比，其主要的變化在于閾值判斷公式上的變化，表達如式(1)所示。

T=argt∈{0，1，…，L-1}max(Kω0(t)(μ0(t)-μ(t))2+ω1(t)(μ1(t)-μ(t))2)

(1)

式中：K為目標約占圖像的權重，本文權重K取0.2；t為所遍歷的分割閾值；L為灰度級， 取值范圍為[0,256]；μ(t)為圖像總灰度均值；μ0(t)為背景類灰度均值；μ1(t)為目標類灰度均值；ω0(t)為背景類的比例；ω1(t)為目標類的比例；T為最佳分割閾值。

2.2 形態學優化算法

連通區標定算法是形態學中常用的方法之一，該方法在四連通和八連通區域對同質區域進行標定，同時記錄標定的坐標(x，y)和標記號i(i=1，2，…，n)，把不符合連通閾值的標記號刪除。相干斑造成的雜散點通常由若干像素組成，而水體則是大面積連通區，因此通過設定較小的連通閾值可以把雜散點有效去除。

在Otsu和改進的Otsu的提取結果中，把真實水體、山體陰影和雜散點都分割成了水體，并且山體陰影形成的偽水體比部分真實水體連通面積大，使用連通標定算法把山體陰影去除，部分真實水體也會被隨之刪除，為解決該問題，本文提出了高程連通區標定算法。

圖3(a)為四連通像元表示，圖3(b)為八連通像元表示。設改進Otsu算法分割結果用F表示，那么像元值可表示為f(i，j)。F是一幅二值圖，用0和1表示，其中0為非水體，1為水體。當f(i，j)=1時，使用四連通或者八連通方法對水體區域進行標定，同時增加了高程h進行判斷，即當連通區域的像元滿足式(2)時，記錄其坐標(i，j)，將F中對應坐標的值進行刪除，賦于0值，表達如式(2)所示。

f(i，j)≥h

(2)

式中：i為F的行；j為F的列；h為高程閾值。

圖3 連通方法

(3)

當滿足式(2)時，連通像元被認為是山體陰影造成誤分的偽水體，利用式(3)更新二值圖中對應坐標的值，即將值賦給0值，偽水體變更為非水體。

3 結果與分析

基于閾值分割法的水體提取，通常分為兩步：第一步，水體粗提取；第二步，水體精提取。文中使用改進Otsu閾值分割法粗提取水體，獲得初始水體分布圖，然后利用形態學方法優化初始水體分布圖，獲得最終的水體分布圖。

3.1 基于改進Otsu的水體粗提取

圖4和圖5分別為Otsu提取的水體結果和改進Otsu提取的水體結果。

圖4 Otsu水體粗提取

圖5 改進Otsu法水體粗提取

由圖4和圖5可知，兩種方法對大面積水體流域都有效地進行了提取，但是在東南和東北山區，Otsu提取的偽水體顯然多于改進Otsu。由于水體和非水體在整個研究區域面積分布并不均勻，造成分割的閾值將會偏向非水體，所以Otsu法在東南和東北山區提取的偽水體不僅僅包括了山體陰影還包括了誤分的非水體值。直接觀察圖4和圖5可知，改進Otsu提取水體的結果要優于Otsu提取結果。統計兩種方法提取水體在研究區的占比，Otsu法提取的水體占整個研究區的59.53%，改進Otsu占19.69%。通過圖2中SDWI圖像的直方圖可知，背景類和目標類是不平衡的，背景類的方差要遠大于水體類的方差，此時Otsu算法的計算閾值會偏向背景類，見圖2中T1，根據T1閾值進行分割，會出現大量的偽水體。如圖1所示，觀察SDWI圖像可判斷出水體在研究區的占比不會超過50%，顯然改進Otsu比Otsu提取水體更為合理。

谷鑫志等[16]把目視解譯閾值作為水體與非水體的真實閾值，將閾值法計算的閾值與目視解譯閾值相比較，差值越小閾值分割精度越高。

目視解譯閾值利用ENVI 5.3軟件中的 Quick stats 工具和New Raster Color Slice 工具，人機交互式選取了真實閾值，即目視解譯的灰度值，目視解譯的灰度值見表2以及圖2中T3。表2中目視解譯閾值是由量化到256灰度級公式通過目視解譯灰度值反解計算得出。Otsu和改進Otsu的分割閾值通過遍歷圖2的直方圖，計算出相關的參數后，利用Otsu閾值判斷公式和式(1)計算得出。

表2 閾值比較

由表2可知，Otsu法計算的閾值與目視解譯閾值相差較大，改進Otsu法計算的閾值比Otsu更接近目視解譯閾值，顯示了本文改進Otsu閾值分割方法的有效性。綜上所述，把改進Otsu水體粗提取圖作為初始水體分布圖，并使用形態學方法進行優化。

3.2 基于形態學優化的水體精提取

圖6顯示了采用結合高程連通標定的形態學優化方法對水體的精細提取流程。如圖6所示，把改進Otsu法提取的初始水體分布圖和DEM作為輸入，利用高程連通標定算法，記錄大于設置高程閾值的坐標，在粗提取結果中刪除被標定坐標位置的偽水體。精提取流程圖中，虛線矩形內是文獻[4]優化方法，兩次連通標定與兩次取反的意義：第一次連通標定的作用是刪除雜散點和殘留的山體陰影；第一次取反把水體轉為了背景，非水體轉為了目標；第二次連通標定的目的是把水面上分布的船只和雜散點進行刪除；第二次取反把真實水體和非水體恢復。

圖6 水體精提取流程圖

1)高程連通區標定與刪除。高程連通標定刪除的核心在于高程閾值h的設置。對圖5矩形區域，設置了六個高程閾值探索山體陰影、真實水體與設置高程閾值h之間的規則，這六個高程閾值分別為80、90、100、110、120和130 m。處理結果如圖7所示，其中圖7(a)至圖7(f)所示分別為高程閾值80、90、100、110、120和130 m的截取結果，圖7(g)為改進Otsu法在矩形區水體粗提取結果。圖7(a)至圖7(g)橢圓區通過目視解譯有三個湖泊，矩形區有兩個湖泊。圖7(a)橢圓區內，山體陰影和部分雜散點被大量剔除，隨著高程閾值的增加,橢圓區域的偽水體數量也逐漸地增加，如圖7(b)至圖7(f)所示。圖7(b)至圖7(g)矩形區兩個湖泊分別位于左上和右下方，高程閾值為80 m時，兩個湖泊被剔除，隨著高程閾值的升高，湖泊輪廓也越加明顯，提取效果較好的是圖7(e)和圖7(f)。圖7(e)與圖7(g)對比可見，矩形右下方湖泊部分被剔除，湖泊輪廓缺失。圖7(f)與圖7(g)對比可見，矩形區的兩個湖泊被完整地保留，并且區域內大部分偽水體被剔除。

圖7 規則探索實驗

綜上所述，設置的閾值越小山體陰影的剔除能力越強，但是不可避免真實水體會被刪除或者輪廓不完整，隨著高程閾值的增加，山體陰影剔除能力減弱，水體保留效果增強。因此，所設置的高程閾值需要滿足兩個規則：真實水體能夠有效地保留、對山體陰影能有效剔除。

為了獲取適合研究區的高程閾值，基于SDWI圖像、原始SAR影像，使用ArcGIS軟件矢量化圖1中目視解譯的小中型湖泊并提取了山體陰影。矢量化的湖泊和山體陰影部分區域，如圖8所示。

圖8 矢量化湖泊與提取的山體陰影

研究區高程范圍在-4.546～837.152 m之間。通過遍歷研究區高程范圍選取高程閾值，對矢量化水體和提取的山體陰影執行高程連通標定刪除。為了簡化算法的執行時間，高程閾值范圍設置為[-5，837]之間，以1為步長，對高程連通標定刪除算法進行全局搜索。基于有效保留水體和有效剔除山體陰影兩個規則，統計了每一個高程閾值保留水體像素數量和存在山體陰影像素數，如圖9所示。圖9中，橫坐標為高程閾值，縱坐標為像素數，即水體和山體陰影在某個高程閾值被保留的像素個數。從圖9可知，真實水體和山體陰影隨著高程閾值的增加，水體和山體陰影保留像素數也隨之增加。經實驗，高程閾值為219 m時，水體被完整保留，見圖9中紅虛線h1。

圖9 高程閾值全局搜索

借鑒機器學習中錯誤率公式，可以得到水體的損失率以及山體陰影的剔除率。損失率是指真實水體被刪除的像素數與總水體像素數之比。剔除率是指山體陰影被刪除的像素數與總山體陰影像素數之比。損失率和剔除率可以共用一個公式表示，但是意義卻不相同。損失率越大真實水體被保留得越少，剔除率越大山體陰影被刪除的能力越強。式(4)即損失率和剔除率的公式。

(4)

式中：S為真實水體或者山體陰影總的像元數；W為處于某個高程閾值時，水體被保留的像元數或者山體陰影存留像元數。

每個高程閾值對應的損失率和剔除率，如圖10所示。圖10中，p1、h1、h2的值分別為0.1、219、109，當高程閾值為h1時，即高程閾值為219 m時，損失率為0.00%，剔除率為30.33%，可見當取得最小的損失率時，剔除率并不高。當高程閾值取h2時，剔除了率達到了73%，其損失率接近于10%，由此可見，在放棄一定損失率時，能夠提高剔除率。表3列出了具有代表性的12個高程閾值與損失率、剔除率的關系。從表3可以看到，當高程閾值為64 m時，剔除率超過了90%，損失率為40.41%;當高程閾值為162 m時，山體陰影的剔除了超過了50%;當高程大于等于832 m時，高程連通標定刪除算法失效。

圖10 損失率與剔除率

表3 高程閾值與損失率、剔除率

可見，當高程小于832 m時，高程連通標定刪除算法對真實水體有一定的保護作用，也兼顧一定的山體陰影剔除能力，但最為有效的高程閾值并不能從表3中得到。利用從表3中的12個高程閾值執行兩次連通標定與兩次取反的結果計算與矢量化的水體交并比(intersection over union，IoU)。

當IoU取最大值時，即為最佳的高程閾值，表達如式(5)至式(6)所示。

(5)

H=argh∈[h1，h2，…，hn]max(IoU)

(6)

式中：A為矢量化小中型湖泊的矢量化像元；M(·)為水體形態學優化結果；h為高程閾值；l、d為損失率和剔除率；H為最佳高程閾值。

2)兩次連通標定和兩次取反。如圖6所示，兩次連通標定與兩次取反是水體精提取的處理環節之一。從表3可知，當高程閾值為832 m時，高程連通標定刪除算法失效，此時進行兩次連通標定和兩次取反處理等同于文獻[4]優化方法。對表3中12個高程閾值進行兩次連通標定與刪除，并計算IoU，計算結果如表4所示。表4給出了兩次連通標定與兩次取反的參數設置。第一次連通值是連通標定算法的連通像元數，當被標記的連通像元小于等于該連通像元時，刪除該連通區域。第二次連通值的作用與第一次連通值的作用相同，但是第二次的目的是為了刪除水面上的復雜目標，并且水面上復雜目標物，不管處于哪個高程閾值下都不會發生變化，因此該連通值在不同的高程閾值中是不變值。對IoU而言，如圖11所示，當高程閾值約為162 m時，得到了最高的IoU。從表3可知，高程閾值為162 m時，損失率為0.72%，剔除率為 50.19%，因此，犧牲較小的損失率，可以極大提高山體陰影的剔除率；剔除率的提升可以改變偽水體的連通區結構，使得可以使用較小的參數刪除偽水體，可以有效地保留真實水體。

表4 各高程閾值的IoU與參數設置

圖11 交并比曲線

3.3 結果精度評價

本文提出高程連通標定和連通標定算法相結合的形態學優化方法本質上不改變改進Otsu法粗提取水體的邊界和大流域連通水體的面積，其作用是把山體陰影、雜散點和水面上船只去除，得到精細的水體分布圖。綜合Sentinel-1A影像與SDWI圖像目視解譯了研究區域139個小中型潛在湖泊，并對其進行了矢量化，湖泊位置如圖1所示。

從圖11可知，隨著高程閾值變大，IoU呈先增加后減少的趨勢，因此IoU曲線會存在一個極大值點，當高程閾值約為162 m時，IoU達到了最大值，并且第一次連通標定參數設置與剔除率有一定的相關性，如表3、表4所示，隨著剔除率的增加，偽水體連通性變差，需要連通刪除的閾值也隨之變小。通過表4可知，高程閾值選取在108～219 m之間的IoU可以達到0.8以上。可視化高程閾值為162 m時的最終結果和文獻[4]結果，如圖12所示。

圖12 水體提取最終結果

圖12(a)為本文選取162 m的高程閾值最終優化結果，圖12(b)為文獻[4]優化結果。從圖12可知，本文方法相比文獻[4]優化方法，在大范圍水體提取時能夠保留更多的細節，對小中型潛在的湖泊有較高提取率。并且從表3可知，高程閾值在74～219 m時，優化結果都高于文獻[4]優化方法，這也表明高程連通標定刪除是一種有效的優化方法。

4 結束語

本文提出了改進Otsu法和高程連通標定的優化算法，結果表明，改進Otsu法閾值分割精度高于傳統Otsu法閾值分割精度，在進行優化時，使用高程連通標定可以很好地保留小中型湖泊，得出結論如下。

使用改進Otsu算法和傳統Otsu法對SDWI圖像的分割閾值進行了計算。由表2可知，目視解譯的分割閾值為226，傳統Otsu閾值分割法分割閾值為210，改進Otsu閾值分割方法分割閾值為225。改進Otsu閾值分割結果與目視解譯閾值相差1，而傳統Otsu閾值分割法相差16，顯然改進Otsu法對目標和背景不平衡的情況下，閾值分割效果更顯著。

設置了六個高程連通閾值，分別為80、90、100、110、120和130 m，討論高程連通閾值應滿足兩個規則。通過矢量化的真實水體和提取的山體陰影，設置步長為1，對高程閾值進行全局搜索，得到當高程閾值大于等于219 m時，可以完整地保留真實水體，但山體陰影剔除能力變弱，當高程大于等于832 m時，高程連通標定算法失效。

所提出的結合高程連通標定算法的形態學優化方法，區別于文獻[4]優化方法，是一種新的組合優化方法，通過計算最終結果的IoU，得到了IoU曲線。IoU曲線表明，研究區最佳的高程閾值約為162 m，108～219 m之間的高程閾值IoU在0.8以上。當高程閾值為162 m時，IoU為0.882，文獻[4]優化后的IoU為0.608。

本文結合了改進Otsu法和形態學優化方法，提取了彭澤縣區域于2020年7月20日的精細水體分布圖。本文方法在提取大范圍水體時，可以有效解決山體陰影與小型湖泊被一同刪除的問題，但仍然有少量小面積湖泊被刪除。