專利名稱：利用星載sar多時相圖像識別地表變化的方法
技術(shù)領(lǐng)域：
本發(fā)明屬于空間遙感與空對地觀測信息處理技術(shù)領(lǐng)域，具體涉及一種利用星載SAR多時相圖像識別地表變化的方法。
背景技術(shù)：
城市與城市群體的發(fā)展是本世紀(jì)人類文明發(fā)展的重大內(nèi)容。城市建筑、綠地、道路的發(fā)展，水體、植物生態(tài)、氣候的變化，人口的遷移等等，所引發(fā)產(chǎn)生的城市地域與內(nèi)涵的變化都需要科學(xué)的評估。多時相(不同時間)的空間遙感觀測為監(jiān)測與評估這一發(fā)展提供了便捷可行的技術(shù)手段，它的一個直接應(yīng)用就是地表變化信息的識別與分類(Rignot andvan Zyl 1993，Singh 1989，Bruzzone and Prieto 2000，Merril and Jiajun 1998，Kasetkasem andVarshney 2002)。全天候高分辨率SAR成像近二十年的發(fā)展與空間遙感計劃(如SIR-CSAR，ERS-1，2，RADARSAT-1，2)的執(zhí)行能提供多年觀測的數(shù)據(jù)圖像，為地表變化信息的獲取與評估提供了條件。
比較不同時相的海量的遙感圖像數(shù)據(jù)，從直觀的圖像灰度值變化作人工定性的評估是十分粗糙的，難以準(zhǔn)確地確定區(qū)域的變化，判斷的閾值不可靠，極為容易失去重要的信息。因此，如何自動獲取與分析地表變化信息，一直是遙感研究的重要課題。
地表變化信息的檢測與分析有兩類方法有監(jiān)督的方法和無監(jiān)督的方法。前者需要多時相的大量實(shí)況數(shù)據(jù)作為分類過程的訓(xùn)練樣本，后者則不需要附加信息，直接對多時相的遙感圖像進(jìn)行分析。本發(fā)明采用無監(jiān)督方法。
Rignot和van Zyl(1993)曾用二幅不同時間的SAR圖像之間的差值作變化分析。Singh(1989)用紅外多通道數(shù)據(jù)進(jìn)行變化矢量分析，由各對應(yīng)像素計算得到的變化矢量幅值的統(tǒng)計來檢測變化。這些方法相對簡單易用，但不能對差值圖像進(jìn)行自動或非啟發(fā)式的分析，對差值圖像設(shè)定的閾值通常是基于經(jīng)驗(yàn)或人為的，而閾值選取適當(dāng)與否在很大程度上決定了檢測結(jié)果的可靠與準(zhǔn)確。
Bruzzone和Prieto(2000)提出了基于Bayes理論和馬爾科夫隨機(jī)場(MRF，MarkovRandom Fields)模型的自動檢測區(qū)域變化的方法，通過期望極大化(EM，ExpectationMaximum)算法來自動獲取一個閾值。他們將差值圖像分成兩類變化區(qū)和無變化區(qū)。但是，他們并沒有用真實(shí)SAR圖像數(shù)據(jù)來說明這一算法的可靠性與可行性。同時，即使是變化區(qū)，也可以有完全不同的變化。比如樓群變成草地，后向散射可能變?nèi)酰徊莸厣嫌袠浠蜃兂呻p向反射增強(qiáng)的樓群，后向散射可能增強(qiáng)。后向散射同一db的增強(qiáng)或減弱完全可以有不同的地表變化的信息內(nèi)涵。隨著多通道、多極化、多類遙感器多時相SAR圖像的獲取，自動檢測、識別、分類與評估地表變化更多的信息一定會成為重要的研究方向。

發(fā)明內(nèi)容
本發(fā)明的目的在于提出一種利用星載SAR多時相圖像正確、有效識別地表變化的方法。
本發(fā)明提出的利用星載SAR多時相圖像識別地表變化的方法，是一種基于星載SAR多時相微波遙感圖像的EM-MRF方法。其基本步驟為(1)選取不同時間的兩幅星載SAR微波遙感圖像，計算其差值圖像，該差值圖像是兩幅圖像的相減，或者是兩張圖像的相除；(2)采用EM算法將差值圖像分成三類增強(qiáng)區(qū)、減弱區(qū)和不變區(qū)，并獲得區(qū)分的兩個閾值；(3)根據(jù)圖像上相鄰像素的空間結(jié)構(gòu)相關(guān)性，將EM算法的結(jié)果作為初值，利用MRF-ICM算法，進(jìn)行迭代計算，使差值圖像的Gibbs能量達(dá)到局部最小值，從而得到關(guān)于城市變化的散射增強(qiáng)、減弱、不變的分類結(jié)果。
本發(fā)明用1996和2002年上海市的ERS-2SAR圖像作了實(shí)例研究，用實(shí)地情況佐證。本發(fā)明同樣可應(yīng)用于我國長江三角洲、珠江流域等城市鏈發(fā)展的自動識別變化信息評估。在多通道、多極化、多時相SAR數(shù)據(jù)獲取后，可很容易地推廣到更多分類區(qū)域的分類檢測。


圖1.差值圖像直方圖。其中，(a)為差值圖像直方圖上閾值劃分示意，(b)為實(shí)際差值圖像歸一化直方圖上三類分布。
圖2.MRF的二階鄰域系統(tǒng)和組系。其中，(a)二階鄰域系統(tǒng)，(b)單像素組系L1及其權(quán)值α，(c)雙像素組系L2及其權(quán)值β。
圖3.1996年6月4日ERS-2SAR圖像。
圖4.2002年4月9日ERS-2SAR圖像。
圖5.二者的差值圖像。
圖6.不同迭代次數(shù)下各參數(shù)值(γ＝0.5)。
圖7.雙閾值EM算法地表變化分類。
圖8.MRF-ICM算法自動檢測地表變化分類。
圖9.上海市變化圖。
圖10為本發(fā)明流程框圖。
具體實(shí)施例方式
1、EM算法估計變化分類的先驗(yàn)、條件概率考慮不同時間(t1，t2)的兩張已經(jīng)對準(zhǔn)的大小為I×J圖像X1和X2。其差值圖像X0＝(X2-X1)(也可用X0＝(X2/X1))，設(shè)X是一個隨機(jī)變量，代表了差值圖像中I×J個像素的值，XD＝{X(i，j)，l≤i≤I，l≤j≤J}。I、J分別為圖像的高度和寬度。
檢測與分析差值圖像是要區(qū)分各像素屬于不變化類ωn或變化類ωc。差值圖像XD的像素值的概率分布函數(shù)p(X)包含了這兩類的條件概率分布，即有p(X)＝p(X|ωn)P(ωn)+p(X|ωc)P(ωc)(1)這里P(ωn)，P(ωc)是“未變化”和“變化”兩類的先驗(yàn)概率分布函數(shù)，p(X|ωk)(k＝n，c)是條件概率分布函數(shù)。以Bruzzon和Prieto(2000)的工作作為用期望極大化(EM)算法的基礎(chǔ)，對p(X|ωn)、p(X|ωc)、p(ωn)和P(ωc)進(jìn)行無監(jiān)督估計。
EM算法是估計不完整性數(shù)據(jù)問題的最大似然估計法的一種普遍的方法(Takajo andTakahashi 1998，Ghiglia and Romero 1994，Pritt and Shipman 1994，Pritt 1996)。它包括了一個估計步驟和一個極大化步驟，兩個步驟都需要迭代到收斂。估計步驟用受到觀測值控制的參數(shù)估計作為未知變量來進(jìn)行計算。極大化步驟提供了對參數(shù)的新的估計。假設(shè)p(X|ωn)和p(X|ωc)服從高斯(Gauss)分布，ωn對應(yīng)的概率密度函數(shù)由均值μn和方差σn2來描述，同樣，與ωc相關(guān)的概率密度函數(shù)由均值μc和方差σc2描述?？梢宰C明(Render andWalker 1984)，估計ωn的各參數(shù)與先驗(yàn)概率分布有公式如下Pt+1(ωn)=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))IJ···(2)]]>μnt+1=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))X(i,j)ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))···(3)]]>
(σn2)t+1=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j)[X(i,j)-μnt]2ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))···(4)]]>上標(biāo)t和t+1代表當(dāng)前t和下一次t+1的迭代。類似地，這些方程也用來估計ωc的先驗(yàn)概率和條件概率密度函數(shù)的均值μc和方差σc2。即將式(2)、(3)、(4)中的n改為C即可，依次記分(2)’、(3)’、(4)’。
Pt+1(ωc)ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))IJ···(2),]]>μct+1=ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))X(i,j)ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))···(3),]]>(σc2)t+1=ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j)[X(i,j)-μct]2ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))···(4),]]>參數(shù)估計先從其初始值開始通過迭代上述公式到收斂來獲得?？梢宰C明(Takajo andTakahashi 1998，Ghiglia and Romero 1994)，每迭代一次，估計的參數(shù)獲得一個對數(shù)函數(shù)形式L(θ)＝ln p(XD|θ)的增加量，其中θ={P(ωn),P(ωc),μn,μc,σn2,σc2},]]>收斂時將獲得對數(shù)函數(shù)形式的本地最大值。
參數(shù)估計的初始值可以通過分析差值圖像直方圖h(X)來獲得。在直方圖上確定兩個閾值Tn和Tc來獲得一個比較接近ωn的像素子集Sn和另一個比較接近ωc的像素子集Sc。我們將Tn和Tc寫成T＝MD(1-γ)和T＝MD(1+γ)，其中MD是h(X)的中值。(即MD＝[max{XD}+min{XD}]/2)，γ∈(0，1)是一個用來定義那些不容易識別的像素在MD周圍變化范圍的初始化參數(shù)。那么集合Sn＝{X(i，j)|X(i，j)＜Tn}和Sc＝{X(i，j)|X(i，j)＞Tc}就可以用來計算兩類ωn和ωc各自的初始統(tǒng)計參數(shù)。
在各像素值相互獨(dú)立的假設(shè)下，由Bayes最小誤差原理，要使差值圖像XD中的像素X(i，j)屬于類別ωk則須使得其后驗(yàn)條件概率最大，即寫為ωk=argmaxωi∈{ωn,ωc}{P(ωi|X(i,j))}=argmax{P(ωi)p(X(i,j)|ωi)}ωi∈{ωn,ωc}···(5)]]>
這等價于在差值圖像上對兩類ωn和ωc求解最大似然邊界To。最佳閾值To由如下方程給出P(ωc)P(ωn)=p(X|ωn)p(X|ωc)···(6)]]>在Gauss概率密度分布的假設(shè)下，由式(6)可得(σn2-σc2)To2+2(μnσc2-μcσn2)To+μc2σn2-μn2σc2-2σn2σc2ln[σcP(ωn)σnP(ωc)]=0···(7)]]>求得ToTo=-b±b2-4ac2a,]]>a=σn2-σc2,b=2(μnσc2-μcσn2),]]>C=μc2σn2-μn2σc2-2σn2σc2ln[σcP(ωn)σnP(ωc)]···(7),]]>2、雙閾值的EM算法SAR圖像各像素值表示雷達(dá)回波(后向散射)的強(qiáng)弱，兩幅不同時間的SAR圖像的差值可表明各像素點(diǎn)上后向散射增強(qiáng)、不變和減弱三種情況。增強(qiáng)與減弱表現(xiàn)了地表變化的不同信息。我們將差值圖像各點(diǎn)分為三類增強(qiáng)ωc1、不變ωn和減弱ωc2。在XD中即分別為正的點(diǎn)集、0周圍的點(diǎn)集和負(fù)的點(diǎn)集。
于是XD的概率密度函數(shù)p(X)為p(X)＝p(X|ωc1)P(ωc1)+p(X|ωn)P(ωn)+p(X|ωc2)P(ωc2)(8)采用EM算法，首先對增強(qiáng)類ωc1和非增強(qiáng)類ωn1(ωn和ωc2)得到閾值To1，再對減弱類ωc2和非減弱類ωn2(ωn和ωc1)得到閾值To2，即ωn1＝ωn∪ωc2，ωn2＝ωn∪ωc1。將這兩個結(jié)果合并得到三種分類。
如圖1(a)，直方圖h(X)中增大區(qū)(x≥0)內(nèi)用兩個閾值Tn1和Tc1來獲得一個比較接近ωn1的像素子集Sn1和另一個比較接近ωc1像素子集Sc1。將Tn1和Tc1寫成Tn1＝MD1(1-γ)和Tn1＝MD1(1+γ)，其中MD1是h(X)中增強(qiáng)區(qū)的中值(即MD1＝[max{XD}/2])，本文中γ取0.5。同樣，在h(X)中減弱區(qū)(X＜0)內(nèi)找到比較接近ωn2和ωc2的像素子集Sn2和Sc2，其中MD1＝[min{XD}/2]。
圖1(b)是以后要用到的圖5ERS-2SAR對上海市觀測圖像實(shí)際數(shù)據(jù)中截取的延中綠地后向散射差值圖像的灰度直方圖。為便于比較，在圖1(b)中用曲線大致標(biāo)出了三類的概率密度分布。直方圖進(jìn)行歸一化，使得直方圖可以與概率密度分布相一致。
在增強(qiáng)區(qū)內(nèi)迭代計算時，不管減弱區(qū)的像素；在減小區(qū)內(nèi)迭代時，不管增強(qiáng)區(qū)的像素。那么得到集合Sn1＝{X(i，j)|0＜X(i，j)＜Tn1}，Sc1＝{X(i，j)|X(i，j)＞Tc1}和Sn2＝{X(i，j)|Tn2＜X(i，j)＜0}，Sc2＝{X(i，j)|x(i，j)＜Tc2}，用來計算兩類ωn1ωc1和ωn2，ωc2各自相關(guān)的初始統(tǒng)計參數(shù)。
確定了初始參數(shù)后，用式(2～4)進(jìn)行上述四類的迭代，收斂后得到四類的先驗(yàn)概率密度函數(shù)及其統(tǒng)計參數(shù)[P(ωn1)，μn1，σn12]，[P(ωc1)，μc1，σc12]，[P(ωn2)，μn2，σn22]和[P(ωc2)，μc2，σc22]。再先后求解式(7)，分別得到閾值To1和閾值To2。
3、Markov隨機(jī)場(MRF)對空間結(jié)構(gòu)相關(guān)的變化分類Bruzzone和Prieto(2000)使用了兩種方法來檢測多時相遙感圖像上的區(qū)域變化，第一種假設(shè)差值圖像各像素之間是統(tǒng)計獨(dú)立的，第二種則認(rèn)為各像素是空間結(jié)構(gòu)相關(guān)的。事實(shí)上，一個像素屬于ωk類，那么它相鄰的像素也很有可能屬于ωk類。使用第二種方法能改善分類的性能，得到更可靠的結(jié)果。
用集合C＝{Cl，1≤l≤L}表示差值圖像上所有分類結(jié)果的集合，Cl＝{Cl(i，j)，1≤i≤I，1≤j≤J}，Cl(i，j)∈{ωc1，ωc2，ωn}其中L＝3IJ，I，J分別為圖像的高度和寬度。
根據(jù)Bayes最小誤差理論，分類結(jié)果應(yīng)使后驗(yàn)概率為最大，即Ck=argmaxCl∈C{P(Cl|XD)}=argmaxCl∈C{P(Cl)p(XD|Cl)}···(9)]]>其中XD為差值圖像，P(C1)是類別標(biāo)注的先驗(yàn)概率，p(XD|Cl)是給定類別標(biāo)注條件下的差值圖像的像素值的聯(lián)合概率密度函數(shù)。式(9)的求解需要同時估計P(Cl)和p(XD|Cl)，計算上會過于復(fù)雜。采用Markov隨機(jī)場(MRF)來建立差值圖像上像素間的空間結(jié)構(gòu)相關(guān)的模型，問題就可以簡化。
為了用MRF來描述問題，需要引進(jìn)空間鄰域系統(tǒng)N?？紤]一個坐標(biāo)為(i，j)的像素，它的鄰域可以用下式表達(dá)N(i，j)＝{(i，j)+(v，ζ)，(v，ζ)∈N}(10)
在本發(fā)明中，N＝{(±1，0)，(0，±1)，(1，±1)，(-1，±1)}，即一個像素的鄰域由它周圍的8個像素構(gòu)成。
用Cl(i，j)表示像素點(diǎn)(i，j)的類別標(biāo)注，ClN表示像素點(diǎn)(i，j)的鄰域像素的類別標(biāo)注，Cl|(i，j)表示差值圖像XD中除了像素點(diǎn)(i，j)之外的其余像素的類別標(biāo)注。把Cl看成一個MRF，則Cl由如下兩條性質(zhì)P(Cl)＞0，_Cl(i，j)∈{ωc1，ωc2，ωn}(11)P(Cl(i，j)|Cl|(i，j))＝P(Cl(i，j)|ClN) (12)MRF的二階鄰域系統(tǒng)和組系如圖2所示。
組系是指鄰域系統(tǒng)中的包含中心像素點(diǎn)的鄰近像素的連通子圖，這里考慮單像素和雙像素的組系。如圖2(b，c)所示，單像素的組系有1種形狀，雙像素的組系有4種形狀，其數(shù)學(xué)表示分別為L1＝(i，j)，L2＝{{(i，j)，(g，h)}|(g，h)∈Ni(i，j)}。
根據(jù)Hammersley-Clifford定理(Besag 1974)，MRF可以等效為Gibbs分布，因此在給定其它像素的類別標(biāo)注的前提下，坐標(biāo)(i，j)處的像素的類別標(biāo)注的條件概率分布可表示為P(Cl(i，j)|{Cl(g，h)，(g，h)≠(i，j)})＝P(Cl(i，j))|{Cl(g，h)，(g，h)∈N(i，j)})=1Zexp[-U(Cl(i,j)|{Cl(g,h),(g,h)∈N(i,j)})]···(13)]]>其中U(·)是Gibbs能量函數(shù)，Z是歸一化因子，Z=Σ(g,h)∈{(i,j)∪N(i,j)}exp[-U(Cl(i,j)|{Cl(g,h)})]···(14)]]>U(Cl|XD)=Σ1≤i≤IΣ1≤j≤JU(Cl(i,j)|Cl|(i,j),XD)=Σ1≤i≤IΣ1≤j≤JU(Cl(i,j)|ClN,XD)···(15)]]>在差值圖像的直方圖為Gauss分布的情況下，U(Cl(i,j)|ClN,XD)=12ln|2πσCl(i,j)2|+(X(i,j)-μCl(i,j))22σCl(i,j)2]]>+v1(Cl(i,j))+Σ(g,h)∈N(i,j)v2(Cl(i,j),Cl(g,h))---(16)]]>其中v1(Cl(i，j))＝α，(i，j)∈L1。這里α為單像素組系的權(quán)值，一般為-0.7-0.5之間，實(shí)例計算中取α＝0。
v2(Cl(i,j),Cl(g,h))=-β,Cl(i,j)=Cl(g,h),{(i,j),(g,h)}∈L2β,Cl(i,j)≠Cl(g,h),{(i,j),(g,h)}∈L2···(17)]]>這里β是雙像素組系的權(quán)值，一般在0.7-1.5之間，實(shí)例計算中取β＝1.0。
根據(jù)式(9)，區(qū)域變化圖的生成就是給差值圖像上各像素點(diǎn)標(biāo)上類別，使得后驗(yàn)概率最大。在Markov模型下，就是要使式(15)的能量函數(shù)最小(Bruzzone and Prieto 2000)，重寫式(15)為U(Cl|XD)=Σ1≤i≤IΣ1≤j≤J[UI(X(i,j)/Cl(i,j))]]>+UC(Cl(i,j)/{Cl(g,h),(g,h)∈N(i,j)}])···(18)]]>其中UC(·)表示考慮了像素之間類別相關(guān)的能量函數(shù)，由式(16)的后兩項(xiàng)給出。UI(·)表示假定像素之間類別獨(dú)立的能量函數(shù)，由式(16)的前兩項(xiàng)給出，即UI(X(i,j)/Cl(i,j))=12ln|2πσCl(i,j)2|+12(X(i,j)-μCl(i,j))2[σCl(i,j)2]-1···(20)]]>其中，σCl(i,j)2∈{σcl2,σc22,σn2},]]>μCl(i,j)∈{μcl,μc2,μn},]]>由EM算法給出。
一般地，式(18)的求解需要用模擬退火等迭代算法，本文采用ICM算法(Besag1986)，此算法計算量小，可以迅速收斂到能量函數(shù)的局部極小值。ICM算法的數(shù)學(xué)表達(dá)式可寫成Cl=ICM(αβμCl(i,j),σCl(i,j)2,XD,Cl0)···(21)]]>其中Cl0為Cl的初始值。
求解(18)得到Cl的步驟如下(a)通過EM算法，得到Cl的初始值Cl0。
(b)計算出新的Cl(i，j)，使得式(18)中的能量函數(shù)最小。
(c)重復(fù)步驟(b)直到收斂。
ICM算法收斂的條件是迭代次數(shù)超過30次或少于 的像素點(diǎn)發(fā)生變化。ICM算法的具體描述可參見Besag(1986)。
4、城市變化識別實(shí)例圖3、4分別為1996年6月4日和2002年4月9日中國上海地區(qū)的ERS-2 SAR遙感圖像，圖像的分辨率是12.5m×12.5m。計算前圖像已經(jīng)校準(zhǔn)，并對初始數(shù)據(jù)進(jìn)行了3×3的均值濾波以消除噪聲起伏。我們從原圖中分別切下580像素×980像素的兩塊區(qū)域，中心是上海市的陸家嘴地區(qū)。圖5是二者的差值圖像(也有將二者相除作為差值圖像)。差值圖像中灰度級大的點(diǎn)，代表該點(diǎn)上的后向散射的差值大。從差值圖中可以看到上海市這5年中的變化，但是如何確定變化的分類則是很難自動確定的。
對差值圖像進(jìn)行EM迭代求解，得到各類在迭代過程中先驗(yàn)概率和統(tǒng)計參數(shù)在各步驟上的值。第一類迭代過程參見圖6(a～c)，第二類迭代過程類似。
在完成EM求解后，分別得到散射增強(qiáng)區(qū)和減弱區(qū)內(nèi)兩個類的先驗(yàn)概率和概率密度分布的均值和方差，根據(jù)公式(7)獲得兩個閾值To1，To2用于分類。
經(jīng)計算，To1＝0.946×10-4，To2＝-1.013×10-4(后向散射差值最大值為5.202×10-4，最小值為-3.381×10-4)。
圖7是采用EM算法所得的區(qū)域變化的檢測結(jié)果。圖中綠色為后向散射增強(qiáng)區(qū)，紅色代表后向散射減弱區(qū)，藍(lán)色代表后向散射變化不大的區(qū)。這一分類結(jié)果比差值圖像更加清晰。
為納入空間結(jié)構(gòu)相關(guān)信息，圖8利用基于MRF模型的ICM算法得到結(jié)構(gòu)相關(guān)區(qū)域變化的自動檢測。與圖7相比，同一類別的像素點(diǎn)更加集中，孤立的像素點(diǎn)更少，符合MRF模型的理論。比如，圖8中奇數(shù)(1，3，5)標(biāo)注的地區(qū)為平坦表面后向散射減弱的代表區(qū)域，其中1處為陸家嘴綠地，3處為浦東世紀(jì)公園，5處為延中綠地。偶數(shù)(2，4)標(biāo)注的地區(qū)為多類建筑物產(chǎn)生雙向反射和散射增強(qiáng)的代表區(qū)域，其中2處為浦東陸家嘴的建筑樓群，4處為上海新國際展覽中心。圖8所顯示的上海市1-5地區(qū)的變化與上海市從1996到2002年的建設(shè)情況是一致的，說明了本文雙閾值EM-MRF檢測算法是有效的。將圖8的結(jié)果覆蓋在上海市旅游圖上，如圖9，是一個有趣的比較。
利用多時相SAR微波遙感圖像提出的EM-MRF模型的算法，可自動檢測多年間城市地區(qū)的變化。這種方法可以對差值圖像進(jìn)行自動分析，不需要作經(jīng)驗(yàn)估計或人工干預(yù)。在假設(shè)像素之間統(tǒng)計獨(dú)立的前提下，按照后向散射變化的強(qiáng)弱，用EM算法將差值圖像分成三類增強(qiáng)區(qū)、減弱區(qū)和持不變區(qū)，獲得兩個閾值。再進(jìn)一步根據(jù)圖像上相鄰像素的空間結(jié)構(gòu)相關(guān)性，將EM算法的結(jié)果作為初始值，利用MRF-ICM算法進(jìn)行迭代計算，使差值圖像的Gibbs能量到達(dá)局部最小值，從而得到關(guān)于城市變化的散射增強(qiáng)、減弱、不變的分類結(jié)果。用1996和2002年的上海地區(qū)的ERS-2的多時相遙感圖像進(jìn)行了算法驗(yàn)證，對比實(shí)地發(fā)展變化，證實(shí)了這種方法的有效性。本方法在獲得多極化、多通道、多類遙感器多時相遙感圖像后，可進(jìn)一步推廣到多閾值的EM-MRF算法，得到更多的信息融合的地表變化信息自動分類。
參考文獻(xiàn)Besag J.(1974)，Spatial interaction and the statistical analysis of lattice systems(withdiscussion).Journal of Royal Statistics Society，Series B，36(2)192-326.
Besag J.(1986)，On the statistic analysis of dirty pictures.Journal of Royal Statistics Society，Series B，48259-302.
Bruzzone L.and D.F.Prieto(2000)，Automatic analysis of the difference image forunsupervised change detection.IEEE Transactions on Geoscience and Remote Sensing，38(3)1171-1182.
Dempster A.P.，N.M.Laird，and D.B.Rubin(1977)，Maximum linkelihood from incompletedata via the EM algorithm.Journal of Royal Statistics Society，39(1)1-38.
Fukunaga K.(1990)，Introduction to Statistical Pattern Recognition，2nded.London，U.K.Academic.
Jin Y.Q.，F(xiàn).Chen and L.Luo(2004)，Automatic analysis of change detection of multi-temporalERS-2SAR images by using two-thresholds EM and MRF algorithms，The Imaging ScienceJournal，52234-241.
Kasetkasem T.and P.Varshney(2002)，An image change detection algorithm based on Markovrandom field model.IEEE Transactions on Geoscience and Remote Sensing，40(8)1815-1823.
Merril K.R.and L.Jiajun(1998)，A comparison of four algorithms for change detection in anurban environment.Remote Sensing of Environment，6395-100.
Moon T.K.(1996)，The expectation-maximization algorithm.Signal Processing Magazine，13(6)47-60.
Render A.P.and H.F.Walker(1984)，Mixture densities，maximum likelihood and the EMalgorithm.SIAM Review，26(2)195-239.
Rignot E.and J.van Zyl(1993)，Change detection techniques for ERS-1SAR data.IEEETransactions on Geoscience and Remote Sensing，31(4)896-906.
Shahshahani B.M.and D.A.Landgrebe(1994)，The effect of unlabeled samples in reducing thesmall size problem and mitigating the Hughes phenomenon.IEEE Transactions onGeoscience and Remote Sensing，321087-1095.
Singh A.(1989)，Digtial change detection techniques using remotely-sensed data.InternationalJournal of Remote Sensing，10(6)989-1003.
權(quán)利要求
1.一種利用星載SAR多時相圖像識別地表變化的方法，其特征在于基本步驟如下(1)選取不同時間的兩幅星載SAR微波遙感圖像，計算其差值圖像，該差值圖像是兩幅圖像的相減，或者是兩張圖像的相除；(2)采用EM算法將差值圖像分成三類增強(qiáng)區(qū)、減弱區(qū)和不變區(qū)，并獲得區(qū)分的兩個閾值；(3)將EM算法的結(jié)果作為初值，利用MRF-ICM算法，進(jìn)行迭代計算，使差值圖像的Gibbs能量達(dá)到局部最小值，從而得到關(guān)于城市變化的散射增強(qiáng)、減弱、不變的分類結(jié)果。
2.根據(jù)權(quán)利要求1所述的利用星載SAR多相圖像識別地表變化的方法，其特征在于所述兩個閾值由下式計算獲得(σn2-σc2)To2+2(μnσc2-μcσn2)To+μc2σn2-μn2σc2-2σn2σc2ln[σcP(ωn)σnP(ωc)]=0---(7)]]>而σn、μn、p(wn)由Pt+1(ωn)=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))IJ---(2)]]>μnt+1=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))X(i,j)ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))---(3)]]>(σn2)t+1=ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))[X(i.j)-μnt]2ΣX(i,j)∈XDPt(ωn)pt(X(i,j)|ωn)pt(X(i,j))---(4)]]>迭代求得，σc、μc、p(wc)由Pt+1(ωc)=ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))IJ---(2)′]]>μct+1=ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))X(i,j)ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))---(3)′]]>(σc2)t+1=ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))[X(i.j)-μct]2ΣX(i,j)∈XDPt(ωc)pt(X(i,j)|ωc)pt(X(i,j))---(4)′]]>迭代求得，其中ωc為差值圖像XD的不變化類像素，ωc為差值圖像XD的變化類像素，p(X)為差值圖像XD的像素值的X的概率分布函數(shù)p(X)＝p(X|ωn)P(ωn)+p(X|ωc)P(ωc)(1)而p(wn)、p(wc)為wn、wc的先驗(yàn)概率分布函數(shù)，p(X|ωk)(k＝n，c)是條件概率分布函數(shù)。
全文摘要
本發(fā)明屬空間遙感和空對地觀測信息處理技術(shù)領(lǐng)域，具體為一種利用星載SAR多時相圖像識別地表變化的方法。本發(fā)明用不同時間SAR的圖像差值，通過求解雙閾值期望極大化和馬爾科夫隨機(jī)場EM－MRF算法，識別不同時間的地表發(fā)生散射增強(qiáng)、減弱與不變?nèi)愖兓?，?996和2002年上海市的歐洲遙感衛(wèi)星ERS－2 SAR圖像作了實(shí)例研究，證明本發(fā)明方法的有效性。
文檔編號G06K9/00GK1760888SQ200510030999
公開日2006年4月19日 申請日期2005年11月3日 優(yōu)先權(quán)日2005年11月3日
發(fā)明者金亞秋 申請人:復(fù)旦大學(xué)