基于Frangi濾波的長江中下游河流信息提取

李化良

(鎮江市勘察測繪研究院，江蘇 鎮江 212004)

1 引 言

河流在全球水循環、工業熱源和生物生長及人類生活需求中發揮了重要的作用[1，2]。隨著全球氣候的變化，河流水體的生態特征和生物、物理性質都出現了改變，影響了河水流動、泥沙傳輸和生物組成[3]。準確有效地提取河流信息對我們的生活和發展至關重要。隨著遙感技術的不斷發展，越來越多的中高分辨率遙感數據能夠免費獲取。這些數據在時間分辨率和空間分辨率上基本滿足河流實時監測的需求，但是若要準確提取地表水體，需要高精度、魯棒性強的河流算法。水體指數算法是最常見的河流提取方法，它根據水體光譜特征差異對目標水體特征進行增強，方法簡單但精度相對較低[4～6]；光譜特征分類在河流水體中提取也很常見，例如監督分類和非監督分類算法，該算法容易忽略較少的水體[7，8]；深度學習算法在目前應用十分廣泛，識別精度高但需要較多樣本[9]。在利用遙感技術提取河流信息時，河流不同大小的尺度特征以及復雜的背景特征是影響河流提取精度的重要因素。基于以上兩點，我們提出了一種基于Frangi濾波的河流信息提取算法[10]，考慮到近些年我國政府加大對長江流域的監管與保護，我們將研究區選取在長江中下游流域。

2 數據源及研究區

長江流域作為產業和城鎮布局的重要載體長期以來都是我國重要經濟區，準確、高效提取長江流域水體對于監測長江岸線和保護長江生態意義重大[11，12]。隨著遙感技術的發展，近些年國內外提供了一系列分辨率較高、重訪周期較短的衛星數據，如美國的Landsat系列衛星(Landsat 5、Landsat OLI)、歐空局的哨兵系列衛星(Sentinle-1，Sentinel-2A/B等)以及中國環境衛星(HJ)、資源衛星數據(ZY-1，ZY-3等)和高分系列數據(GF)。多源衛星數據為監測長江水體提供了豐富的基礎數據源，通過高精度河流提取算法能夠對遙感影像中的河流信息進行精確提取。為了監測全球氣候變化對人類和環境影響，歐洲航空局實施了“伽利略”計劃，目前共發射Sentinel-1，Sentinel-2A/B，Sentinle-3和Sentinel-5等一系列衛星。作為美國Landsat系列衛星的后續衛星，Sentinel-2A/B在空間定位和影像信噪比上相較于Landsat OLI都有較大改進[13]。Sentinel-2A/B共包括13個波段，分別為4個10m分辨率、6個 20 m分辨率和3個60 m分辨率波段。鑒于長江東西流向很長(約 6 300 km)，本文使用 10 m分辨率(波段2：490 nm，波段4：665 nm和波段8：842 nm)的Sentinel-2A/B影像來提取長江中下游區域。

3 基于Frangi濾波的河流信息提取算法

遙感影像識別地面河流水體時，河流形狀和復雜的背景是影響算法提取的兩個重要因素。因此，首先通過水體指數來初步增強河流水體特征；再使用Frangi濾波來進一步突出線狀河流水體特征；最后利用水平集分割對河流進行提取，算法流程如圖1所示。

圖1 基于Frangi濾波的河流提取算法流程圖

3.1 水體指數算法

水體指數是河流識別中常用方法，其原理為近紅外或中紅外波段水體的強吸收而導致的吸收谷，其簡單有效且對面積較大的水體識別效果較好。目前，最常用的水體指數算法包括歸一化差異水體指數(Normalized Difference Water Index)、改進的歸一化差異水體指數(Modified NDWI，MNDWI)和增強水體指數(Enhanced Water Index，NWI)等[14]。本文利用Sentinel-2A/B MSI中3和8波進行水體指數運算：

(1)

其中band3和band8分別對應Sentinel-2A/B MSI影像的第3和第8波段。本文以江蘇境內的長江流域為例來闡述河流信息提取的流程(圖2(a))。經過波段運算，圖2中的長江流域水體變得高亮顯著。但從圖2(b)中河流邊界處水陸界線十分模糊，簡單的閾值分割方法無法準確提取干流附近細小支流。

圖2 (a)長江流域真彩色圖(b)NDWI圖

3.2 Frangi濾波算法

Frgangi濾波在醫學影像中常常被用于視網膜、血管等線性幾何結構的識別。視網膜和血管的形態與河流形狀本質上是一致的，本文創新性地將Frangi濾波推廣到陸地河流提取中。對原始影像做高斯運算，能夠得到特定尺度下高斯濾波圖像，公式如下：

Iσ(p)=Iσ(x，y)=I0(x，y)?Gσ(x，y)

(2)

其中，Iσ表示尺度為σ濾波處理后的圖像，p=(x，y)為圖像中像元位置，為卷積運算符，Gσ(x，y)代表標準偏差為σ的二維高斯核函數。其中，σ∈[σmin，σmax]，σmin和σmax分別表示最小和最大像元尺度值。

(3)

差分運算是識別數字圖像中亮度變化差異明顯鄰域像元最常用算法，差異較大像元表現為物理性質變化、灰度特征不連續。不同尺度差分運算的基本形式由公式(3)表示：

(4)

利用二階高斯濾波運算對水體指數圖像進行處理，計算每一像元在特定尺度下高斯矩陣，為式(4)：

(5)

由于Hessian矩陣是對稱矩陣，故存在兩個實數特征值λ1和λ2及對應的特征向量e1和e2，特征值λ1、λ2表示像元灰度值變化大小和趨勢。水體像元在圖像中通常表現為高亮特征且背景低暗，此時λ1≈0且λ2<0[15]。

(6)

其中參數β，γ用來調節濾波敏感性的Rb和S。

根據尺度空間理論只有像元在最佳尺度因子時Lσ才能夠達到最大值。因此，在多尺度線狀特征提取過程中將不斷進行多尺度運算來得到最大尺度特征值。

(7)

圖3是不同尺度σ下線性河流的增強情況，可以發現σ=1時只能增強細小的水體(圖3(a))；隨著尺度的增加河流越來越明顯清晰，在σ=6時已經能識別出來較寬的河流。由于長江太寬，使用σ=9無法對其進行完全識別，需要更大的尺度因子(圖3(d))。為了節省計算時間，我們對水體指數圖像(NDWI)取閾值選取較大閾值后直接提取明顯河流，然后將得到的顯著河流與Frangi濾波得到的結果相疊加，得到最終的圖像。

圖3 不同尺度σ下線狀河流的增強情況

3.3 水平集分割

(8)

xN表示像元x的鄰域，K(·)代表靈活窗口方程來保證臨近像元非負性，b用來評估影像密度與真實情況反射率偏差。

為簡化計算，該方程能夠通過水平集轉換進行區域離散化：

(9)

Mi(·)是聚類中心Ci類別方程，φ(·)為水平集分割算法。

對圖3得到的結果做水平集分割，得到線狀水體的分布，然后將所得結果與NDWI影像中提取的高亮河流合并，然后進行噪聲去除得到最終結果。

4 河流提取結果及評價

圖4為最終提取結果，圖中不同尺度的河流都被提取出來，細小的河流也被有效地識別。寬度較寬的區域河流呈現高亮特征，細小的區域出現比較暗絲狀結構，對最終得到的結果進行精度分析。進行精度分析時，對圖2(a)真彩色圖像中的水體進行目視解譯做矢量化結果圖，作為精度分析結果的“真值”，然后在ENVI軟件中使用混淆矩陣對河流提取結果進行精度分析總體精度為93%，Kappa系數為0.86。與監督分類算法中支持向量機(Support vector machine，SVM)提取長江下游區域河流(圖4(a))，本算法能夠提取更為細小的河流和溝渠，河流連通性也更好。最終結果基本滿足政府部分對農業灌溉、城鎮居民用水和長江流域水體生態保護等方面需要。

圖4 (a)支持向量機提取河流(b)基于Frangi濾波的提取河流

5 結 語

河流對全球變化和地方環境影響十分敏感，對其長度和寬度的提取一直都是比較困難的。本研究基于Frangi濾波發展了一種適用于長江中下游流域線狀水體的算法，并在該區域得到了較好的結果。由于我們的算法對線狀的特征十分敏感，該算法的應用不應當僅限于線狀河流，其對道路和房屋建筑的識別也有一定的啟示意義。