基于SLIC的哨兵1號雷達數據水體信息提取

湛 南 渝,李 小 濤,路 京 選,許 文 波

(1.電子科技大學 資源與環境學院，四川 成都 611731； 2.中國水利水電科學研究院,北京 100038)

星載合成孔徑雷達(Synthetic Aperture Radar，SAR)因具有全天時、全天候、大覆蓋面積、高分辨率、成像不受云霧影像等特點，被廣泛應用于農業、海洋、災害監測和軍事領域。

在SAR圖像中，由于鏡面反射的物理特性，水體的后向散射系數較低。SAR圖像固有的斑點噪聲以及波浪和植被反射的影響，使得水體表面出現一定程度的明暗變化。山區陰影和城市建筑物陰影，呈現出跟水體一樣的低后向散射系數，這些干擾信息給水體的精確提取造成了較大的困難[1]。

目前常用于SAR圖像水體信息提取的方法是基于閾值的分割方法[2-3]，其中Otsu閾值法[4-7]最為常見，該方法具有簡單快速的優點，但對于復雜的影像，純粹的單閾值法提取效果并不好，因此如何自動確定合適的閾值一直是一個難題。監督分類方法[8-9]是另外一種提取SAR圖像水體信息的方法，但該方法需要手動選擇訓練樣本、較復雜，提取精度與訓練樣本和分類算法密切相關。

SLIC(Simple Linear Iterative Cluster)[10]算法根據色彩和距離相似程度對圖像進行分割，具有耗時短，生成的超像素塊大小均勻、輪廓規則等特點，被廣泛應用在彩色圖像、光學遙感圖像分割中。本文針對Sentinel-1 SAR圖像的波段特征，簡化了SLIC[11]應用于SAR圖像分割，并對分割的后超像素塊進行分層區域合并，然后通過二值化得到水體信息。

1 SAR圖像SLIC超像素分割

1.1 SLIC超像素分割算法

SLIC以顏色和距離特征對圖像進行聚類，首先將圖像的顏色空間由RGB轉為LAB，將LAB色彩空間的3個分量和距離x,y組成5維特征向量，然后使用K-Means按照一定的空間度量對圖像進行聚類，生成一個個大小較為均勻的超像素，該算法時間復雜度為O(N)(N為圖像的像素個數)，計算效率較高。SLIC算法步驟如下[12]。

(2)調整種子點位置。為了避免將初始種子點定位落在物體邊緣上，以及防止初始種子點落在噪聲上，計算初始種子點n×n鄰域(一般n=3)的梯度值，調整初始種子點的位置為其梯度最小值所處位置，并為每一個初始種子點賦上不同的標簽。

(3)計算距離度量。計算每個像素點pi與距離其最近的種子點qk的距離，距離度量由顏色距離和空間距離組成，計算公式如下：

(1)

(2)

(3)

式中，dc為在LAB顏色空間的歐式距離；ds為在x，y空間上的歐式距離；D為像素點與種子點的總距離，D越小，說明像素點與該種子點之間越相似；m為緊湊度因子，用來調整顏色信息和距離信息的相對重要程度，對于LAB顏色空間來說，一般m=[1,40]；S為種子點之間的距離。

(4)局部迭代聚類。為了增加算法的運算效率，SLIC在對每一個種子點聚類時，將算法的計算范圍限制在了2S×2S的鄰域范圍內。在范圍內計算所有種子點的距離，如果該點的距離小于之前計算的距離，則認為該點在這個超像素中，并對該點的標簽值和距離值重新賦值。重復迭代(1)～(4)步驟，直到誤差收斂到某一閾值或者最大的迭代次數后，停止計算。實驗中發現，在多數情況下，10次迭代對于絕大多數影像都能取得比較好的分割效果。

(5)增強連通性。經過上述步驟后，會出現多連通、超像素尺寸過小、單個超像素被切割成多個不連續超像素等情況，需通過增強連通性來解決此問題。主要思路為：新建一張標記表，表內元素均為-1，按照“Z”型走向將面積小于理想超像素一半大小的超像素重新分配給鄰近的超像素，同時賦予相應的標簽，直到所有點遍歷完成。

1.2 基于SAR的超像素分割算法

傳統的SLIC超像素分割方法基于彩色圖像RGB，雖然理論上可以把RGB 3個波段都賦值為灰度信息，但由于運算時將RGB顏色空間轉LAB顏色空間并不適用，因此本文采用灰度信息替代彩色信息的[13]。簡化后的算式如下：

(4)

(5)

式中，dgray為在灰度顏色空間的歐式距離；D為像素點與種子點的總距離。

1.3 區域合并

由于SLIC只是把相似的塊進行了聚類，并不能直接將圖像進行二分化，所以可以先根據一定合并規則逐步將相似度較高的塊進行合并，從而只剩下相似度差值較大的區域，再使用閾值分割實現水體信息的提取。常采用的方法是層次區域合并(Hierarchical Stepwise Region Merging，HSRM)[14]，它基于區域鄰接圖(Region Adjacency Graph，RAG)進行區域合并，用相似度來判斷兩個鄰近區域是否融合。區域就是初始分割后，屬于同一標記的像素集合。K個區域的RAG表示為G=(V,E)，其中頂點集V=1,2,…，K，對應的邊集合E?V×V。每個頂點節點代表一個區域，邊表示和鄰近區域相連，邊權重就是兩個區域之間的相似度，如圖1所示。

相似度的計算公式為

(6)

計算所有區域之間的相似度，找到RAG的邊集合中最相似的邊(Cij最小)，融合最相似邊對應的兩個區域的頂點，得到K-1的RAG，更新頂點和邊集合，直到最小的相似度大于一個閾值，停止合并。

圖1 區域鄰接示意

2 實驗結果與分析

2.1 數據源與研究區域

Sentinel-1號衛星是歐空局發射的第一個哥白尼計劃衛星，由A、B兩顆衛星組成，兩個衛星的組合，可以達到每6 d對同一地點進行影像獲取。其中C波段衛星，有4種成像模式：超精細模式(Strip Map Mode，SM)、干涉寬幅模式(Interferometric Wide Swath，IW)、超寬幅模式(Extra-Wide Swath Mode，EW)、微波模式(Wave-Mode，Wave)。對地球陸地進行覆蓋的主要為IW模式，SM模式主要用于應急事件，EW和Wave模式主要用于海洋監測。Sentinel-1號衛星數據可以免費獲取，IW模式的寬幅為250 km，分辨率為5m×20 m，從衛星拍攝到數據分發至數據庫只要3～6 h左右，對地觀測能夠穿透云霧，不受天氣影像，這些特性使其非常適合應用于洪澇災害遙感監測。

為了驗證方法的有效性，基于2018年7月29日Sentinel-1A GRDH的IW的影像，選擇了湘江湖南省段的4個實驗區域進行水體信息提取，這4個區域的水體占比依次降低。并用同一天的光學影像Landsat8提取的水體作為真實水體進行精度評估。

2.2 超像素分割參數分析

在超像素分割算法中，參數的設置對于超像素的分割效果起著十分重要的作用，參數主要有兩個：超像素個數K和緊湊度因子m。緊湊度因子主要用于調節距離信息和顏色信息的相對重要程度，對于灰度值為0～255 RGB的影像而言，m一般取值為[1,40]。為探究緊湊度因子對分割的影響，采用固定超像素個數，改變緊湊度因子的方式。當m太小時，分割就會失??；當m太大時，顏色信息占比很小，距離基本由空間距離決定，圖像會被分割成均勻的方塊，結果如圖2所示。

圖2 不同緊湊度因子對于超像素分割的影響

根據超像素分割的原理，超像素的個數由超像素的邊長N決定，通過分層區域合并，使用閾值二值化，得到最終水體提取結果。不同邊長N對應的水體信息提取結果見圖3。

圖3 不同超像素邊長提取結果

從圖3中可以看出，分割出的水域面積主要取決于分割目標的大小，也就是超像素的數量。由于小于這個面積的超像素塊在增強聯通性的時候被合并，所以理論上能提出來的水域大于超像素大小一半以上的水體。即當邊長較大時，可以提出大面積的水域，當邊長較小時，可以提出細小的水體和河流，但分割時間會增加。所以對于不同分辨率的影像和不同的分割目標，應設置不同的超像素邊長。

考慮sentinel-1中的水體顯示特征和效率的因素，本文的參數設置為：超像素的邊長N=10，m=10。

2.3 與Otsu閾值法對比

為了驗證本方法的有效性，選擇了4個水體占比依次減少的實驗區域，將本方法的提取效果與Otsu閾值法的提取效果進行對比分析，結果如圖4～7所示。

圖4 區域1水體信息提取結果及灰度直方圖

區域1為寬河道水體，水域面積占比約為12%，此時圖像直方圖顯示為顯著的雙峰特性，兩種方法在主河道的提取效果上均很好，Otsu方法提取出了細小的水體，本文提出的方法提取結果比較純凈，提取的準確率高于Otsu閾值法。

圖5 區域2水體信息提取結果及灰度直方圖

區域2為寬河道與細小河道并存，此時水體占比約為7%，此時直方圖仍具有相對顯著的雙峰特性，兩種方法在主要的河道的提取效果上均保持不錯，但Otsu閾值法仍存在較多的誤提區域，提取的準確率相對較低。

圖6 區域3水體信息提取結果及灰度直方圖

圖7 區域4水體信息提取結果及灰度直方圖

區域3和區域4均為蜿蜒的細小河道，水體占比約為3%和2%，此時直方圖中雙峰特性不顯著，Otsu閾值法的提取效果非常差，誤提現象十分嚴重，而本文方法則仍能提取出大部分水體，保證了河道的完整性。

綜上所述，當圖像水域面積占比較大時，圖像具有明顯雙峰特性，兩種方法的水體提取效果均不錯，但當水體占比較少時，雙峰特性不顯著，甚至呈現單峰特性，Otsu水體提取效果很差。而本文提出的算法由于同時考慮了距離和色彩信息，在提取之前先對水域聚類，在合并時能使得聯通的水體得到保留，所以仍能提取出大部分水體。因此本文方法在水體提取的穩定性上強于Otsu閾值法。

2.4 精度評估與分析

將本文方法和Otsu閾值法提取的水體范圍和同一時間landsat8確定水體范圍進行比較，采用查全率和查準率以及綜合評價指標來評估其精度。查全率R、查準率P以及綜合評價指標F1定義為

(7)

(8)

(9)

如表1所示，Otsu算法在4個區域中的查全率都比較高，但這源于它提取了大量的錯誤水體，尤其是在區域3和區域4水體占比較小的情況下，準確率相當的低。而本文算法則在查準率上全面優于Otsu算法，綜合性能也強于Otsu算法。當圖像中水體占比較少，河道較窄時，雙峰特性不明顯，此時Otsu基本不能識別水體，但本文的方法仍能識別出大部分河道，說明本文方法有很好的魯棒性。

表1 不同方法水體提取精度對比

3 實例應用

洪水監測是水體提取的一個重要應用[15-16]。以2019年7月份株洲的洪水監測為例，驗證本文提出方法的有效性。

基于Sentinel-1A雷達影像，對受災較為嚴重的湘江和淥水的交匯地段進行動態監測數據如表2所示。

表2 監測數據情況

監測結果如圖8～10所示。監測結果顯示，2019年7月12日監測區域內湘江和淥水沿岸受災嚴重，監測范圍內，淹沒面積達到10.06 km2；到7月19日，沿岸洪水基本全部退去，河道恢復正常。

圖8 株洲洪水前后水體信息提取情況

圖9 株洲災前災中對比

圖10 株洲災中災后對比

4 結 語

本文基于Sentinel-1A SAR數據，提出了采用簡化的超像素分割方法SLIC先對圖像進行分割，再用分層區域合并準備對圖像進行合并，最后進行二值化得到水體信息的方法。并與Otsu閾值法水體信息提取進行了對比，結果表明，在水域面積占比較高和較低的場景中，本文提出的方法在綜合性能上都高于Otsu方法，能夠較好地提取出水體信息。同時以2019年7月份株洲洪水為例進行了實例應用，得到了較好的監測結果，驗證了本方法的適用性。