專利名稱：：基于非參數密度估計的遙感圖像變化檢測方法
技術領域：
：本發(fā)明屬于數字圖像處理
技術領域：
，涉及多時相遙感圖像的變化檢測，具體的說是一種基于非參數密度估計的遙感圖像變化檢測。
背景技術：
：變化檢測技術是指通過分析在同一地區(qū)但不同時間獲得的兩幅圖像來辨識其變化信息。隨著遙感圖像獲取技術和手段的日益先進以及遙感圖像數據的海量積累，變化檢測技術在環(huán)境監(jiān)測、土地利用/覆蓋、森林/植被變化分析、災情監(jiān)測、農業(yè)調查、城市變化分析、軍事偵察和打擊效果評估等方面的應用越來越廣泛。在已發(fā)表的文獻中，基于非監(jiān)督的變化檢測技術主要基于以下3個步驟1)圖像的預處理，包括輻射校正、幾何配準、圖像去噪等；2)差異圖像的構建，具體指的是將兩幅圖像進行逐個像素的比較；3)變化區(qū)域的提取，主要包括閾值法和分類法，其中基于MRF(MarkovRandomFields)模型的分類方法，由于顧及了上下文關系，有較強的抗噪性，得到了一些學者的關注。Bruzzone禾口Prieto(2000)在文章"Automaticanalysisofthedifferenceimageforunsupervisedchangedetection"中提出了基于Bayes理論禾口MRF模型的非監(jiān)督變化檢測方法，假設差異圖像中與變化類和非變化類相關的統(tǒng)計項符合高斯混合模型(G匪，GaussianMixtureModels)，并采用期望最大化(EM，ExpectationMaximum)算法來估計其模型參數，最后分別采用貝葉斯最小錯誤率閾值和MRF對差異圖像進行分類。作為該方法的進一步改進，2002年Bruzzone禾口Prieto又在文章"Anadaptivesemiparametricandcontext—basedapproachto皿supervisedchangedetectioninmultitemporalremote-sensingimages"采用了簡化的Parzen估計和EM算法來估計差異圖像中與變化類和非變化類像素灰度級相關的統(tǒng)計項，但由于統(tǒng)計項的非監(jiān)督估計和MRF空間正則化的過程是分離的，所以變化檢測處理效率低。江利明，廖明生(2006)等在文章"顧及空間鄰域關系的多時相SAR圖像變化檢測"提出基于EM-MPM模型的變化檢測方法，并和雙閾值EM算法進行了比較，有效地提高了變化區(qū)域提取的可靠性和準確性。孫強(2007)在其博士論文"基于統(tǒng)計模型的SAR圖像處理與解譯"中提出了一種基于廣義高斯混合模型的SAR圖像變化檢測方法。在GGM(GeneralGaussMixture，GGM)的先驗下，通過基于Gibbs采樣估計方法的模型統(tǒng)計推斷，對兩幅相關SAR圖像的對數比圖像進行最大似然分類，并在此基礎上基于MRF進行自適應的空間約束，完成檢測結果的更新。宋妍，袁修孝(2009)等在文章"基于混合高斯密度模型和空間上下文信息的遙感圖像變化檢測方法及擴展"中提出了一種運用遺傳K均值算法與EM算法聯(lián)合解算高斯混合密度模型參數的方法，該方法可以自動地解算出模型的統(tǒng)計參數；然后，比較概率松弛迭代法以及MRF模型法的圖像變化檢測效果；最后，對傳統(tǒng)的基于模擬退火法的MRF方法進行改進，提出了一種變權MRF方法，檢測結果能更好地保持圖像的結構性，并有效地去除了孤立噪聲。以上方法假設了差異圖像中與變化類和非變化類相關的統(tǒng)計項符合具體的模型如高斯混合模型、廣義高斯混合模型等，需要進行復雜的參數估計過程，且參數估計的精確程度會影響變化檢測的結果，而實際中差異圖像的統(tǒng)計項不一定符合這些具體的模型，使得這些方法對差異圖像中與變化類和非變化類相關的統(tǒng)計項的估計存在偏差，進而影響變化檢測精度。
發(fā)明內容本發(fā)明的目的在于克服上述已有的遙感圖像變化檢測技術的不足，提出一種基于非參數密度估計的遙感圖像變化檢測，以減小差異圖像中與變化類和非變化類相關統(tǒng)計項的估計偏差，提高變化檢測精度。實現本發(fā)明目的的技術方案是采用非參數密度估計方法估計差異圖像與變化類和非變化類相關的統(tǒng)計項，并結合宋妍，袁修孝(2009)等人提出的變權馬爾科夫隨機場(MarkovRandomFields)進行自適應的空間約束，對遙感圖像的變化進行檢測，其實現步驟包括如下(1)輸入兩幅不同時相的遙感圖像，并對每幅圖像的每個通道分別進行窗口大小為3X3像素的中值濾波，得到兩時相的去噪后圖像；(2)將去噪后的兩幅圖像應用變化矢量分析，得到一幅差異圖像，并根據該差異圖像計算變權馬爾科夫隨機場的權值因子W;(3)應用K-均值聚類算法將差異圖像聚成變化類和非變化類，得到初始分類結果；(4)利用初始分類結果，采用非參數密度估計方法估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數得到變化類和非變化類的似然能量；(5)對初始分類結果利用馬爾科夫隨機場計算變化類和非變化類的先驗能量；(6)利用權值因子W、變化類和非變化類的似然能量及變化類和非變化類的先驗能量計算變化類的總能量和非變化類的總能量，將總能量較小的那一類作為當前類別，得到類別更新后的結果；(7)對類別更新后的結果，采用非參數密度估計方法重新估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數得到變化類和非變化類的似然能量，并利用馬爾科夫隨機場重新計算變化類和非變化類的先驗能量；(8)重復步驟(6)及步驟(7)直至迭代終止，并存儲每次類別更新后的結果，得到每個像素點的類別更新集合，該迭代終止條件為迭代次數不超過50次及兩次迭代之間相異的像元數目比例小于給定的閾值；(9)利用每個像素點的類別更新集合估計變化類的后驗概率和非變化類的后驗概率，將后驗概率較大的那一類作為該像素點的最終變化檢測結果。本發(fā)明與現有技術相比具有如下優(yōu)點(1)本發(fā)明由于采用非參數密度估計方法估計差異圖像的類條件概率密度，克服了現有技術采用高斯混合模型和廣義高斯混合模型假設的缺陷，不需要事先對遙感影像的5類條件概率密度做出假設，能夠得到精確的估計結果。(2)本發(fā)明由于結合了變權馬爾科夫隨機場進行自適應的空間約束來迭代更新變化檢測結果，使得檢測結果能更好地保持圖像的結構信息，并有效地去除孤立噪聲。(3)本發(fā)明由于將類別統(tǒng)計項的估計和自適應的空間約束融為一體，提高了變化檢測處理效率。圖1是本發(fā)明的實現流程示意圖；圖2是本發(fā)明第一組實驗的變化檢測結果圖圖3是本發(fā)明第二組實驗的變化檢測結果圖圖4是本發(fā)明第三組實驗的變化檢測結果圖具體實施例方式參照圖l，本發(fā)明的實施如下步驟l，輸入兩幅不同時相的遙感圖像，并對每幅圖像的每個通道分別進行窗口大小為3X3像素的中值濾波，得到兩時相的去噪后圖像&和X2。步驟2，將去噪后的兩幅圖像&和X2應用變化矢量分析，得到一幅差異圖像Xd，并根據該差異圖像計算變權馬爾科夫隨機場的權值因子W，具體步驟如下(2a)利用變化矢量分析法計算差異圖像Xd，即3(1)其中，Xn、X。及X13為圖像&的三個通道圖像;X21、X22及X23為圖像X2的三個通道(2b)計算權值因子W:首先，計算每個像素點的特征值，即W,力=J]2Ix0l，")—;(2)=1rt=l其中p為像素點局部窗口的大小，x(m，n)為局部窗口內每個像素點的灰度值，u(i，j)為局部窗口像素的均值；然后，利用整幅圖像中像素特征值t(i，j)的最大值及最小值將t(i，j)映射到區(qū)間上，得到每個像素點的權值因子W(i，j)，Vmin=0.5，Vmax=8。步驟3，應用K-均值聚類算法將差異圖像聚成兩類，將均值較大的一類作為變化類，均值較小的一類作為非變化類，得到初始分類結果。步驟4，利用初始分類結果，采用非參數密度估計方法估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數，得到變化類和非變化類的似然能量，具體步驟如下(4a)采用非參數密度估計方法估計差異圖像中變化類的類條件概率密度A;IwJ和非變化類的類條件概率密度》(J^I化)，即6<formula>formulaseeoriginaldocumentpage7</formula>其中，Sn和S。分別表示非變化類和變化類的像素集合，Nn和N。分別表示非變化類和變化類的像素數目，K(□)為高斯核函數，Hn和H。分別表示非變化類和變化類的自適應窗寬平滑參數，與像素數目及像素點的頻數f(Xij)有關，通過如下公式計算Hn=H0(a/Nn_f(X^P);(5)Hc=h0(a/Nc_f(XjP);(6)其中，H。，a禾Pp均為經驗常數，H。=1，a=40000，P=10;(4b)對非變化類和變化類的類條件概率密度取負自然對數，得到的非變化類的似然能量LEu(i，j)及變化類的似然能量LEc(i，j)，艮卩丄五"(;,/)=—ln(》(A^lw"》；(7)丄Ec(/,_;)=-ln(》(^^Ic))。(8)步驟5，對初始分類結果利用馬爾科夫隨機場計算變化類和非變化類的先驗能量，并對馬爾科夫隨機場采用各向同性的二階馬爾科夫隨機場鄰域，則變化類的先驗能量PEc(i，j)及非變化類PEu(i，j)的先驗能量為尸&("f(c"力，c(p，《))；(9)PEu(i，j)=-8-PEc(i，j)。(10)其中，C(i，j)為像素點(i，j)處的類別，S為C(i，j)的二階馬爾科夫隨機場鄰域，C(p，q)為S中的類別，V(C(i，j)，C(p，q))為鄰域勢函數，通過狄拉克函數計算:f一l，C(,，力-C(m)<formula>formulaseeoriginaldocumentpage7</formula>步驟6，利用權值因子W、變化類和非變化類的似然能量及變化類和非變化類的先驗能量計算非變化類的總能量TEu(i，j)和變化類的總能量TEc(i，j):TEu(i，j)=LEu(i，j)+W(i，j)XPEu(i，j);(12)TEc(i，j)=LEc(i，j)+W(i，j)XPEc(i，j)，(13)若TEu(i，j)<TEc(i，j)，則將像素點(i，j)處的類別更新為非變化類，否則為變化類，得到類別更新后的結果。步驟7，對類別更新后的結果，采用非參數密度估計方法重新估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數得到變化類和非變化類的似然能量，并利用馬爾科夫隨機場重新計算變化類和非變化類的先驗能量。步驟8，重復步驟(6)及步驟(7)直至迭代終止，并存儲每次類別更新后的結果，得到每個像素點的類別更新集合，該迭代終止條件有兩種一種是迭代次數不超過50次，另一種是兩次迭代之間相異的像元數目比例小于給定的閾值T，T=5X10一8。步驟9，利用每個像素點的類別更新集合估計變化類的后驗概率和非變化類的后驗概率，將后驗概率較大的那一類作為該像素點的最終變化檢測結果。本發(fā)明的效果可以通過以下實驗進一步說明本發(fā)明的對比實驗為宋妍和袁修孝(2009)等在文章"基于混合高斯密度模型和空間上下文信息的遙感圖像變化檢測方法及擴展"中提出的變化檢測方法，變化檢測結果的性能采用虛警數、漏檢數及總錯誤數三個指標進行評價。本發(fā)明所設計的三組實驗；第一組為ATM(AirborneThematicM即per)3波段圖像和模擬變化圖像構成的模擬數據集，分別如圖2(a)和圖2(b)所示。其中ATM圖像位于英國Feltwell村莊的一個農田區(qū)，模擬變化圖像是通過模擬地球的天氣變化和電磁波的輻射特性等因素影響并人工地嵌入一些變化區(qū)域得到，圖像大小均為470X335,256灰度級，兩幅圖像的配準誤差為1.5個像元左右。圖2(c)為變化參考圖。對圖2(a)和圖2(b)應用變化矢量分析法得到的差異圖像，如圖2(d)所示。圖2(e)為采用對比實驗方法得到的變化檢測結果，圖2(f)為采用本發(fā)明方法得到的變化檢測結果。第二組為2000年4月和2002年5月的墨西哥郊外的兩幅Landsat7ETM+4波段遙感圖像，分別如圖3(a)和圖3(b)所示。圖像大小均為512X512,256灰度級，圖像配準誤差為1.5個像元左右，變化區(qū)域主要為大火破壞了大面積的當地植被所致，變化參考圖如圖3(c)所示。對圖3(a)和圖3(b)應用變化矢量分析法得到的差異圖像，如圖3(d)所示。圖3(e)為采用對比實驗方法得到的變化檢測結果，圖3(f)為采用本發(fā)明方法得到的變化檢測結果。第三組為1995年9月和1996年7月Landsat-5衛(wèi)星TM(ThematicM即per)傳感器接收的兩幅多光譜圖像，分別如圖4(a)和圖4(b)所示。圖像大小均為300X412,256灰度級。試驗區(qū)為意大利撒丁島包含湖泊的一部分，變化前后湖中水位上升，變化參考圖如圖4(c)所示。對圖4(a)和圖4(b)應用變化矢量分析法得到的差異圖像，如圖4(d)所示。圖4(e)為采用對比實驗方法得到的變化檢測結果，圖4(f)為采用本發(fā)明方法得到的變化檢測結果。表1為第一組實驗結果，從表中可以看出與對比實驗相比，本發(fā)明方法的虛警數減少了1280個像元，漏檢數增加了779個像元，但總的錯誤數減少了501個像元。從圖2(e)和圖2(f)可以看出與對比實驗方法結果相比，本發(fā)明方法減少了孤立噪聲，有效地保持變化區(qū)域的結構信息，總體上說本發(fā)明方法是有效的。表2為第二組實驗結果，從表中可以看出與對比實驗相比，本發(fā)明方法的虛警數減少了623個像元，漏檢數增加了500個像元，但總的錯誤數減少了123個像元。從圖3(e)和圖3(f)可以看出與對比實驗方法結果相比，本發(fā)明方法減少了孤立噪聲，有效地保持變化區(qū)域的結構信息，總體上說本發(fā)明方法是有效的。表3為第三組實驗結果，從表中可以看出與對比實驗相比，本發(fā)明方法的虛警數減少了571個像元，漏檢數增加了153個像元，但總的錯誤數減少了418個像元。從圖4(e)和圖4(f)可以看出與對比實驗方法結果相比，本發(fā)明方法減少了孤立噪聲，有效地保持變化區(qū)域的結構信息，總體上說本發(fā)明方法是有效的。表1第一組實驗結果8<table>tableseeoriginaldocumentpage9</column></row><table>權利要求一種基于非參數密度估計的遙感圖像變化檢測方法，包括如下步驟(1)輸入兩幅不同時相的遙感圖像，并對每幅圖像的每個通道分別進行窗口大小為3×3像素的中值濾波，得到兩時相的去噪后圖像；(2)將去噪后的兩幅圖像應用變化矢量分析，得到一幅差異圖像，并根據該差異圖像計算變權馬爾科夫隨機場的權值因子W；(3)應用K-均值聚類算法將差異圖像聚成變化類和非變化類，得到初始分類結果；(4)利用初始分類結果，采用非參數密度估計方法估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數得到變化類和非變化類的似然能量；(5)對初始分類結果利用馬爾科夫隨機場計算變化類和非變化類的先驗能量；(6)利用權值因子W、變化類和非變化類的似然能量及變化類和非變化類的先驗能量計算變化類的總能量和非變化類的總能量，將總能量較小的那一類作為當前類別，得到類別更新后的結果；(7)對類別更新后的結果，采用非參數密度估計方法重新估計差異圖像中變化類和非變化類的類條件概率密度，再對該類條件概率密度取負自然對數得到變化類和非變化類的似然能量，并利用馬爾科夫隨機場重新計算變化類和非變化類的先驗能量；(8)重復步驟(6)及步驟(7)直至迭代終止，并存儲每次類別更新后的結果，得到每個像素點的類別更新集合，該迭代終止條件為迭代次數不超過50次及兩次迭代之間相異的像元數目比例小于給定的閾值；(9)利用每個像素點的類別更新集合估計變化類的后驗概率和非變化類的后驗概率，將后驗概率較大的那一類作為該像素點的最終變化檢測結果。2.根據權利要求l所述的遙感圖像變化檢測方法，其中步驟(2)所述計算變權馬爾科夫隨機場的權值因子，按如下步驟計算首先，計算每個像素點的特征值，即m=l=1其中P為像素點局部窗口的大小，x(m，n)為局部窗口內每個像素點的灰度值，u(i，j)為局部窗口像素的均值；然后，利用整幅圖像中像素特征值t(i，j)的最大值及最小值將t(i，j)映射到I,Vmax]區(qū)間上，得到每個像素點的權值因子W(i，j)，Vmin=0.5，Vmax=8。3.根據權利要求l所述的遙感圖像變化檢測方法，其中步驟(4)所述采用非參數密度估計方法估計差異圖像中變化類和非變化類的類條件概率密度，通過如下公式進行<formula>formulaseeoriginaldocumentpage2</formula>其中，,(;^I^)和》(;^iA)分別表示非變化類和變化類的類條件概率密度，Sn和sc分別表示非變化類和變化類的像素集合，nn和n。分別表示非變化類和變化類的像素數目，K(□)為高斯核函數，4和H。分別表示非變化類和變化類的自適應窗寬平滑參數Hn=H。(a/Hn-f(Xij)13);Hc=H。(a/Nc-f(Xij)13);其中，H。，a禾Pp均為經驗常數，H。=1，a=40000，P=10。全文摘要本發(fā)明公開了一種基于非參數密度估計的遙感圖像變化檢測方法，主要解決現有技術對差異圖像中與變化類和非變化類相關的統(tǒng)計項的估計存在偏差的問題。其實現過程是輸入兩幅不同時相的遙感圖像，并對每幅圖像的每個通道去噪，得到兩時相的去噪后圖像，并采用變化矢量分析法構造差異影像；應用K-均值聚類算法將差異圖像聚成變化類和非變化類，得到初始分類結果，并采用非參數密度估計的方法估計差異影像中與變化類和非變化類相關的統(tǒng)計項；結合變權馬爾科夫隨機場模型進行自適應的空間約束，得到最終的變化檢測結果。實驗表明本發(fā)明能夠有效地保持圖像的結構信息，并去除孤立噪聲，提高變化檢測處理效率，可用于災情監(jiān)測、土地利用、農業(yè)調查領域。文檔編號G01S7/48GK101694719SQ20091002429公開日2010年4月14日申請日期2009年10月13日優(yōu)先權日2009年10月13日發(fā)明者侯彪,公茂果,劉芳,焦李成,王桂婷,范元章,鐘樺,馬文萍申請人:西安電子科技大學;