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

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

0 引 言

洪澇災(zāi)害以其發(fā)生的頻率高、影響的范圍廣、造成的破壞強(qiáng)等特點(diǎn)一直受到科研人員的高度關(guān)注。由于洪災(zāi)成因復(fù)雜,對(duì)其進(jìn)行準(zhǔn)確預(yù)測(cè)較為困難。因此,在洪災(zāi)發(fā)生后,及時(shí)地進(jìn)行洪災(zāi)區(qū)域提取對(duì)救災(zāi)工作有著重要的意義。

目前,遙感衛(wèi)星技術(shù)被廣泛應(yīng)用于洪災(zāi)監(jiān)測(cè)中。其中,光學(xué)衛(wèi)星的圖像直觀、分辨率高,但受觀測(cè)區(qū)域的氣象條件影響較大,在復(fù)雜的氣象條件下難以勝任檢測(cè)工作,且僅能在白天工作。而星載合成孔徑雷達(dá)(synthetic aperture radar, SAR)憑借其全天時(shí)、全天候工作,分辨率高,不受晝夜、天氣影響,檢測(cè)能力強(qiáng)、范圍廣等獨(dú)特的優(yōu)勢(shì),在近些年來被廣泛應(yīng)用到對(duì)地環(huán)境監(jiān)測(cè)中。基于上述的優(yōu)點(diǎn),本文開展了基于星載SAR圖像的洪水災(zāi)害檢測(cè)處理研究。

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

1 基于SAR圖像的洪災(zāi)區(qū)域傳統(tǒng)提取方法

目前,有許多學(xué)者對(duì)遙感技術(shù)在洪災(zāi)檢測(cè)中的應(yīng)用進(jìn)行研究,并取得了一些成果。由于光學(xué)衛(wèi)星成像易受天氣及水汽云層的影響,因此多采用SAR圖像或光學(xué)圖像與SAR圖像融合的方式進(jìn)行水體提取與洪災(zāi)區(qū)域檢測(cè)。

對(duì)SAR來說,由于不同地表的粗糙度相差較大,雷達(dá)信號(hào)在其上的散射方式不同,導(dǎo)致不同地表對(duì)應(yīng)的后向散射系數(shù)不同。雷達(dá)信號(hào)在水體上以鏡面反射為主,而在陸地上以散射為主。根據(jù)這一點(diǎn),可利用圖像分割或圖像分類進(jìn)行水體提取,再通過圖像變化檢測(cè)提取洪災(zāi)區(qū)域。傳統(tǒng)的水體提取方法有Otsu閾值分割算法、最大熵閾值分割算法、K-Means算法、期望最大化(expectation-maximum, EM)算法等方法。近年來,隨著深度學(xué)習(xí)等領(lǐng)域的發(fā)展, U型網(wǎng)絡(luò)(U-Net)、全卷積網(wǎng)絡(luò)(fully convolutional networks, FCN)等深度學(xué)習(xí)模型也在SAR的水體提取中得到了應(yīng)用。

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

(1)

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

(2)

式中:為圖像的總熵值;、分別代表閾值分割后兩部分圖像的熵值;、分別為兩部分圖像中各個(gè)像素值所占比例;為閾值對(duì)應(yīng)的像素值等級(jí);為像素值等級(jí)個(gè)數(shù)。利用最大熵閾值分割算法提取水體,同樣存在時(shí)間復(fù)雜度高、易受山體陰影影響等缺點(diǎn)。

在水體提取中,常用的圖像分類方法有K-Means算法和EM算法。K-Means算法通過多次迭代,修正聚類中心,將原始數(shù)據(jù)分為若干不同類別,直到準(zhǔn)則函數(shù)達(dá)到最優(yōu)值,從而實(shí)現(xiàn)圖像的分類。EM算法的原理是在概率模型中找到參數(shù)最大似然估計(jì)或最大后驗(yàn)估計(jì)的算法。通過多次迭代,實(shí)現(xiàn)對(duì)目標(biāo)的分類。

(3)

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

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

此外,相關(guān)學(xué)者還基于SAR圖像與光學(xué)圖像融合結(jié)果開展洪災(zāi)區(qū)域提取研究,取得了較好的效果,但超出本文研究范圍,不再詳細(xì)說明。

由于洪水檢測(cè)的特殊性,對(duì)算法的時(shí)效性和準(zhǔn)確性均有較高的要求,上述算法均需較長(zhǎng)時(shí)間。針對(duì)這一問題,本文將開展基于星載SAR圖像的洪災(zāi)區(qū)域快速提取方法研究。

2 基于改進(jìn)迭代閾值分割的SAR圖像洪水區(qū)域快速提取方法

迭代閾值分割算法在光學(xué)圖像處理中得到了較多的應(yīng)用,但在SAR圖像處理中應(yīng)用較少。由于SAR在成像過程中會(huì)受到相干斑噪聲及各類畸變的影響,因此部分像素點(diǎn)會(huì)受到影響,像素值因此被改變,進(jìn)而影響到閾值的選取。為了減少異常點(diǎn)對(duì)提取效果的影響,需結(jié)合SAR圖像特點(diǎn),對(duì)噪聲進(jìn)行抑制?；诘撝捣指钏惴ǖ脑砑癝AR圖像處理的相關(guān)理論,對(duì)迭代閾值分割算法進(jìn)行改進(jìn),以實(shí)現(xiàn)對(duì)SAR圖像洪水區(qū)域的快速提取,如圖1所示。

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

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

2.1 預(yù)處理

閾值分割是完全根據(jù)圖像像素值進(jìn)行分割的一項(xiàng)技術(shù),因此圖像中存在的各類畸變及異常點(diǎn)都會(huì)對(duì)分割結(jié)果帶來影響。SAR在成像過程中會(huì)受到各類畸變的影響而產(chǎn)生部分異常點(diǎn),進(jìn)而影響后續(xù)水體的提取及洪災(zāi)區(qū)域的檢測(cè)。因此,在水體提取和圖像變化檢測(cè)之前,需要對(duì)原始圖像進(jìn)行預(yù)處理,以達(dá)到抑制各類噪聲與畸變的目的。預(yù)處理的流程如圖2所示。

圖2 預(yù)處理流程圖Fig.2 Preprocessing flow chart

對(duì)災(zāi)前與災(zāi)后圖像進(jìn)行裁剪,使洪災(zāi)區(qū)域位于圖像中央,方便觀測(cè)的同時(shí)減小圖像大小,加快處理速率,如圖3所示。再依次對(duì)圖像進(jìn)行多視處理、輻射校正、幾何校正,分別對(duì)圖像的相干斑噪聲、輻射畸變、幾何畸變進(jìn)行抑制,得到預(yù)處理后的結(jié)果。

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

2.2 改進(jìn)的迭代閾值分割算法

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

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

圖4中,比較器對(duì)輸入圖像像素與當(dāng)前閾值進(jìn)行比較。選擇器根據(jù)比較器的結(jié)果向轉(zhuǎn)換器發(fā)送信號(hào),轉(zhuǎn)換器在接受信號(hào)后,將對(duì)應(yīng)的輸入圖像像素點(diǎn)轉(zhuǎn)換為背景圖像像素點(diǎn)或目標(biāo)圖像像素點(diǎn);閾值平均模塊負(fù)責(zé)分別計(jì)算兩個(gè)區(qū)域的閾值并進(jìn)行平均得到新的閾值。當(dāng)圖像中的每個(gè)像素點(diǎn)都經(jīng)過一次處理后,即完成一次迭代。通過多次迭代,當(dāng)閾值不再變化或兩次閾值的差值小于規(guī)定值時(shí),即完成迭代閾值分割。

Trussel對(duì)該計(jì)算框架進(jìn)行數(shù)學(xué)歸納,將閾值計(jì)算公式歸納為

(4)

式中:表示第次迭代的閾值;表示灰度等級(jí);()表示對(duì)應(yīng)灰度等級(jí)的像素點(diǎn)個(gè)數(shù)。Magid等又對(duì)Trussel的歸納進(jìn)行了證明。

假設(shè)閾值滿足:

(5)

式中:、分別為背景部分和目標(biāo)部分的某一特殊數(shù)學(xué)表示。定義誤差函數(shù):

(6)

式中:為各個(gè)像素點(diǎn)的灰度值;()為每個(gè)像素值等級(jí)所占比例。為使誤差函數(shù)最小,對(duì)其求導(dǎo),得

(7)

計(jì)算可得

(8)

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

=

(9)

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

迭代閾值分割算法通過多次迭代不斷優(yōu)化閾值,當(dāng)相鄰兩次閾值相同或其差值小到預(yù)定范圍時(shí),確定最后一次迭代的結(jié)果為最終閾值。迭代閾值分割算法流程圖如圖5所示。

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

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

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

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

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

具體步驟如下:

首先,將預(yù)處理后的SAR圖像中每個(gè)像素點(diǎn)的像素值轉(zhuǎn)化為分貝形式：

=10lg

(10)

式中:為原始像素點(diǎn)的像素值;為轉(zhuǎn)化后dB圖像對(duì)應(yīng)像素點(diǎn)的值。對(duì)轉(zhuǎn)化后圖像的像素點(diǎn)值進(jìn)行分級(jí)、量化,分為個(gè)等級(jí),并統(tǒng)計(jì)每個(gè)等級(jí)的像素點(diǎn)個(gè)數(shù)。之后,對(duì)其進(jìn)行擬合。目前,在變化檢測(cè)技術(shù)發(fā)展過程中,基于參數(shù)模型的方法應(yīng)用較廣。常用的模型有高斯模型、廣義高斯模型、瑞利-高斯模型等?？紤]到高斯分布的靈活性和穩(wěn)定性,且多種涉及被動(dòng)傳感器的數(shù)據(jù)分布模型均為高斯模型或廣義高斯模型,本文中的SAR圖像像素值分布圖采用高斯函數(shù)進(jìn)行擬合,本文選用高斯函數(shù)對(duì)數(shù)據(jù)進(jìn)行擬合,擬合函數(shù)如下:

(11)

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

迭代閾值分割算法的總體思路為:首先設(shè)置初始閾值,再不斷對(duì)其進(jìn)行調(diào)整、迭代,每一次調(diào)整的輸入?yún)?shù)為上一次迭代中得到的結(jié)果,經(jīng)過數(shù)次迭代之后,最后產(chǎn)生一個(gè)理想的閾值。

其算法實(shí)現(xiàn)步驟如下:

首先,將圖片中的像素點(diǎn)依據(jù)其幅值進(jìn)行排序,將初始閾值設(shè)置為所有像素點(diǎn)像素值的平均值:

(12)

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

(13)

(14)

式中:代表像素點(diǎn)的值;為未轉(zhuǎn)化前像素點(diǎn)的值;表示其對(duì)應(yīng)的概率;表示閾值對(duì)應(yīng)的像素等級(jí)。

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

(15)

對(duì)新的閾值與之前計(jì)算的閾值進(jìn)行對(duì)比。若與的值相同,則為所求閾值;若不相同,則繼續(xù)計(jì)算,不斷迭代,直至相等或兩次閾值的差值小于設(shè)定值時(shí),最后一次迭代的閾值即為所求閾值。根據(jù)閾值,將原圖像轉(zhuǎn)化為二值圖像。所得二值圖像即為洪水區(qū)域的提取結(jié)果:

(16)

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

2.3 形態(tài)學(xué)濾波與洪水區(qū)域提取

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

洪水區(qū)域提取的最后一步,對(duì)災(zāi)前、災(zāi)后圖像進(jìn)行變化檢測(cè)。在之前的處理中,實(shí)現(xiàn)了對(duì)災(zāi)前與災(zāi)后水體的分別提取。在此基礎(chǔ)上,利用相減變化檢測(cè)可得到災(zāi)后相比災(zāi)前擴(kuò)張的水體部分,此部分即為洪災(zāi)區(qū)域。通過先提取水體再進(jìn)行變化檢測(cè)的方式,可抑制零星噪點(diǎn)對(duì)變化檢測(cè)結(jié)果的影響。具體流程為:利用災(zāi)后水體提取圖像減去災(zāi)前水體提取圖像,對(duì)應(yīng)相隨點(diǎn)值相減。對(duì)于未變化地區(qū),相減后的結(jié)果為0,對(duì)于遭受洪澇災(zāi)害的地區(qū),相減后結(jié)果為255。所得二值圖像即為遭受洪災(zāi)區(qū)域的提取結(jié)果。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

2020年7月,鄱陽(yáng)湖流域發(fā)生特大洪水,給當(dāng)?shù)厝罕姷纳?cái)產(chǎn)安全帶來了嚴(yán)重威脅。本文基于星載SAR圖像對(duì)2020年7月鄱陽(yáng)湖流域的洪水進(jìn)行研究。本文選用歐空局哨兵-1衛(wèi)星的數(shù)據(jù)進(jìn)行相關(guān)研究。選用的SAR圖像為L(zhǎng)evel-1級(jí)別,干涉寬幅(interferometric wide swath, IW)模式、地距多視影像(ground range detected, GRD)類型的單極化成像。本次使用的數(shù)據(jù)信息如表1所示。

表1 原始圖像數(shù)據(jù)信息Table 1 Original image data information

3.2 實(shí)驗(yàn)結(jié)果

經(jīng)過預(yù)處理的圖像如圖7所示。

圖7 預(yù)處理后的SAR圖像Fig.7 Preprocessed SAR image

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

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

由圖8可知,數(shù)據(jù)1與數(shù)據(jù)3應(yīng)采用六階高斯函數(shù),數(shù)據(jù)2采用五階高斯函數(shù)。高斯擬合后的結(jié)果如圖9所示。常量取值如表2所示。

圖9 原始數(shù)據(jù)與擬合結(jié)果Fig.9 Original data and fitting results

表2 高斯擬合函數(shù)參數(shù)Table 2 Parameters of Gaussian fitting function

利用迭代閾值分割算法分別計(jì)算3幅圖像的閾值,計(jì)算結(jié)果如下:2019年7月20日的SAR圖像閾值等級(jí)為119,圖像中對(duì)應(yīng)像素點(diǎn)的值為-16.18 dB;2020年3月16日的SAR圖像閾值等級(jí)為109,圖像中對(duì)應(yīng)像素點(diǎn)的值為-16.07 dB;2020年7月14日的SAR圖像閾值等級(jí)為122,圖像中對(duì)應(yīng)像素點(diǎn)的值為-15.36 dB。依據(jù)閾值對(duì)圖像進(jìn)行分割,得到初步提取的結(jié)果,再對(duì)其進(jìn)行形態(tài)學(xué)濾波,濾波前后結(jié)果如圖10所示。

圖10 形態(tài)學(xué)濾波前后水體提取結(jié)果Fig.10 Water area extraction before and after morphological filtering

對(duì)提取水體后的結(jié)果進(jìn)行統(tǒng)計(jì):數(shù)據(jù)1的SAR圖像中,提取水體所占比例為23.31%;數(shù)據(jù)2的SAR圖像中,提取水體所占比例為18.45%;數(shù)據(jù)3的SAR圖像中,提取水體所占比例為27.81%。與數(shù)據(jù)2相比,鄱陽(yáng)湖流域的水體面積在7月增加了50.73%;即使與數(shù)據(jù)3,即2019年同期相比,鄱陽(yáng)湖流域的水體面積也增加了19.30%。

分別對(duì)數(shù)據(jù)3與數(shù)據(jù)2、數(shù)據(jù)3與數(shù)據(jù)1的提取結(jié)果進(jìn)行相減變化檢測(cè),并將結(jié)果與災(zāi)前、災(zāi)后疊加形成的RGB (red green blue)圖像進(jìn)行對(duì)比,其結(jié)果如圖11所示。

圖11 洪災(zāi)區(qū)域提取結(jié)果對(duì)比Fig.11 Comparison of extraction of flood area

在RGB圖像中,紅色區(qū)域即為遭受洪澇災(zāi)害地區(qū),由圖11可知,通過水體提取和相減變化檢測(cè),洪災(zāi)區(qū)域基本實(shí)現(xiàn)了提取。根據(jù)相減變化檢測(cè)結(jié)果可知,與2020年3月相比,在2020年7月,鄱陽(yáng)湖流域的水體面積有明顯的增大,鄱陽(yáng)湖流域的西南部分?jǐn)U張尤為明顯,鄱陽(yáng)湖的北部區(qū)域擴(kuò)張相對(duì)較少;與2019年同期相比,鄱陽(yáng)湖也擴(kuò)張了一部分,擴(kuò)張部分集中在鄱陽(yáng)湖西部區(qū)域和東南區(qū)域。

為評(píng)估水體提取結(jié)果,利用Landsat-8光學(xué)衛(wèi)星開展驗(yàn)證實(shí)驗(yàn)。由于2020年7月鄱陽(yáng)湖流域雨水較多、云層較厚,光學(xué)衛(wèi)星難以觀測(cè),因此選用2019年7月Landsat-8的圖像與提取結(jié)果進(jìn)行對(duì)比,如圖12所示。

圖12 Landsat-8光學(xué)衛(wèi)星2019年7月鄱陽(yáng)湖區(qū)域圖像Fig.12 Landsat-8 optical satellite image of Poyang Lake area in July 2019

除了與迭代閾值分割算法進(jìn)行對(duì)比外,本文還選用Otsu閾值分割算法、K-Means非監(jiān)督性學(xué)習(xí)、EM算法的提取結(jié)果進(jìn)行對(duì)比,如圖13所示。

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

為驗(yàn)證本文所提方法的有效性,通過采樣評(píng)估對(duì)上述3種提取方式的提取效果進(jìn)行評(píng)估。具體方法為:在光學(xué)衛(wèi)星圖像中分別隨機(jī)選取300個(gè)點(diǎn),分別統(tǒng)計(jì)3種提取方式的準(zhǔn)確率、虛警率、漏警率、計(jì)算時(shí)間等。其結(jié)果如表3所示。

表3 不同提取方式評(píng)估結(jié)果Table 3 Evaluation of different extraction methods

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

4 結(jié) 論

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

后續(xù)在這方面的主要工作內(nèi)容為:① 對(duì)算法的適用性進(jìn)行進(jìn)一步檢驗(yàn),考慮復(fù)雜情況下的洪澇區(qū)域提取,如城市內(nèi)澇等;② 考慮結(jié)合光學(xué)圖像信息,改進(jìn)算法結(jié)構(gòu),提升準(zhǔn)確率。