利用紋理融合與廣義高斯模型的高分辨率SAR影像變化檢測

劉本強(qiáng)，趙 爭，盛玉婷，張忠芳

(1.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266590；2.中國測繪科學(xué)研究院，北京 100830；3.北京航天世景信息技術(shù)有限公司，北京 100089)

合成孔徑雷達(dá)(Synthetic Aperture Radar, SAR)具有全天時、全天候的監(jiān)測能力。利用SAR影像進(jìn)行變化檢測是當(dāng)前研究的熱點(diǎn)問題。變化檢測是通過對同一地區(qū)不同時期的遙感影像比對分析，獲取地物的變化信息，在災(zāi)害監(jiān)測[1-2]、土地利用[3]、城市化發(fā)展[4]等方面應(yīng)用廣泛。

近年來，Cosmo-SkyMed、TerraSAR-X、Radarsat-2和ALOS-2等商業(yè)合成孔徑雷達(dá)衛(wèi)星的發(fā)射，極大滿足獲取高分辨率SAR影像的需求。不同于中低分辨率SAR圖像，高分辨率的SAR圖像具有更加豐富的紋理、幾何和細(xì)節(jié)等特征，得到研究人員的廣泛關(guān)注[5-7]。文獻(xiàn)[8-9]針對高分辨率的機(jī)載SAR數(shù)據(jù)，綜合考慮灰度共生矩陣紋理、LBP紋理、Gabor紋理等紋理特征，描述震后建筑物的損壞等級。在變化檢測中引入紋理信息不僅能夠抑制斑點(diǎn)噪聲，而且在檢測性能上也會有所提高。陳志鵬等[10]采用紋理信息代替灰度信息，使用差值法獲取城區(qū)變化檢測結(jié)果，并比較不同紋理的檢測性能。李坤等[11]針對地震災(zāi)害導(dǎo)致的地表變化問題，分析SAR圖像灰度和紋理特征的敏感程度，結(jié)合主成分分析和K均值聚類技術(shù)提取震前震后的變化信息。Wang等[12]利用方向?yàn)V波和分形紋理進(jìn)行SAR影像變化檢測。趙凌君等[13]構(gòu)建9種常見的紋理特征，然后利用Bhattacharyya距離進(jìn)行特征選擇，通過K近鄰分類器實(shí)現(xiàn)建筑區(qū)的檢測。

廣義高斯模型(Generalized Gaussian Model，GGM)能夠較好的描述SAR圖像統(tǒng)計特性，文獻(xiàn)[14-16]均利用廣義高斯模型描述差異圖的變化類像元和未變化類像元的概率分布，并基于KI閾值法實(shí)現(xiàn)多時相單極化SAR影像的非監(jiān)督變化檢測。Zhao等[17-18]結(jié)合統(tǒng)計模型的方法分別利用Kittler-Illingworth(KI)算法和改進(jìn)的KI閾值算法對相應(yīng)的差異影像進(jìn)行閾值提取，實(shí)現(xiàn)多時相全極化SAR數(shù)據(jù)的非監(jiān)督變化檢測，獲得較好的變化檢測結(jié)果。

本文為了充分利用高分辨率SAR影像的灰度信息和紋理信息，提高變化檢測精度，提出一種紋理信息融合與廣義高斯模型相結(jié)合的SAR影像變化檢測方法。該方法以灰度共生矩陣(gray level co-occurrence matrix, GLCM)計算的紋理特征構(gòu)造紋理差值差異圖，以對數(shù)比值運(yùn)算構(gòu)造灰度差異圖，以小波變換為工具融合兩類差異圖。然后通過廣義高斯模型對融合后的差異圖統(tǒng)計建模，估計變化類像元和未變化類像元的概率分布，利用KI閾值準(zhǔn)則獲取最佳分割閾值，得到最終變化檢測結(jié)果。

1 算法介紹

1.1 灰度共生矩陣與紋理特征

紋理描述像素的鄰域灰度空間分布規(guī)律，是特征提取、影像分類、目標(biāo)識別的重要手段。常用的紋理提取方法有統(tǒng)計方法、幾何方法、結(jié)構(gòu)方法、模型方法和信號處理方法。灰度共生矩陣是最常用的一種基于統(tǒng)計的紋理提取方法，由Haralick等[19]提出，定義14種基于GLCM的紋理特征。

灰度共生矩陣定義為

P(i,j,d,θ)=#{[(x,y),(x+Δx,y+Δy)|d,

θ,f(x,y)=i,f(x+Δx,y+Δy)=j]}.

(1)

歸一化方式為

(2)

式中：#{?}表示大括號中成立的像元對數(shù)；(x,y)為像元坐標(biāo)，且x=0,1,…,M-1,y=0,1,…,N-1;M，N為影像的高度和寬度；i,j=0,1,…,L-1表示灰度級；Δx，Δy為像素坐標(biāo)偏移量；d為步長，θ為方向，一般取0°、45°、90°、135°；f(?)表示像素坐標(biāo)的灰度值。

根據(jù)灰度共生矩陣可以定義如下紋理特征：角二階矩(angular second moment, ASM)、熵(entropy)、逆差矩(inverse difference moment, IDM)。

角二階矩：

(3)

熵：

(4)

逆差矩：

(5)

1.2 基于小波變換的差異圖融合

小波變換作為一種數(shù)學(xué)工具，在圖像融合處理領(lǐng)域有著很重要的地位。Gong等[20]和慕彩虹等[21]分別利用離散小波變換(discrete wavelet transform, DWT)和離散平穩(wěn)小波變換(stationary wavelet transform, SWT)將對數(shù)比值差異圖和均值比差異圖融合，取得顯著的效果。

首先利用式(6)和式(7)構(gòu)造紋理差異圖和灰度差異圖，其中Ti(x,y)為i時刻SAR影像的紋理特征影像，Ii(x,y)為i時刻SAR強(qiáng)度影像，i=1,2表示兩個時相，x和y表示影像的行列。

TDI(x,y)=|T1(x,y)-T2(x,y)|,

(6)

LR(x,y)=log(I1(x,y)/I2(x,y)).

(7)

然后，對兩種差異圖分別進(jìn)行3層SWT小波分解，每一層分解分別得到一個低頻系數(shù)矩陣和3個高頻系數(shù)矩陣。對于低頻系數(shù)，采用式(8)的加權(quán)平均融合規(guī)則，對于高頻系數(shù)，采用式(9)的局部能量最小融合規(guī)則。

(8)

(9)

Ek(x,y)=∑Ω(Hk(x,y))2

(10)

表示以(x,y)為中心的局部能量，Ω為計算局部能量時窗口大小，本文中取3×3。

最后，將融合后的低頻系數(shù)和高頻系數(shù)執(zhí)行SWT逆變換，得到融合后的差異影像。

1.3 廣義高斯模型下的KI閾值(GGKI)

對于融合后差異圖Fswt(x,y)上存在著的變化類ωc和未變化類ωu，根據(jù)貝葉斯決策理論，差異圖上的任意一個像元Fswt(x,y)可由式(11)表示，其中Xf表示融合后差異圖的像元灰度值。

p(Xf)=P(ωu)p(Xf|ωu)+P(ωc)p(Xf|ωc).

(11)

在廣義高斯模型的假設(shè)下，兩類像元的條件概率密度可以表示為

p(Xf|ωi)=aiexp(-(bi|Xf|-mi)βi),

i∈{u,c},

(12)

其中

(13)

參數(shù)mi，σi和βi分別表示廣義高斯分布的均值、標(biāo)準(zhǔn)差和形狀控制參數(shù)，Γ(?)為Gamma函數(shù)，表達(dá)式為

(14)

KI閾值選擇法是一種以最小錯誤率貝葉斯理論為基礎(chǔ)的閾值選擇方法。圖像直方圖h(Xf)的最佳分割閾值T由一個能夠描述整幅圖像的準(zhǔn)則函數(shù)J(T)確定。

(15)

式中c(Xf,T)為代價函數(shù)，其表達(dá)式為

(16)

其中，p(ωi|Xf,T)(i=u,c)為閾值T下的未變化類和變化類的后驗(yàn)概率密度。使得圖像分類誤差最小的閾值為最佳閾值，即T*為最佳閾值時需使得準(zhǔn)則函數(shù)J(T)最小。

(17)

根據(jù)貝葉斯理論，后驗(yàn)概率密度函數(shù)可以由先驗(yàn)概率和條件概率密度表示為

p(ωi|Xf,T)=

(18)

綜上所述，可以推導(dǎo)出準(zhǔn)則函數(shù)的表達(dá)式為

H(Ω,T)-(Pu(T)lnau(T)+Pc(T)lnac(T)).

(19)

式中H(Ω,T)表示類別集合Ω∈{ωu,ωc}的熵，表達(dá)式為

(20)

對于廣義高斯模型中相關(guān)參數(shù)的估計，可以采用

(21)

2 實(shí)驗(yàn)與分析

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

實(shí)驗(yàn)選用3 m分辨率的TerraSAR-X數(shù)據(jù)，覆蓋廣州市和貴陽市某區(qū)域。為方便實(shí)驗(yàn)研究，本文選取其中1 500 m×1 300 m和1 000 m×1 000 m的區(qū)域進(jìn)行實(shí)驗(yàn)，見圖1和圖2。其中圖1(c)和圖2(c)為人工標(biāo)定的變化參考圖，用來定量評價檢測效果。

圖1 廣州市TerraSAR-X數(shù)據(jù)

圖2 貴陽市TerraSAR-X數(shù)據(jù)

2.2 實(shí)驗(yàn)方法

圖3是本文方法的流程圖。主要包括3個部分內(nèi)容：

1)對實(shí)驗(yàn)影像做預(yù)處理(輻射校正、幾何校正、影像配準(zhǔn)、影像濾波)，其中影像濾波采用3×3 Lee濾波，以削弱斑點(diǎn)噪聲的影響；

2)GLCM計算中，紋理窗口選擇11×11，d=1，θ分別取0°、45°、90°、135°，計算后的4個方向紋理特征取均值，再參照式(6)構(gòu)造不同紋理特征的紋理差異圖；

3)按照式(7)構(gòu)造對數(shù)比值差異圖，并與不同紋理特征的紋理差值差異圖做基于SWT變換的融合處理；

4)采用廣義高斯分布對融合后的差異圖建模，利用KI閾值準(zhǔn)則獲得最佳分割閾值，提取變化信息并評價檢測精度。

圖3 本文方法流程

2.3 實(shí)驗(yàn)結(jié)果及分析

將對數(shù)比值圖與紋理差值差異圖經(jīng)SWT變換，得到的融合差異圖見圖4和圖5，其分別是對數(shù)比值圖與角二階距紋理(ASM)、熵紋理(ENT)、逆差距紋理(IDM)的融合結(jié)果、對數(shù)比值差異圖?？梢钥闯?，在對數(shù)比值差異圖中，背景中的噪點(diǎn)特別多，變化區(qū)域也“破破碎碎”。而融合紋理信息的差異圖相對平滑，噪點(diǎn)明顯減少，且變化區(qū)域更加完整。這是因?yàn)樵朦c(diǎn)為頻率域的高頻部分，而高頻系數(shù)采用局部能量最小原則進(jìn)行融合，對噪點(diǎn)有一定的抑制作用。同時，每一層的SWT分解也會降低噪聲的影響。

采用廣義高斯分布對融合差異圖建模，擬合變化類像素和未變化類像素的概率分布，利用KI閾值準(zhǔn)則獲取最佳分割閾值提取變化信息，得到變化二值圖。圖6和圖7分別是廣義高斯模型對兩組數(shù)據(jù)集融合差異圖的擬合結(jié)果。可以看出對于融合后的差異圖，廣義高斯模型均能較好的擬合其變化類像元和未變化類像元的概率密度函數(shù)，特別注意的是，廣義高斯模型對于融合逆差矩(IDM)紋理的差異圖擬合效果最好。

圖4 廣州數(shù)據(jù)集差異圖

圖5 貴陽數(shù)據(jù)集差異圖

圖6 廣州數(shù)據(jù)集廣義高斯模型擬合結(jié)果

圖7 貴陽數(shù)據(jù)集廣義高斯模型擬合結(jié)果

基于廣義高斯模型KI閾值算法(GGKI)的變化檢測結(jié)果見圖8和圖9?？梢灾庇^地看出融合3種紋理信息均能夠較好的提取出變化區(qū)域，但融合ASM紋理檢測結(jié)果的背景中含有較多的噪點(diǎn)。

為了定量評價檢測精度，驗(yàn)證所提方法的有效性，分別對不同方法得到的結(jié)果計算檢測指標(biāo)，進(jìn)行精度評定，結(jié)果見表1。

錯檢率：實(shí)際未變化像元被錯分為變化像元的個數(shù)與實(shí)際未變化像元總數(shù)的比率。

漏檢率：實(shí)際變化像元被錯分為未變化像元的個數(shù)與實(shí)際變化像元總數(shù)的比率。

總正確率：兩類檢測正確的像元和與總像元的比率。

Kappa系數(shù)：其值越大，檢測精度越高。

圖8 廣州數(shù)據(jù)集變化檢測結(jié)果

圖9 貴陽數(shù)據(jù)集變化檢測結(jié)果

數(shù)據(jù)集方法錯檢率/%漏檢率/%總正確率/%Kappa系數(shù)廣州融合ASM+GGKI2.3720.1096.730.6941融合ENT+GGKI1.8324.9197.000.7011融合IDM+GGKI1.7035.4897.300.7068貴陽融合ASM+GGKI4.5720.7292.940.7339融合ENT+GGKI3.3318.4394.340.7829融合IDM+GGKI2.5718.6894.950.8025

從表1中可以看出，針對高分辨率SAR影像：①融合紋理信息與廣義高斯模型的變化檢測方法其Kappa系數(shù)均達(dá)到70%左右，能夠滿足變化檢測需求。②不同紋理信息的檢測性能不同，在具體應(yīng)用時可以針對特定需求選擇不同的紋理特征參與融合。③融合IDM紋理信息的檢測方法具有最低的錯檢率，并且其總正確率和Kappa系數(shù)最高，檢測效果最好。

針對水域界限，本文算法能夠較好的提取其變化信息。水體多為面狀和條帶狀，在SAR影像上呈暗色調(diào)，灰度值較小，與周圍其他地物具有很高的對比度，容易區(qū)分與識別。本文在計算GLCM時窗口大小取11×11，由于水體區(qū)域的灰度值變化幅度較小，所以計算出的紋理特征更加趨于同質(zhì)，更加便于區(qū)分水體與非水體，通過紋理差值就能發(fā)現(xiàn)變化的區(qū)域。再通過小波變換融合紋理差異圖和灰度差異圖，進(jìn)一步增強(qiáng)變化信息，抑制背景中的噪點(diǎn)，提高變化檢測性能。

3 結(jié)束語

本文充分利用高分辨率SAR影像豐富的紋理信息，以SWT變換為工具，將紋理差異圖與對數(shù)比值圖融合。融合后的差異圖有效地利用兩類差異圖的優(yōu)勢信息，充分反映變化信息的位置與范圍，性能更加優(yōu)異。同時，利用基于廣義高斯模型的KI閾值分割方法，提取變化信息，得到變化檢測二值圖。本文選用兩組3 m分辨率的TerraSAR-X影像，對比分析融合不同紋理信息的檢測性能，實(shí)驗(yàn)結(jié)果證明廣義高斯模型能夠較好地擬合出融合紋理信息差異圖的概率分布。同時紋理融合與廣義高斯模型相結(jié)合的變化檢測方法具有較可靠的檢測精度，應(yīng)用于高分辨率SAR影像變化檢測可行且有效。

參考文獻(xiàn)：

[1] ZHAO J,YANG J,LU Z,et al.A Novel Method of Change Detection in Bi-Temporal PolSAR Data Using a Joint-Classification Classifier Based on a Similarity Measure[J].Remote Sensing,2017,9(8): 846.

[2] 劉瑩,陶超,閆培,等.圖割能量驅(qū)動的高分辨率遙感影像震害損毀建筑物檢測[J].測繪學(xué)報,2017,46(7):910-917.

[3] 楊育麗,馬明國,周兆葉.基于地形梯度的山地流域土地利用格局及時空演變研究-以云南小江流域?yàn)槔齕J].測繪工程,2017,26 (12): 8-16.

[4] SCHLAFFER S,MATGEN P,HOLLAUS M,et al.Flood Detection from Multi-temporal SAR Data Using Harmonic Analysis and Change Detection[J].International Journal of Applied Earth Observation and Geoinformation,2015(38): 15-24.

[5] LEFEBVRE A,CORPETTI T,HUBERT-MOY L.Object-oriented Approach and Texture Analysis for Change Detection in Very High Resolution Images[C]//Geoscience and Remote Sensing Symposium,2008,IV-663-IV-666.

[6] TISON C,NICOLAS J M,TUPIN F,et al.A New Statistical Model for Markovian Classification of Urban Areas in High-resolution SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(10): 2046-2057.

[7] WANG G G,HUANG X T,ZHOU Z M,et al.A New SAR Image Change Detection Algorithm Based on Texture Feature [C]//Synthetic Aperture Radar (APSAR),2011 3rd International Asia-Pacific Conference on.IEEE,2011: 1-4.

[8] SUN W,SHI L,YANG J,et al.BuildingCollapse Assessment in Urban Areas Using Texture Information from Postevent SAR Data[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2016,9(8): 3792-3808.

[9] SHI L,SUN W,YANG J,et al.Building Collapse Assessment by the Use of Postearthquake Chinese VHR Airborne SAR[J].IEEE Geoscience and Remote Sensing Letters,2015,12(10): 2021-2025.

[10] 陳志鵬,鄧鵬,種勁松,等.紋理特征在SAR圖像變化檢測中的應(yīng)用[J].遙感技術(shù)與應(yīng)用,2002,17(3): 162-166.

[11] 李坤,楊然,王雷光,等.基于敏感特征向量的SAR圖像災(zāi)害變化檢測技術(shù)[J].計算機(jī)工程與應(yīng)用,2012,48(3):24-28.

[12] WANG G,HUANG X,ZHOU Z,et al.A New SAR Image Change Detection Algorithm Based on Texture Feature[C]//Synthetic Aperture Radar (APSAR),2011 3rd International Asia-Pacific Conference on.IEEE,2011: 1-4.

[13] 趙凌君,秦玉亮,高貴,等.利用 GLCM 紋理分析的高分辨率 SAR 圖像建筑區(qū)檢測[J].遙感學(xué)報,2009 (3): 475-490.

[14] BAZI Y,BRUZZONE L,MELGANI F.An Unsupervised Approach Based on The Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(4): 874-887.

[15] 袁增云,胡召玲.基于廣義高斯模型的多時相 SAR 圖像變化檢測[J].江蘇師范大學(xué)學(xué)報(自然科學(xué)版),2013,31(1): 74-78.

[16] 胡召玲.廣義高斯模型及 KI 雙閾值法的 SAR 圖像非監(jiān)督變化檢測[J].測繪學(xué)報,2013 (1): 116-122.

[17] ZHAO J Q,YANG J,LI P X,et al.An Unsupervised Change Detection Based on Test Statistic and KI from Multi-temporal and Full Polarimetric SAR Images[J].ISPRS-International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2016: 611-615.

[18] ZHAO J Q,YANG J,LU Z,et al.An Unsupervised Method of Change Detection in Multi-Temporal PolSAR Data Using a Test Statistic and An Improved K&I Algorithm[J].Applied Sciences,2017,7(12): 1297.

[19] HARALICK R M,SHANMUGAM K.Textural Features for Image Classification[J].IEEE Transactions on systems,man,and cybernetics,1973 (6): 610-621.

[20] GONG M G,ZHOU Z Q,MA J J.Change Detection in Synthetic Aperture Radar Images Based on Image Fusion and Fuzzy Clustering[J].IEEE Transactions on Image Processing,2012,21(4):2141-2151.

[21] 慕彩紅,霍利利,劉逸,等.基于小波融合和PCA-核模糊聚類的遙感圖像變化檢測[J].電子學(xué)報,2015,43(7):1375-1381.