張瀟予,趙鳳軍,李 寧

(1.中國科學(xué)院電子學(xué)研究所,北京100190;2.中國科學(xué)院大學(xué),北京100190)

0 引言

圖像變化檢測是研究不同時段同一場景圖像之間發(fā)生的變化[1-2],在災(zāi)害評估、土地利用、城市規(guī)劃及軍事偵察等方面具有廣闊的應(yīng)用前景。合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)是一種主動式微波成像傳感器[3-6],與傳統(tǒng)光學(xué)遙感和高光譜遙感相比,SAR具有全天時、全天候、大區(qū)域?qū)Φ赜^測的特點(diǎn),因而在變化檢測領(lǐng)域的應(yīng)用越來越廣。

在SAR圖像變化檢測技術(shù)中,無監(jiān)督變化檢測由于不需要先驗變化信息的支持,成為國內(nèi)外研究的重點(diǎn)[7]。目前,SAR圖像變化檢測主要有兩類思想:1)分類后比較法,即先將兩幅SAR圖像分別進(jìn)行分類,然后比較分類結(jié)果,從而得到變化信息。2)先獲取差異圖,然后對差異圖進(jìn)行分析,從而獲取最終的變化信息。在這兩類思想中,第一類思想的變化檢測精度受分類結(jié)果的影響很大,且容易導(dǎo)致分類誤差積累問題,而第二類思想的思路簡單明了,且具有較大的研究空間,被廣泛采用。但第二類思想仍存在以下3個問題:1)如何抑制斑點(diǎn)噪聲。SAR是相干系統(tǒng),斑點(diǎn)噪聲是其固有特性,SAR圖像中斑點(diǎn)噪聲的存在無疑會影響變化檢測的結(jié)果。SAR圖像的這一特性使SAR圖像變化檢測成為變化檢測領(lǐng)域極具挑戰(zhàn)的難點(diǎn)和熱點(diǎn)問題。2)如何構(gòu)造性能良好的差異圖。差異圖可分離性的好壞直接影響變化檢測的最終結(jié)果,因此構(gòu)造性能良好的差異圖成為這類方法的研究重點(diǎn)之一。3)如何設(shè)計一種有效的分類方法,用以區(qū)分變化信息和非變化信息。這是SAR圖像變化檢測中的重要步驟之一,會直接影響變化檢測的結(jié)果。

傳統(tǒng)的SAR圖像變化檢測方法僅使用單一特征,未能充分利用圖像的細(xì)節(jié)信息,因此檢測效果并不理想。近年來,多通道特征如ICF(Integral Channel Features),ACF(Aggregate Channel Features)等被成功應(yīng)用于行人檢測[8]及人臉檢測[9-10]領(lǐng)域。多通道特征從不同角度集合了多種通道特征,圖像信息更加豐富。基于上述內(nèi)容,本文提出一種基于多通道特征的無監(jiān)督SAR圖像變化檢測方法。首先,引入引導(dǎo)濾波來抑制相干斑噪聲并盡可能地保留SAR圖像的邊緣及局部信息;其次,提取8個通道特征來獲得性能良好的差異圖;最后,利用主成分分析(PCA)和K-means聚類得到最終的變化信息。

1 本文算法

本文算法可以分為3個部分:引導(dǎo)濾波、基于多通道特征的差異圖獲取和差異圖分析。算法整體流程如圖1所示。

1.1 引導(dǎo)圖像濾波

為了在抑制相干斑噪聲的同時盡可能多地保留SAR圖像的邊緣及局部信息,本文引入引導(dǎo)圖像濾波[11]對SAR圖像進(jìn)行預(yù)處理。

引導(dǎo)圖像濾波的本質(zhì)是將空域濾波模型推廣到基于圖像信息的濾波,它包括引導(dǎo)圖像I、輸入圖像p、輸出圖像q,其中I和p是根據(jù)具體應(yīng)用事先就給定的。

假設(shè)輸出圖像q與引導(dǎo)圖像I在以k為窗后中心的窗口wk中滿足線性關(guān)系,如下所示:

圖1 本文算法流程圖

式中,(ak,bk)是當(dāng)窗口為wk時的線性系數(shù)。假定窗口大小為r,這一個線性的關(guān)系就保證了只有當(dāng)I有邊緣時,輸出q才有邊緣,這是因為?q=a?i。這一模型在圖像摳圖[12]、超分辨率圖像[13]、圖像去霧[14]等應(yīng)用中得到有效的證實。

為了確定線性系數(shù)的大小,使得輸出q和輸入p之間的差異最小,即將下式最小化:

式中,ε是用來防止ak過大的正則化因子,則上式得解為

式中,uk和分別為I在窗口wk下的均值和方 差,|w|為窗口中像素的個數(shù),接下來計算整幅圖像每個窗口的線性系數(shù),但是由于一個像素會被多個窗口包含,使得式(1)中qi的值并不一樣,因此在計算每一個輸出時可以采用計算包含該點(diǎn)的所有線性系數(shù)的平均值。計算公式如下:

實際上,可以將式(3)重新寫成p的加權(quán)和:ak=∑j Akj(I)pj,其中Aij是僅依賴于I的權(quán)重,同理:bk=∑j Bkj(I)pj,qi=∑j Wij(I)pj。當(dāng)把引導(dǎo)濾波用作邊緣保持濾波器時,往往有I=p,如果ε=0,顯然ak=1,bk=0是E(ak,bk)為最小值的解,這時的濾波器沒有任何作用,輸出原始圖像。如果ε＞0,在像素強(qiáng)度變化小的區(qū)域(或單色區(qū)域),有a近似于(或等于)0,而b近似于(或等于),相當(dāng)于對圖像作了一個加權(quán)均值濾波;而在像素強(qiáng)度變化大的區(qū)域,a近似于1,b近似于0,對圖像的濾波效果很弱,有助于保持邊緣。ε的作用是決定濾波效果,在窗口大小不變的情況下,隨著ε的增大,濾波效果越明顯。

綜上所述,引導(dǎo)濾波具有與雙邊濾波那樣保持邊緣平滑的特性,不受梯度反轉(zhuǎn)偽影影響,所以它在細(xì)節(jié)處理上比雙邊濾波具有更好的效果,被廣泛應(yīng)用于光學(xué)圖像的去噪、細(xì)節(jié)平滑、圖像增強(qiáng)和圖像去霧等領(lǐng)域。引導(dǎo)濾波最大的優(yōu)勢在于,可以寫出時間復(fù)雜度與窗口大小無關(guān)的算法,因此在使用大窗口處理圖片時,其效率更高。由于其良好的邊緣保持特性和較快的計算效率(時間復(fù)雜度為O(N)),本文將其應(yīng)用到SAR圖像去噪,并取得了良好的效果。

1.2 基于多通道特征的差異圖獲取

在變化檢測中,差異圖獲取是一個關(guān)鍵環(huán)節(jié),差異圖的好壞影響最終的變化檢測結(jié)果。為獲取性能良好的差異圖,本文提出了基于多通道特征的差異圖獲取。首先對預(yù)處理后的SAR圖像提取8個通道特征,然后利用對數(shù)比值法獲取8個通道特征的差異圖,最后對各通道差異圖進(jìn)行融合,得到最終的差異圖。采用多個通道特征獲取差異圖,可以從不同角度集成多種特征,圖像信息更豐富,更有利于區(qū)分變化信息和非變化信息。

1.2.1 通道特征

圖像的通道廣義上可以看作是原始圖像的某種映射,它的像素值是由原始圖像的對應(yīng)區(qū)域通過某種計算得到的[8]。給定一個輸入圖像I,其通道圖像C的計算公式如下:

式中,Ω代表圖像的某種通道計算函數(shù)。

通道類型有多種,如顏色通道、灰度、線性濾波、非線性變換、逐點(diǎn)變換、積分直方圖、梯度直方圖等[8]。本文選取8個通道特征用于獲取差異圖,包括1個灰度通道、1個歸一化的梯度幅值通道和6個梯度方向直方圖通道。在選取的8個通道特征中:

1)灰度是最簡單的通道。對于SAR圖像而言,由于其本身就是一張灰度圖,所以,SAR圖像的灰度通道C=I,即SAR圖像的原圖,如圖2(a)所示。

2)梯度幅值通道屬于非線性變換通道的一種,如圖2(b)所示,梯度幅值通道可以捕獲邊緣強(qiáng)度信息。

3)梯度方向直方圖是一個加權(quán)直方圖,如圖2(c)～(h)所示。梯度方向取0°～360°,然后將梯度分布平均分成6個方向角度,即每60°取一個方向角度。梯度方向直方圖計算公式如下:

式中,G(x,y)代表圖像I(x,y)處的梯度幅值,Θ(x,y)代表圖像I(x,y)處量化的梯度方向。通過對原始SAR圖像I進(jìn)行不同尺度的模糊,可以計算出不同尺度的梯度信息。此外,借助于梯度幅值信息,可以對計算出來的梯度直方圖進(jìn)行歸一化。

本文所使用的多通道特征提取流程如圖2所示。

從圖2可以看出,灰度通道和梯度通道的大小和原圖相同,而6個梯度方向直方圖的大小僅為原圖的1/3,這是由于在計算梯度方向直方圖時使用了3×3的單元格。為了得到和原圖大小相同的梯度方向直方圖,方便后期進(jìn)行差異圖融合,本文對其進(jìn)行插值處理。從而最終提取到所使用的8個通道特征。

之所以選取這8個通道特征,是因為它們不僅提取方法操作簡單,可以快速高效地計算出來,而且包含豐富的圖像細(xì)節(jié)信息,可以從不同角度收集輸入圖像的信息,有利于抑制斑點(diǎn)噪聲,獲得性能良好的差異圖,從而提高檢測精度。

圖2 多通道特征提取流程圖

1.2.2 差異圖融合

對數(shù)比值法在獲取差異圖方面有獨(dú)特的優(yōu)勢,它不僅能夠很大程度上消除乘性噪聲的影響,減少定標(biāo)處理引起的額外誤差,而且能夠突顯出SAR圖像的相對變化區(qū)域[15]。為充分利用不同的通道特征信息,生成基于多通道特征的融合差異圖,本文利用對數(shù)比值法獲取每一個通道特征的差異圖,計算公式如下:

式中,Ri表示各通道特征的差異圖,I0i表示t0時刻第i個通道圖像的幅度,I1i表示t1時刻第i個通道圖像的幅度。

在獲取各通道特征的差異圖后,取各通道差異圖的平均值作為最終的差異圖,計算公式如下:

1.3 差異圖分析

在獲取性能良好的差異圖后,為了有效地區(qū)分變化信息和非變化信息,得到精確、可靠的變化檢測結(jié)果,本文采用PCA和K-means聚類對差異圖進(jìn)行分析[16]。具體步驟如下:

1)將差異圖R分成n×n的互不重疊的塊(假設(shè)分成了m個小塊),將每一小塊轉(zhuǎn)化為行向量Xi(i=1,2,…,m),用這些行向量構(gòu)建矩陣F,矩陣大小為m×n2。分塊的思想使得算法具有一定的抗噪性能。

2)用PCA將步驟1)中矩陣F的行向量維度從n2維降成S維,得到對應(yīng)的變換矩陣Z。

3)取差異圖R中每個像素所在的n×n鄰域小塊,將每個小塊轉(zhuǎn)化成行向量,共有H×W個行向量(H×W為差異圖R的大小)。用得到的變換矩陣Z將所有的行向量變換到特征向量空間Q上(每個像素就由S維的向量表示),Q的維度為S×HW。

4)用K-means聚類算法將特征向量空間Q聚為兩類:變化類和非變化類,從而得到最終的變化檢測結(jié)果。

2 實驗結(jié)果

2.1 實驗數(shù)據(jù)集

本文選取兩組配準(zhǔn)后的數(shù)據(jù)集進(jìn)行算法驗證及實驗分析。

數(shù)據(jù)集I是由RADARSAT-1衛(wèi)星拍攝的加拿大渥太華地區(qū)的兩景SAR圖像(分別如圖3(a)、圖3(b)所示),拍攝時間分別為1997年7月和1997年8月,圖像大小為290×350,圖3(c)為變化信息參考圖。

數(shù)據(jù)集II是由C波段ERS-2衛(wèi)星拍攝的伯爾尼地區(qū)的兩景SAR圖像(分別如圖4(a)、圖4(b)所示),拍攝時間分別為1999年4月和1999年5月,圖像大小為301×301,圖4(c)為變化信息參考圖。

圖3 數(shù)據(jù)集Ⅰ:渥太華地區(qū)SAR圖像

圖4 數(shù)據(jù)集Ⅱ:伯爾尼地區(qū)SAR圖像

之所以選取這兩組數(shù)據(jù)集,首先是因為這兩組數(shù)據(jù)集都存在較多的斑點(diǎn)噪聲,其次是因為這兩組數(shù)據(jù)集具有對比性,在數(shù)據(jù)集I的圖像中有較多的變化像素,而數(shù)據(jù)集II的圖像中變化像素較少。

2.2 評價指標(biāo)

在評價指標(biāo)方面,本文選取錯檢數(shù)(FP)、漏檢數(shù)(FN)、總體錯誤數(shù)(OE)、百分比修正(PCC)、Kappa系數(shù)這5個檢測指標(biāo)來評價算法性能,其中:

1)FP:未變化像素被錯檢為變化像素的個數(shù),即錯檢數(shù)(或稱虛警數(shù))。

2)FN:變化像素被檢測為未變化像素的個數(shù),即漏檢數(shù)。

3)OE:總體錯誤數(shù),其值為錯檢數(shù)與漏檢數(shù)之和。

4)PCC:百分比修正,計算公式如下:

式中,TP為正確檢測到的變化像素的數(shù)目,TN為正確檢測到的未變化像素的數(shù)目。PCC的值越高,說明算法檢測效果越好。

5)Kappa系數(shù):計算公式如下:

式中,

式中,Mc為實際的變化像素的數(shù)目,Mu為實際的非變化像素的數(shù)目。

Kappa系數(shù)的值越高,說明算法的檢測精度越高。

2.3 實驗結(jié)果及分析

2.3.1 檢測結(jié)果分析

為驗證算法的性能,本文算法GF_MCF將與以下3種算法進(jìn)行對比:文獻(xiàn)[16]中的算法命名為PCAKM;文獻(xiàn)[17]中的算法命名為MRFFCM;文獻(xiàn)[18]中的算法命名為MOHO。

兩組數(shù)據(jù)集的變化檢測結(jié)果如圖5、圖6所示。其中,圖5(a)和圖6(a)顯示的是不同算法提取的變化信息。圖5(b)和圖6(b)展示了更具體的檢測結(jié)果。具體評價指標(biāo)如表1所示。

圖5 渥太華地區(qū)不同算法檢測結(jié)果

圖6 伯爾尼地區(qū)不同算法的檢測結(jié)果

表1 算法性能指標(biāo)

1)數(shù)據(jù)集I(渥太華地區(qū))實驗結(jié)果分析

從圖5可以看出,本文算法GF_MCF在對數(shù)據(jù)集I的檢測效果上要明顯優(yōu)于另外3種算法,不僅有較少的漏檢和錯檢像素,且保持了更多的圖像細(xì)節(jié)。這是由于GF_MCF算法結(jié)合了多個通道特征的優(yōu)點(diǎn),圖像細(xì)節(jié)信息更加豐富,因此檢測效果更好。從表1可以看出,GF_MCF算法在OE,PCC和Kappa系數(shù)上要明顯優(yōu)于另外3種算法,進(jìn)一步說明了GF_MCF算法的檢測優(yōu)勢。

2)數(shù)據(jù)集II(伯爾尼地區(qū))實驗結(jié)果分析

數(shù)據(jù)集II與數(shù)據(jù)集I的特點(diǎn)不同,數(shù)據(jù)集II的變化區(qū)域較小且比較集中,圖像細(xì)節(jié)更容易被錯檢和漏檢。從圖6可以看出,由于斑點(diǎn)噪聲的影響,PCAKM,MRFFCM和MOHO這3種算法的檢測結(jié)果中都出現(xiàn)了孤立點(diǎn),而GF_MCF算法的檢測結(jié)果無孤立點(diǎn)產(chǎn)生,對斑點(diǎn)噪聲的抑制明顯強(qiáng)于另外3種算法。這不僅是因為GF_MCF算法中使用了引導(dǎo)濾波,還因為多通道特征的思想對斑點(diǎn)噪聲也有一定的抑制作用,從而提高了算法的抗噪性能。同時,從表1可以看出,GF_MCF的FN和OE最低,PCC和Kappa系數(shù)最高,檢測效果最好。

3)通道特征分析

本文采用了3種類型的通道特征:灰度通道、梯度幅值通道和梯度直方圖通道(6個)。圖7展示了不同通道特征之間相互組合后對檢測結(jié)果的影響,其中n代表差異圖分析部分中的窗口大小。從圖7可以看出,結(jié)合了3種類型通道特征的檢測效果最好。此外,從圖中還可以看出,當(dāng)窗口n=3時,PCC的值最高,檢測效果最好。此組實驗再次證明了通道特征在SAR圖像變化檢測領(lǐng)域的有效性。

圖7 通道特征組合分析

2.3.2 引導(dǎo)圖像濾波性能分析

為了驗證引導(dǎo)圖像濾波在SAR圖像變化檢測領(lǐng)域的適用性,本文在PCAKM,MRFFCM和MOHO算法中加入了引導(dǎo)濾波,分別命名為GF_PCAKM,GF_MRFFCM和GF_MOHO,檢測性能指標(biāo)如表1所示。從表1可以看出,加入引導(dǎo)濾波后,PCAKM,MRFFCM和MOHO算法的檢測性能都有明顯的提升。

同時,為進(jìn)一步說明引導(dǎo)濾波的性能及有效性,本文將引導(dǎo)圖像濾波與其他濾波算法在檢測性能和濾波時間等指標(biāo)方面進(jìn)行對比,具體結(jié)果如表2所示。

表2 各濾波算法性能指標(biāo)

從表2可以看出,引導(dǎo)濾波和其他濾波算法相比,OE最低,PCC和Kappa系數(shù)最高,并且濾波所耗時間最短,說明引導(dǎo)圖像濾波的檢測性能明顯優(yōu)于另外幾種濾波算法。

在引導(dǎo)濾波中,可調(diào)參數(shù)有2個:窗口大小r和正則化因子ε。圖8展示了引導(dǎo)濾波的2個參數(shù)對檢測效果的影響。從圖8可以看出,當(dāng)窗口大小r=13,正則化因子ε=0.007時,PCC最高,檢測效果最好。

圖8 引導(dǎo)圖像濾波參數(shù)對檢測效果的影響

綜上所述,引導(dǎo)圖像濾波對斑點(diǎn)噪聲有一定的抑制作用,適用于SAR圖像變化檢測,且檢測效果較好。

2.3.3 算法抗噪性能分析

以上兩部分的實驗分析表明,本文算法具有良好的抗噪性能。為了進(jìn)一步驗證本文算法GF_MCF對斑點(diǎn)噪聲的抑制能力,本文給出了不同算法的抗噪性能曲線,如圖9所示。

給定一個輸入SAR圖像X,大小為H×W,加入斑點(diǎn)噪聲后的圖像計為。加入的噪聲水平用峰值信噪比(PSNR)來衡量,PSNR計算公式如下(假設(shè)輸入SAR圖像的最大像素值為1):

給定兩個SAR圖像X1,X2,用變化檢測算法得到檢測結(jié)果圖CM1,然后往X1中加入斑點(diǎn)噪聲,X2保持不變,用同一種變化檢測算法得到檢測結(jié)果圖CM2。CM1和CM2之間的差異可用參數(shù)τ表示,計算公式如下:

τ的大小代表算法對斑點(diǎn)噪聲的抑制能力,其值越高則算法的抗噪性能越好。

圖9為不同算法的抗噪性能曲線。從圖9可以看出,本文算法GF_MCF對斑點(diǎn)噪聲具有較好的抑制能力,這也是本文算法的優(yōu)勢之一。

圖9 算法抗噪性能曲線

3 結(jié)束語

本文提出了一種基于多通道特征的SAR圖像變化檢測方法。實驗結(jié)果表明,該方法對SAR圖像的變化具有良好的檢測效果,檢測精度進(jìn)一步提高且具有較好的抗噪性能。在未來的工作中,擬進(jìn)一步對通道特征進(jìn)行探索,并考慮用兩級聚類分析差異圖,以得到更為精確的檢測結(jié)果并進(jìn)一步提高算法的魯棒性。

[1]GONG Maoguo,LI Yu,JIAO Licheng,et al.SAR Change Detection Based on Intensity and Texture Changes[J].ISPRS Journal of Photogrammetry and Remote Sensing,2014,93(1):123-135.

[2]LI H,CELIK T,LONGBOTHAM N,et al.Gabor Feature Based Unsupervised Change Detection of Multitemporal SAR Images Based on Two-Level Clustering[J].IEEE Geoscience and Remote Sensing Letters,2015,12(12):2458-2462.

[3]鄧云凱,趙鳳軍,王宇.星載SAR技術(shù)的發(fā)展趨勢及應(yīng)用淺析[J].雷達(dá)學(xué)報,2012,1(1):1-10.

[4]NAIDOO L,MATHIEU R,MAIN R,et al.The Assessment of Data Mining Algorithms for Modelling Savannah Woody Cover Using Multi-Frequency(X-,C-and L-Band)Synthetic Aperture Radar(SAR)Datasets[C]∥2014 IEEE International Geoscience and Remote Sensing Symposium,Quebec City,QC,Canada:IEEE,2014:1049-1052.

[5]BHATTACHARYA C,ROY A,NAVNEET S,et al.MicroSAR:Calibration of X-Band High Resolution FMCW Synthetic Aperture Radar(SAR)[C]∥2014 IEEE International Microwave and RF Conference,Bangalore,India:IEEE,2014:377-380.

[6]KRIEGER G,MOREIRA A,FIEDLER H,et al.TanDEM-X:A Satellite Formation for High-Resolution SAR Interferometry[J].IEEE Trans on Geoscience and Remote Sensing,2007,45(11):3317-3341.

[7]WANG Yan,DU Lan,DAI Hui.Unsupervised SAR Image Change Detection Based on SIFT Keypoints and Region Information[J].IEEE Geoscience and Remote Sensing Letters,2016,13(7):931-935.

[8]DOLLáR P,TU Z,PERONA P,et al.Integral Channel Features[C]∥British Machine Vision Conference,London:BMVA Press,2009:1-11.

[9]YANG Bin,YAN Junjie,LEI Zhen,et al.Aggregate Channel Features for Multi-View Face Detection[C]∥2014 IEEE International Joint Conference on Biometrics,Clearwater,FL,USA:IEEE,2014:1-8.

[10]WANG Shuo,YANG Bin,LEI Zhen,et al.A Convolutional Neural Network Combined with Aggregate Channel Feature for Face Detection[C]∥6th International Conference on Wireless,Mobile and Multi-Media,Beijing:IET,2015:304-308.

[11]HE Kaiming,SUN Jian,TANG Xiaoou.Guided Image Filtering[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2013,35(6):1397-1409.

[12]LEVIN A,LISCHINSKI D,WEISS Y.A Closed-Form Solution to Natural Image Matting[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2008,30(2):228-242.

[13]ZOMET A,PELEG S.Multi-Sensor Super-Resolution[C]∥Sixth IEEE Workshop on Applications of Computer Vision,Orlando,FL,USA:IEEE,2002:27-31.

[14]HE Kaiming,SUN Jian,TANG Xiaoou.Single Image Haze Removal Using Dark Channel Prior[J].IEEE Trans on Pattern Analysis and Machine Intelligence,2011,33(12):2341-2353.

[15]吳濤,陳曦,牛蕾,等.非監(jiān)督SAR圖像變化檢測研究最新進(jìn)展[J].遙感信息,2013,28(1):110-118.

[16]CELIK T.Unsupervised Change Detection in Satellite Images Using Principal Component Analysis and k-Means Clustering[J].IEEE Geoscience and Remote Sensing Letters,2009,6(4):772-776.

[17]GONG Maoguo,SU Linzhi,JIA Meng,et al.Fuzzy Clustering with a Modified MRF Energy Function for Change Detection in Synthetic Aperture Radar Images[J].IEEE Trans on Fuzzy Systems,2014,22(1):98-109.

[18]萬紅林,汪洋,江凱,等.基于模糊貼近度和非緊湊鄰域的變化檢測[J].雷達(dá)科學(xué)與技術(shù),2014,12(3):229-234.WAN Honglin,WANG Yang,JIANG Kai,et al.Fuzzy Neartude and Information in a Non-Compact Neighborhood Based Change Detection[J].Radar Science and Technology,2014,12(3):229-234.(in Chinese)