基于迭代閾值分割的星載SAR洪水區域快速提取

苗 添, 曾虹程, 王 賀, 陳 杰

(北京航空航天大學電子信息工程學院, 北京 100191)

0 引 言

洪澇災害以其發生的頻率高、影響的范圍廣、造成的破壞強等特點一直受到科研人員的高度關注。由于洪災成因復雜,對其進行準確預測較為困難。因此,在洪災發生后,及時地進行洪災區域提取對救災工作有著重要的意義。

目前,遙感衛星技術被廣泛應用于洪災監測中。其中,光學衛星的圖像直觀、分辨率高,但受觀測區域的氣象條件影響較大,在復雜的氣象條件下難以勝任檢測工作,且僅能在白天工作。而星載合成孔徑雷達(synthetic aperture radar, SAR)憑借其全天時、全天候工作,分辨率高,不受晝夜、天氣影響,檢測能力強、范圍廣等獨特的優勢,在近些年來被廣泛應用到對地環境監測中。基于上述的優點,本文開展了基于星載SAR圖像的洪水災害檢測處理研究。

利用星載SAR進行洪水區域提取主要是通過水體提取和圖像變化檢測來實現的。當前,較為常見的洪水區域提取方法是:利用Otsu閾值分割算法、K-Means聚類算法等圖像分割或圖像分類方式提取災前與災后水體,再通過相減變化檢測來提取洪災區域。這類方法效果穩定、可靠,提取的準確率較高,但也存在局限:由于需要進行遍歷或多次迭代,其耗時都相對較長,檢測效率較低。針對這一問題,本文提出一種基于改進迭代閾值分割的SAR圖像洪水區域快速提取方法,可在保持較高準確率的同時顯著縮短算法的處理時間,提升洪災區域提取的效率。

1 基于SAR圖像的洪災區域傳統提取方法

目前,有許多學者對遙感技術在洪災檢測中的應用進行研究,并取得了一些成果。由于光學衛星成像易受天氣及水汽云層的影響,因此多采用SAR圖像或光學圖像與SAR圖像融合的方式進行水體提取與洪災區域檢測。

對SAR來說,由于不同地表的粗糙度相差較大,雷達信號在其上的散射方式不同,導致不同地表對應的后向散射系數不同。雷達信號在水體上以鏡面反射為主,而在陸地上以散射為主。根據這一點,可利用圖像分割或圖像分類進行水體提取,再通過圖像變化檢測提取洪災區域。傳統的水體提取方法有Otsu閾值分割算法、最大熵閾值分割算法、K-Means算法、期望最大化(expectation-maximum, EM)算法等方法。近年來,隨著深度學習等領域的發展, U型網絡(U-Net)、全卷積網絡(fully convolutional networks, FCN)等深度學習模型也在SAR的水體提取中得到了應用。

Otsu閾值分割算法在水體提取中的應用較為普遍。其原理為尋找閾值,使該閾值分割所得的兩部分圖像的類間方差達到最大值:

(1)

最大熵閾值分割算法與Otsu閾值分割算法原理類似,區別在于其閾值取值為使兩部分圖像的總熵值達到最大值時的值:

(2)

式中:為圖像的總熵值;、分別代表閾值分割后兩部分圖像的熵值;、分別為兩部分圖像中各個像素值所占比例;為閾值對應的像素值等級;為像素值等級個數。利用最大熵閾值分割算法提取水體,同樣存在時間復雜度高、易受山體陰影影響等缺點。

在水體提取中,常用的圖像分類方法有K-Means算法和EM算法。K-Means算法通過多次迭代,修正聚類中心,將原始數據分為若干不同類別,直到準則函數達到最優值,從而實現圖像的分類。EM算法的原理是在概率模型中找到參數最大似然估計或最大后驗估計的算法。通過多次迭代,實現對目標的分類。

(3)

式中:為待估測值;為觀測變量值;為未知變量值;()為似然函數;()為對似然函數取對數的結果。其原理如下:通過多次迭代,改變,使似然函數最大化,優化分類結果。隨著機器學習的發展,這兩種非監督性學習在SAR圖像水體提取中也得到了較多應用,但是由于需要通過多次迭代來進行提取,因而所需時間往往比閾值分割法更長。

近年來,隨著相關理論與技術的完善,深度學習在SAR圖像處理中得到了廣泛應用。在水體提取中,研究者們根據SAR圖像自身的特點對經典的深度學習分割模型進行改進,提出了針對SAR圖像的基于U-Net、FCN-8S、雙邊分割網絡(bilateral segmentation network, BiSeNet)等分割網絡的水體提取算法,取得了較好的結果。戴牧宸等根據SAR圖像特點,對傳統的雙邊網絡進行改進,減少了原網絡中空間路徑的卷積層數,采用深度殘差網絡(deep residual network, ResNet)作為上下文路徑骨干網絡對SAR圖像進行海陸分割,實驗表明,其分割效果較好,且泛化能力較強。Shamsolmoali等基于常規的U-Net網絡,針對SAR圖像的特點進行改進,將殘差網絡加入傳統的U-Net之中。在上下采樣路徑中,加入緊密連接的殘差網絡塊,對不同尺度的上下文信息進行聚合。實驗表明,該方法較傳統的U-Net網絡而言,提取準確率得到了提高。

此外,相關學者還基于SAR圖像與光學圖像融合結果開展洪災區域提取研究,取得了較好的效果,但超出本文研究范圍,不再詳細說明。

由于洪水檢測的特殊性,對算法的時效性和準確性均有較高的要求,上述算法均需較長時間。針對這一問題,本文將開展基于星載SAR圖像的洪災區域快速提取方法研究。

2 基于改進迭代閾值分割的SAR圖像洪水區域快速提取方法

迭代閾值分割算法在光學圖像處理中得到了較多的應用,但在SAR圖像處理中應用較少。由于SAR在成像過程中會受到相干斑噪聲及各類畸變的影響,因此部分像素點會受到影響,像素值因此被改變,進而影響到閾值的選取。為了減少異常點對提取效果的影響,需結合SAR圖像特點,對噪聲進行抑制。基于迭代閾值分割算法的原理及SAR圖像處理的相關理論,對迭代閾值分割算法進行改進,以實現對SAR圖像洪水區域的快速提取,如圖1所示。

圖1 基于迭代閾值分割算法的SAR圖像洪水區域提取流程圖Fig.1 Flow chart of flood area extraction of SAR image based on iterative threshold segmentation algorithm

首先,對災前和災后的SAR圖像進行預處理,具體包括裁剪、多視處理、輻射校正、幾何校正,保證SAR圖像的相干斑噪聲、輻射畸變、幾何畸變等得到抑制;在此基礎上,對預處理后的數據進行分級、高斯擬合、抽樣,對異常點帶來的影響進行進一步抑制;其次,利用迭代閾值分割算法計算所需閾值,實現對SAR圖像進行分割處理;然后,開展形態學濾波處理,抑制圖像噪聲對水體識別帶來的影響;最后,基于相減變化檢測,實現洪災區域的有效提取。

2.1 預處理

閾值分割是完全根據圖像像素值進行分割的一項技術,因此圖像中存在的各類畸變及異常點都會對分割結果帶來影響。SAR在成像過程中會受到各類畸變的影響而產生部分異常點,進而影響后續水體的提取及洪災區域的檢測。因此,在水體提取和圖像變化檢測之前,需要對原始圖像進行預處理,以達到抑制各類噪聲與畸變的目的。預處理的流程如圖2所示。

圖2 預處理流程圖Fig.2 Preprocessing flow chart

對災前與災后圖像進行裁剪,使洪災區域位于圖像中央,方便觀測的同時減小圖像大小,加快處理速率,如圖3所示。再依次對圖像進行多視處理、輻射校正、幾何校正,分別對圖像的相干斑噪聲、輻射畸變、幾何畸變進行抑制,得到預處理后的結果。

圖3 原始圖像與裁剪后的圖像的范圍Fig.3 Range of the original image and the cropped image

2.2 改進的迭代閾值分割算法

預處理結束之后,基于迭代閾值分割算法原理,分別對災前和災后研究區域的SAR圖像進行水體提取。迭代閾值分割算法是基于無限逼近思想提出的一種閾值分割方法。最初,由Ridler和Calvard提出一種閾值計算框架,通過多次選擇不斷優化閾值,以達到最佳分割效果,如圖4所示。

圖4 用于迭代閾值選擇的示意圖圖像處理器Fig.4 Schematic image processor for iterative threshold selection

圖4中,比較器對輸入圖像像素與當前閾值進行比較。選擇器根據比較器的結果向轉換器發送信號,轉換器在接受信號后,將對應的輸入圖像像素點轉換為背景圖像像素點或目標圖像像素點;閾值平均模塊負責分別計算兩個區域的閾值并進行平均得到新的閾值。當圖像中的每個像素點都經過一次處理后,即完成一次迭代。通過多次迭代,當閾值不再變化或兩次閾值的差值小于規定值時,即完成迭代閾值分割。

Trussel對該計算框架進行數學歸納,將閾值計算公式歸納為

(4)

式中:表示第次迭代的閾值;表示灰度等級;()表示對應灰度等級的像素點個數。Magid等又對Trussel的歸納進行了證明。

假設閾值滿足:

(5)

式中:、分別為背景部分和目標部分的某一特殊數學表示。定義誤差函數:

(6)

式中:為各個像素點的灰度值;()為每個像素值等級所占比例。為使誤差函數最小,對其求導,得

(7)

計算可得

(8)

式中:表示背景部分圖像的平均灰度等級。同理,可得

=

(9)

式中:表示目標部分圖像的平均灰度等級。至此,迭代閾值分割算法的證明完成。

迭代閾值分割算法通過多次迭代不斷優化閾值,當相鄰兩次閾值相同或其差值小到預定范圍時,確定最后一次迭代的結果為最終閾值。迭代閾值分割算法流程圖如圖5所示。

圖5 迭代閾值分割算法流程圖Fig.5 Flow chart of iterative threshold segmentation algorithm

如圖5所示,其計算流程如下:將高斯擬合并采樣后的數據輸入,將初始閾值設置為所有像素點的均值,并進行初次分割。在分割后,分別計算背景部分與目標部分的平均像素值,并將二者的平均值作為新的閾值, 比較與,當二者不相等時,繼續迭代,當二者相等時,完成閾值的計算。

這一算法對于像素值呈雙峰狀分布的圖像有良好的分割效果。常見的Otsu閾值分割算法和最大熵閾值分割算法由于需要對所有像素等級進行遍歷,因此其時間復雜度較高,采用常規方法實現上述兩種算法,復雜度為()。而迭代閾值分割算法不需要對所有像素值進行遍歷,僅通過數次迭代即可達到較好的分割效果,計算速度相對更快,效率更高。

在進行SAR圖像水體提取時,由外界影響導致的異常點往往會影響提取效果,針對這一問題,提出一種改進的迭代閾值分割算法。由于水體的后向散射系數較低,陸地的后向散射系數較高,因此在理想狀態下,包含較大區域水體的SAR圖像后向散射系數分布為雙峰狀,當山體陰影影響較大時,后向散射系數的分布會產生畸變,使個別等級的像素點個數異常。針對SAR圖像自身易受山體陰影影響、存在相干斑噪聲等干擾的問題,結合SAR圖像中像素值的分布特征,可先對SAR圖像的像素值進行分級量化,對量化后的數據進行擬合再抽樣,對該類畸變進行抑制。之后,利用抽樣所得數據進行迭代閾值分割,最后再進行形態學濾波,對噪點進行抑制。其流程如圖6所示。

圖6 改進的迭代閾值分割算法流程圖Fig.6 Flow chart of improved iterative threshold segmentation algorithm

具體步驟如下:

首先,將預處理后的SAR圖像中每個像素點的像素值轉化為分貝形式：

=10lg

(10)

式中:為原始像素點的像素值;為轉化后dB圖像對應像素點的值。對轉化后圖像的像素點值進行分級、量化,分為個等級,并統計每個等級的像素點個數。之后,對其進行擬合。目前,在變化檢測技術發展過程中,基于參數模型的方法應用較廣。常用的模型有高斯模型、廣義高斯模型、瑞利-高斯模型等。考慮到高斯分布的靈活性和穩定性,且多種涉及被動傳感器的數據分布模型均為高斯模型或廣義高斯模型,本文中的SAR圖像像素值分布圖采用高斯函數進行擬合,本文選用高斯函數對數據進行擬合,擬合函數如下:

(11)

式中:()為高斯擬合函數;為原始的離散數據;、、為高斯函數中的常量;為擬合階數。依據擬合情況調整擬合階數,以達到最優擬合效果。完成擬合后,進行抽樣,將抽樣結果存入×1的一維矩陣之中。然后,再對抽樣結果進行迭代閾值分割。

迭代閾值分割算法的總體思路為:首先設置初始閾值,再不斷對其進行調整、迭代,每一次調整的輸入參數為上一次迭代中得到的結果,經過數次迭代之后,最后產生一個理想的閾值。

其算法實現步驟如下:

首先,將圖片中的像素點依據其幅值進行排序,將初始閾值設置為所有像素點像素值的平均值:

(12)

然后,按照初始閾值將原始圖像分成1、2兩部分,分別計算兩部分的后向散射系數差值的平均值:

(13)

(14)

式中:代表像素點的值;為未轉化前像素點的值;表示其對應的概率;表示閾值對應的像素等級。

接著將二者的平均值作為新的閾值:

(15)

對新的閾值與之前計算的閾值進行對比。若與的值相同,則為所求閾值;若不相同,則繼續計算,不斷迭代,直至相等或兩次閾值的差值小于設定值時,最后一次迭代的閾值即為所求閾值。根據閾值,將原圖像轉化為二值圖像。所得二值圖像即為洪水區域的提取結果:

(16)

式中:為原始像素值;為閾值;()為處理后的像素值。

2.3 形態學濾波與洪水區域提取

由于閾值分割法在處理時會忽視圖像中蘊含的空間信息,因此部分沒有得到完全抑制的異常點會給算法的判斷帶來影響,從而產生離散的噪點。針對這一問題,在完成水體初步提取后,采用形態學濾波對噪點進行抑制。由于噪點分布較為分散,為了減少形態學濾波對水體提取效果的影響,盡量選擇尺寸較小的處理元素進行處理。本文中,選取尺寸為3×3的矩形進行形態學濾波。先進行腐蝕運算,再進行膨脹運算。腐蝕運算以處理單元的中心點為錨點,在原圖像中移動處理單元,遍歷圖像的每一個像素。取結構元素覆蓋下的原圖對應區域內的所有像素的最小值,用這個最小值替換當前像素值,得到結果;膨脹運算與腐蝕運算相反,取結構元素覆蓋下的原圖對應區域內的所有像素的最大值,用最大值替換當前像素值,得到結果。經過形態學濾波,圖像的噪聲得到抑制,災前與災后的水體得到提取。

洪水區域提取的最后一步,對災前、災后圖像進行變化檢測。在之前的處理中,實現了對災前與災后水體的分別提取。在此基礎上,利用相減變化檢測可得到災后相比災前擴張的水體部分,此部分即為洪災區域。通過先提取水體再進行變化檢測的方式,可抑制零星噪點對變化檢測結果的影響。具體流程為:利用災后水體提取圖像減去災前水體提取圖像,對應相隨點值相減。對于未變化地區,相減后的結果為0,對于遭受洪澇災害的地區,相減后結果為255。所得二值圖像即為遭受洪災區域的提取結果。

3 實驗結果與分析

3.1 實驗數據

2020年7月,鄱陽湖流域發生特大洪水,給當地群眾的生命財產安全帶來了嚴重威脅。本文基于星載SAR圖像對2020年7月鄱陽湖流域的洪水進行研究。本文選用歐空局哨兵-1衛星的數據進行相關研究。選用的SAR圖像為Level-1級別,干涉寬幅(interferometric wide swath, IW)模式、地距多視影像(ground range detected, GRD)類型的單極化成像。本次使用的數據信息如表1所示。

表1 原始圖像數據信息Table 1 Original image data information

3.2 實驗結果

經過預處理的圖像如圖7所示。

圖7 預處理后的SAR圖像Fig.7 Preprocessed SAR image

提取預處理后SAR圖像像素點數據,按照其像素值大小,均勻分為500個等級,并統計每個等級像素點個數。之后利用高斯分布對數據進行擬合,對擬合效果進行評估。分別計算不同階數擬合函數的擬合系數與均方根誤差(root mean square error, RMSE),如圖8所示。

圖8 不同階數高斯擬合的效果Fig.8 Effect of different orders of Gaussian fitting

由圖8可知,數據1與數據3應采用六階高斯函數,數據2采用五階高斯函數。高斯擬合后的結果如圖9所示。常量取值如表2所示。

圖9 原始數據與擬合結果Fig.9 Original data and fitting results

表2 高斯擬合函數參數Table 2 Parameters of Gaussian fitting function

利用迭代閾值分割算法分別計算3幅圖像的閾值,計算結果如下:2019年7月20日的SAR圖像閾值等級為119,圖像中對應像素點的值為-16.18 dB;2020年3月16日的SAR圖像閾值等級為109,圖像中對應像素點的值為-16.07 dB;2020年7月14日的SAR圖像閾值等級為122,圖像中對應像素點的值為-15.36 dB。依據閾值對圖像進行分割,得到初步提取的結果,再對其進行形態學濾波,濾波前后結果如圖10所示。

圖10 形態學濾波前后水體提取結果Fig.10 Water area extraction before and after morphological filtering

對提取水體后的結果進行統計:數據1的SAR圖像中,提取水體所占比例為23.31%;數據2的SAR圖像中,提取水體所占比例為18.45%;數據3的SAR圖像中,提取水體所占比例為27.81%。與數據2相比,鄱陽湖流域的水體面積在7月增加了50.73%;即使與數據3,即2019年同期相比,鄱陽湖流域的水體面積也增加了19.30%。

分別對數據3與數據2、數據3與數據1的提取結果進行相減變化檢測,并將結果與災前、災后疊加形成的RGB (red green blue)圖像進行對比,其結果如圖11所示。

圖11 洪災區域提取結果對比Fig.11 Comparison of extraction of flood area

在RGB圖像中,紅色區域即為遭受洪澇災害地區,由圖11可知,通過水體提取和相減變化檢測,洪災區域基本實現了提取。根據相減變化檢測結果可知,與2020年3月相比,在2020年7月,鄱陽湖流域的水體面積有明顯的增大,鄱陽湖流域的西南部分擴張尤為明顯,鄱陽湖的北部區域擴張相對較少;與2019年同期相比,鄱陽湖也擴張了一部分,擴張部分集中在鄱陽湖西部區域和東南區域。

為評估水體提取結果,利用Landsat-8光學衛星開展驗證實驗。由于2020年7月鄱陽湖流域雨水較多、云層較厚,光學衛星難以觀測,因此選用2019年7月Landsat-8的圖像與提取結果進行對比,如圖12所示。

圖12 Landsat-8光學衛星2019年7月鄱陽湖區域圖像Fig.12 Landsat-8 optical satellite image of Poyang Lake area in July 2019

除了與迭代閾值分割算法進行對比外,本文還選用Otsu閾值分割算法、K-Means非監督性學習、EM算法的提取結果進行對比,如圖13所示。

圖13 不同算法的提取結果Fig.13 Extraction results of different algorithms

為驗證本文所提方法的有效性,通過采樣評估對上述3種提取方式的提取效果進行評估。具體方法為:在光學衛星圖像中分別隨機選取300個點,分別統計3種提取方式的準確率、虛警率、漏警率、計算時間等。其結果如表3所示。

表3 不同提取方式評估結果Table 3 Evaluation of different extraction methods

根據隨機采樣后的對比結果可知,對于2019年7月鄱陽湖流域的SAR圖像,迭代閾值分割法與Otsu閾值分割法的閾值較為接近,兩者的提取效果差距不大;K-Means非監督性學習提取的結果與前二者相比,準確率較低、虛警率較高、漏警率較為接近。EM算法的準確率最高,虛警率和漏警率都較低;在計算時間方面,本文使用處理器為Intel Core i7-9850H CPU的計算機完成實驗, 其中迭代閾值分割算法完成對2019年7月鄱陽湖流域SAR圖像數據的處理需要1.045 s,而Otsu閾值分割算法需要10.016 s,K-Means非監督性學習需要78.153 s。Otsu閾值分割算法由于需要對所有后向散射系數等級進行一次遍歷,因此算法的時間復雜度較高,導致所需時間較長。而對于K-Means算法和EM算法,均需要進行多次迭代,所需時間遠大于前兩種方法。對于迭代閾值分割算法,在一般情況下,僅需迭代數次即可得出閾值,所需時間大大減少,提取效率也相應得到了提高。

4 結 論

本文提出了一種基于改進迭代閾值分割的SAR圖像洪水區域快速提取方法。該方法首先通過預處理與高斯擬合對噪聲和畸變進行抑制,之后通過迭代閾值分割算法提取水體,再利用形態學濾波抑制噪聲,最后通過變化檢測提取洪水區域。與傳統的洪水區域提取算法相比,該算法在保證較高提取準確率的同時提高了處理速度,使提取效率得到提升。將該方法應用于2020年7月鄱陽湖流域的洪水災害研究中,通過預處理、水體提取、變化檢測等流程,得到了鄱陽湖流域洪水區域的提取結果。通過將結果與Landsat-8光學衛星圖像的對比,該方法的可靠性與準確性得到了驗證;通過與傳統提取方式的對比實驗,表明了迭代閾值分割算法在洪水區域提取中的高效性,為利用SAR圖像提取洪水區域提供了新的研究思路。

后續在這方面的主要工作內容為:① 對算法的適用性進行進一步檢驗,考慮復雜情況下的洪澇區域提取,如城市內澇等;② 考慮結合光學圖像信息,改進算法結構,提升準確率。