專利名稱：基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法
技術領域：
本發明屬于光學遙感圖像處理技術領域，涉及同一多光譜傳感器在不同時間獲取的相同地理位置的遙感圖像間存在較大平移和旋轉錯位情況下的配準，可用于多時相多光譜遙感圖像變化檢測、融合、鑲嵌等方面的前期配準工作。
背景技術：
近年來，隨著多光譜遙感圖像數據的公開化程度和處理技術水平的不斷提高，多光譜遙感圖像數據已經被廣泛用于城市規劃、資源調查、自然災害預防、農作物長勢監測、軍事目標探查等多種應用領域。多光譜遙感圖像的廣泛應用首先需要進行圖像數據的幾何校正、輻射校正、配準、去噪、分割、鑲嵌、融合、目標檢測、變化檢測等多項數字圖像處理技術。其中，在進行多光譜遙感圖像變化檢測、融合以及鑲嵌等工作前，首先需要對不同傳感器獲取的同一地區的圖像或相同傳感器在不同時間獲取的同一地區的圖像進行配準，以消除因獲取圖像的時間、角度、環境和傳感器成像機理的不同造成圖像間的平移、旋轉、伸縮及局部形變問題。1992 年 L.G.Brown 在文獻”A Survey of Image Registration Techniques'ACMComputing Surveys, 1992, vol.24, n0.4, pp.325-376)中指出，所有的圖像配準技術都是由特征空間、搜索空間、搜索策略和相似性測度四個方面構成。特征空間指圖像中用于進行搜索匹配的信息集合，是搜索匹配的對象；而每一次搜索的匹配程度是以相似性測度為評價標準的，因而相似性測度直接關系到配準的精確度；搜索空間是搜索匹配得到的所有圖像變換參數所組成的空間；搜索策略則是搜索最優匹配參數時所采取的優化策略。多光譜遙感圖像配準的過程一般為首先確定特征空間和搜索空間，然后以特征空間為基礎設計并采用某種相似性測度，進而按一定的搜索策略搜索最佳匹配變換參數，并將其帶入圖像變換模型從而達到配準的目的。根據特征空間的不同，現有的多光譜遙感圖像配準方法通常可以分為基于像素灰度特征和基于圖像幾何特征兩大類。大部分基于像素灰度特征的方法一般直接利用整幅圖像的灰度信息，通過建立某種像素間的相似性測度來衡量兩幅圖像重疊部分地表反射屬性的匹配程度，進而尋找到最優匹配時的平移、旋轉和伸縮等變換參數。常見的灰度相似性測度如平均平方差MSD、平均絕對差MAD、歸一化互相關NCC等通常要求待配準圖像對應像素的灰度值具有線性關系，因而在利用它們進行不同時相的圖像配準時容易受到噪聲和圖像間整體亮度差異的影響，而對于不同傳感器或不同波段圖像的配準則基本不適用。最初作為多模態醫學圖像配準相似性測度的最大互信息準則由于不要求待配準圖像對應像素的灰度信息具有嚴格的線性關系，而認為只要兩幅圖像中主要的幾何結構特征對準時即達到最佳匹配，因此已被用于不同時相和不同波段圖像間的配準。在對不同時相的多光譜遙感圖像的配準中，由于基于像素灰度特征的配準方法在搜索最佳匹配變換參數時通常是由全體像素共同決定圖像間的對準程度，計算量大，并且難免會受到圖像中不同位置大小的云、陰影及變化地物的干擾，造成配準精度和穩定性的下降。
基于圖像幾何特征的配準方法則首先需要對待配準圖像提取角點、邊緣等明顯的幾何特征，然后進行對應特征對或特征集之間的搜索匹配，進而得到圖像變換參數。基于邊緣特征的配準方法常見的如利用一些傳統空域邊緣檢測算子進行邊緣特征提取，但往往由于提取的邊緣過于瑣碎，不能很好地體現圖像中的結構特征，對后續匹配精度和速度造成影響，因而需要對提取到的邊緣進行優化處理。近年來，小波變換wavelet、曲波變換curvelet和輪廓波變換contourlet等多尺度幾何分析方法由于在描述和分析圖像中的邊緣特征時具有優異的性能，也常被用于邊緣特征的提取。1997年J.W.Hsieh.H.Y.M.Liao、K.C.Fan 等在文獻，，Image registration using a new edge-based approach，，(ComputerVision and Image Understanding, 1997，vol.67，n0.2，pp.112-130)中首次提出利用小波變換得到系數的模局部極大值進行邊緣特征點提取，然后以互相關作為相似性測度進行邊緣特征點對的搜索匹配，其中還使用了一種非迭代的策略消除誤匹配點對，達到了較好的配準效果。2008年孫淑一.吳勇和吳建民在文獻”一種基于邊緣特征的圖像配準方法”(計算機工程與應用，2008，vol.44，n0.7，pp.94-96)中、2010年孫雅琳和王菲在文獻”基于邊緣和互信息的紅外與可見光圖像配準”(電子科技，2010, vol.23，n0.4，pp.80-82)中也都采用了相同的邊緣特征提取方法，不同的是，這兩種方法在搜索匹配時前者計算了兩幅邊緣圖像重疊部分的交互方差，而后者計算了兩幅邊緣圖像重疊部分的歸一化互信息。如果將這兩種方法用于不同時相多光譜遙感圖像的配準，則容易受到圖像中變化部分對應特征點的干擾而造成誤配準。2009年陳志剛、尹福昌和孫孚在文獻”基于非采樣Contourlet變換高分辨率遙感圖像配準”(光學學報，2009, vol.29，n0.10，pp.2744-2750)中提出首先由圖像通過非下采樣Contourlet變換得到的各方向子帶計算得到梯度矢量模，閾值分割得到邊緣特征點后由歸一化互相關匹配搜索匹配控制點，并利用概率支持算法消除誤匹配點對，降低了變化部分對應特征點錯誤匹配的概率。總體來說，由于不直接采用像素的灰度值作為相似性測量的對象，因而該類方法受不同時相間的亮度差異、不同傳感器以及不同波段圖像間的成像特性差異的影響較小，適合多時相多光譜遙感圖像的配準，但要求待配準圖像中含有比較明顯的邊緣信息，否則提取不到理想的邊緣將很可能導致配準失敗。另外，其中邊緣特征提取和匹配方法的準確性及穩定性也是影響其配準性能的主要因素
發明內容
本發明的目的在于針對上述已有技術的不足，提出一種基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，以消除不同時相圖像間存在的整體亮度差異和變化部分對配準精度及穩定性造成的影響，提高配準的精度和穩定性。為實現上述目的，本發明的技術方案包括如下步驟:(I)輸入在兩個時相獲取的同一地區的多光譜圖像集IdP I2，其中，時相I的波段序列圖像集I1=IA11I，時相2的波段序列圖像集I2= {A2b}，Atb為兩個圖像集中的每一幅單波段圖像，其中，上標b表示波段序號，b=l, 2，…，B, B為總波段數，下標t為時相序號，t={l, 2}，每一幅單波段圖像Atb均由η行η列像素構成，η的取值為32的倍數；(2)對兩時相多光譜圖像集I1和I2分別進行窗口大小為3X3像素的中值濾波去噪,并歸一化處理,得到歸一化圖像集IjPf2;
(3)將歸一化圖像集Ijpf2中的每一幅圖像均進行多層二維離散Haar小波變換，
得到第I層的低頻系數矩陣為AP、垂直方向的高頻系數矩陣為V、水平方向的高頻系數
矩陣為對角線方向的高頻系數矩陣為D:,其中分解層數1=1，...，L，L為最大分解層數，其取值范圍為I 4 ;(4)將時相I和時相2的所有B個波段圖像的第L層垂直方向的高頻系數矩陣V:、
水平方向的高頻系數矩陣和對角線方向的高頻系數矩陣Dfi取絕對值后空間位置對應元素相加，分別得到時相I和時相2的高頻系數圖像HVD1和HVD2 ；(5)對高頻系數圖像HVD1和HVD2采用最大類間方差法分別進行閾值分割，得到時相I和時相2的邊緣圖像0和 2;(6)對邊緣圖像G和02，統計每一條不與其他邊緣有8鄰域連接的邊緣線段在8鄰域下所包含的像素個數，并將其中像素個數小于長度閾值Tc的邊緣線段刪除，對應得到
時相I和時相2的強邊緣圖像C1和C2，其中:
權利要求
1.一種基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，包括步驟如下: (1)輸入在兩個時相獲取的同一地區的多光譜圖像集I1和I2，其中，時相I的波段序列圖像集I1=IA11I，時相2的波段序列圖像集I2= {A2b}，Atb為兩個圖像集中的每一幅單波段圖像，其中，上標b表示波段序號，b=l, 2，…，B, B為總波段數,下標t為時相序號，t={l, 2},每一幅單波段圖像Atb均由η行η列像素構成，η的取值為32的倍數； (2)對兩時相多光譜圖像集I1和I2分別進行窗口大小為3X3像素的中值濾波去噪，并歸一化處理，得到歸一化圖像集IjPf2; (3)將歸一化圖像集^和12中的每一幅圖像均進行多層二維離散Haar小波變換，得到第I層的低頻系數矩陣為Ag、垂直方向的高頻系數矩陣為Vfj、水平方向的高頻系數矩陣為對角線方向的高頻系數矩陣為1^，其中分解層數1=1，...，L，L為最大分解層數，其取值范圍為I 4 ; (4)將時相I和時相2的所有B個波段圖像的第L層垂直方向的高頻系數矩陣%、水平方向的高頻系數矩陣Hf,、對角線方向的高頻系數矩陣私取絕對值后空間位置對應元素相加，分別得到時相I和時相2的高頻系數圖像HVD1和HVD2 ； (5)對高頻系數圖像HVD1和HVD2采用最大類間方差法分別進行閾值分割，得到時相I和時相2的邊緣圖像 和0 (6)對邊緣圖像C,和62,統計每一條不與其他邊緣有8鄰域連接的邊緣線段在8鄰域下所包含的像素個數，并將其中像素個數小于長度閾值Tc的邊緣線段刪除，對應得到時相I和時相2的強邊緣圖像C1和C2，其中%:1.IUu (7)將強邊緣圖像C1和C2進行旋轉和平移匹配，并計算每次匹配的邊緣對齊度Si (a, q，P)，得到邊緣對齊度集合SI= {Si (a，q，P) }，其中a為旋轉角度，q為橫移參量，P為縱移參量，a=0, 1,2, -,180, q=0, 1,2,…，n/2L+l，p=0，1,2,...，n/2L+l ； (8)用邊緣對齊度集合SI中的最大值所對應的旋轉角度aM、橫移參量qM和縱移參量pM對時相I的各波段歸一化圖像I/均進行相應的旋轉和平移，得到時相I的各波段粗配準圖像ΑΛ將時相I的所有波段粗配準圖像的集合與步驟(2)得到的時相2的歸一化圖像集 2合并，得到粗配準圖像集合I; (9)對粗配準圖像集合I中兩個時相各波段圖像分別采用對數殘差修正模型將像素的灰度值轉換為相對地物光譜反射率值，得到相對地物光譜反射率圖像集合R ; (10)對時相I的各波段相對地物光譜反射率圖像R113,分別進行I層二維離散Haar小波變換，得到對應波段的低頻系數矩陣ARb、垂直方向的高頻系數矩陣VRb、水平方向的高頻系數矩陣HRb、對角線方向的高頻系數矩陣DRb ； (11)將時相I各波段低頻系 數矩陣ARb中的元素全部置為零，并與高頻系數矩陣VRb、HRb和DRb進行二維離散Haar小波逆變換，得到時相I的各波段重構圖像； (12)將時相I所有B個波段的重構圖像取絕對值后對應像素值相加，得到一幅模糊邊緣圖像； (13)對時相I的模糊邊緣圖像采用最大類間方差法進行閾值分割，得到粗配準邊緣圖像CR; (14)將粗配準邊緣圖像CR從第I行第I列像素開始，由左至右、由上至下分為互不重疊且大小均為K列K行的uXu個子圖像CRij，其中i和j分別為粗配準邊緣子圖像的行塊序號和列塊序號，i=l，2，...，u，j=l，2，…，u，u為粗配準邊緣圖像CR分塊后的水平方向和垂直方向的子圖像數目，且滿足n=KX U，K為偶數，40彡K彡100 ； (15)將時相I的各波段相對地物光譜反射率圖像R113中與粗配準邊緣子圖像CRij相同空間位置范圍的圖像塊記為時相I的各波段相對地物光譜反射率子圖像.并將時相I的所有波段相對地物光譜反射率子圖像丨丨與時相2的所有波段相對地物光譜反射率圖像IR2bI采用光譜反射率曲線差異最小方法進行匹配，得到時相2的相對地物光譜反射率圖像集{R2}與時相I的相對地物光譜反射率子圖像{RUj}相匹配的光譜反射率曲線總差異量DRCiP匹配旋角參量Θ mi」、匹配橫移參量Qmij以及匹配縱移參量Pmij ； (16)將時相2的相對地物光譜反射率圖像集{RJ與時相I的相對地物光譜反射率子圖像{R j}相匹配的光譜反射率曲線總差異量DRCij按行塊序號和列塊序號的順序作為對應行列的元素構成矩陣，利用最大類間方差法對該矩陣進行閾值分割，得到差異二值矩陣Duc ; (17)分別統計差異二值矩陣du。中O值元素對應的匹配旋角參量emu、匹配水平參量Qmij以及匹配垂直參量Pmij的數目，將各參量的數目最大值記為最優旋角參量Θ mm、最優水平參量Q_以及最優垂直參量P_ ； (18)用最優旋角參量Θ_、最優水平參量Q_以及最優垂直參量P_對粗配準圖像集合！中時相2的各波段歸一化圖像進行相應的旋轉和平移，得到時相2的各波段最終配準圖像 (19)將時相2的所有波段最終配準圖像{A/j與步驟(8)中得到的時相I的所有波段粗配準圖像丨A。合并，構成最終配準圖像集合L
2.根據權利要求1所述的基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，其中步驟(7)所述的將強邊緣圖像C1和C2進行旋轉和平移匹配，并計算每次匹配的邊緣對齊度Si (a, q, p),通過如下步驟進行: (7a)將強邊緣圖像C1和C2相匹配的旋轉角度a、橫移參量q和縱移參量p的初始值均賦為O ;將強邊緣圖像C1的行序號和列序號均為ζ的像素點記為Ce，以像素點Ce為中心點O，水平方向為X軸、垂直方向為Y軸建立XOY直角坐標系，其中ζ=η2 +1 ； (7b)將強邊緣圖像C1在XOY坐標系中以O點為中心逆時針旋轉a-90度，采用最近鄰線性插值方法進行圖像插值，得到旋轉圖像Cla ； (7c)在XOY坐標系中，對橫坐標和縱坐標范圍均在[-2 ζ，2 ζ ]內且旋轉圖像Cla以外的像素點的值賦為0，得到大小為(4 ζ +1) X (4 ζ +1)像素的旋轉擴展圖像CEla ； (7d)將強邊緣圖像C2與旋轉擴展圖像CEla的第1+p行至第η/2ι+ρ行、第Ι+q列至第n/2L+q列范圍內的像素進行對應空間位置的像素值與運算，將旋轉擴展圖像CEla的行序號在[Ι+ρ,η/^+ρ]范圍外且列序號在[l+q’n/^+q]范圍外的像素值賦為O，得到一幅與旋轉擴展圖像CEla大小相同的邊緣交集圖像CAaqp ； (7e)在邊緣交集圖像CAaqp中，對任意一條不與其他邊緣有8鄰域連接的邊緣線段，計算該邊緣線段中包含的像素個數Lc與長度閾值Tc的比值Lc/Tc，若該比值小于1，則將該邊緣線段包含的所有像素的值都賦為1，否則，賦值為Lc/Tc ； (7f)對邊緣交集圖像CAaqp中的每一條獨立的邊緣線段都執行步驟(7e)，得到增強邊緣交集圖像CA' aqp ； (7g)將增強邊緣交集圖像CA' _中所有像素的值求和，得到邊緣對齊度Si (a, q, p)； (7h)若橫移參量q〈2 ζ +1，則縱移參量P保持不變，對橫移參量q的值加1，返回步驟(7d);若橫移參量q=2(+l且縱移參量ρ〈2ζ+1，則將橫移參量q的值置為O，縱移參量p的值加1，返回步驟(7d);若橫移參量q=2ζ+l，縱移參量p=2ζ+l且旋轉角度a〈180，則將旋轉角度a的值加3，橫移參量q和縱移參量P的值都置為O，返回步驟(7b);若橫移參量q=2 ζ +1，縱移參量ρ=2 ζ +1，且旋轉角度a=180，則強邊緣圖像C1和C2的旋轉和平移匹配完成，得到所有的邊緣對齊度{Si(a，q，p)}。
3.根據權利要求1所述的基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，其中步驟(9)所述的對粗配準圖像 集合I中兩個時相各波段圖像分別采用對數殘差修正模型將像素的灰度值轉換為相對地物光譜反射率值，得到相對地物光譜反射率圖像集合R，通過如下步驟進行: (9a)將時相I的第b波段粗配準圖像A,中像素點(x，y)的灰度值記為/ΤΛυ),再將轉換為相對地物光譜反射率值Ab(Xj)的對數:
4.根據權利要求1所述的基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，其中步驟(15)所述的將時相I的所有波段相對地物光譜反射率子圖像與時相2的所有波段相對地物光譜反射率圖像{R2b}采用光譜反射率曲線差異最小方法進行匹配，得到時相2的相對地物光譜反射率圖像集{R2}與時相I的相對地物光譜反射率子圖像{Rnj}相匹配的光譜反射率曲線總差異量DRCij、匹配旋角參量Qmij、匹配橫移參量Qmij以及匹配縱移參量Pmu，通過如下步驟進行: (15a)將粗配準邊緣子圖像CRij的行塊序號i和列塊序號j初始值賦為2，CRu的中心點為該子圖像的第V行第V 列的像素，其中V=Ceil (K/2)，ceil (.)表示取大于或等于括號內數值的最小整數，K為粗配準邊緣子圖像CRij的總行數或總列數； (15b)若粗配準邊緣子圖像CRij中像素值為I的邊緣像素點的總數小于λ ΧΚ，執行步驟(15c)，否則，將時相2的相對地物光譜反射率圖像集{R2}與時相I的相對地物光譜反射率子圖像{RUj}相匹配的旋角參量Θ、水平參量Q以及垂直參量P初始值賦為O,執行步驟(15d)，其中邊緣約束系數λ的取值范圍為2彡λ < 4 ; (15c)若列塊序號j〈u-l時，則行塊序號i保持不變，列塊序號j的值加1，返回步驟(15b);若列塊序號j=u-l且行塊序號i〈u_l時，則行塊序號i的值加1，列塊序號j的值賦為2，返回步驟(15b);若列塊序號j=u-l且行塊序號i=u-l時，則將所有行塊序號i=l或i=u-l,或列塊序號j=l或j=u-l的粗配準邊緣子圖像CRij的光譜反射率曲線總差異量DRCij置為空值，執行步驟(16); (15d)將時相2的各波段相對地物光譜反射率圖像R2b中與粗配準邊緣子圖像CRij中心點空間位置對應的像素點記為Ro，其橫坐標為(1-l)XK+v，縱坐標為(j-l)XK+v;以像素點Ro為中心點，分別將時相2的各波段相對地物光譜反射率圖像R2b逆時針旋轉Θ -5度，采用最近鄰線性插值方法進行圖像插值，得到各波段的相對地物光譜反射率旋轉圖像在該圖像中將以像素點Ro為中心、大小為(3v+l) X (3v+l)像素的圖像塊作為時相2的各波段窗圖像 (15e)分別將時相I的各波段相對地物光譜反射率圖像R113中與粗配準邊緣子圖像CRij相同空間位置范圍的圖像塊作為時相I的各波段搜索圖像64; (15f)將時相I的各波段搜索圖像中的第I行第I列像素點與時相2的各波段窗圖像中的第1+Q行第1+P列像素點空間位置對應，則時相I的各波段搜索圖像^^中的像素點(X' ,1')與時相2的各波段窗圖像中的像素點(X' +Q，y' +P)空間位置對應，計算時相I的各波段搜索圖像Ι 〗,與時相2的各波段窗圖像食 ,之間的光譜反射率曲線差異量 DRC(0，Q，p):
全文摘要
本發明公開了一種基于邊緣和光譜反射率曲線的多時相遙感圖像配準方法，主要解決現有技術對存在變化區域和整體亮度差異的兩時相多光譜圖像集進行配準時穩定性不足、精度低的問題。其實現步驟為對輸入兩時相多光譜圖像集進行預處理后的同時相的小波高頻系數矩陣相加并分割，得到兩時相強邊緣圖像；計算兩時相強邊緣圖像的邊緣對齊度的配準參數，以對兩時相圖像集進行配準，得到粗配準圖像集；將粗配準圖像集轉換為相對地物光譜反射率圖像集，對該圖像集中的時相1的小波高頻系數重構得到粗配準邊緣圖像，用該圖像對時相1和時相2的相對地物光譜反射率圖像進行局部匹配，得到最終配準圖像集。本發明可用于變化檢測、融合、及鑲嵌。
文檔編號G06T5/00GK103198483SQ20131011788
公開日2013年7月10日 申請日期2013年4月7日 優先權日2013年4月7日
發明者王桂婷, 焦李成, 孫一博, 公茂果, 鐘樺, 王爽, 張小華 申請人:西安電子科技大學