楊德賀，陳宜金，楊 溪，呂京國(guó)，張子昕，張 帥

(1.中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083;2.北京建筑工程學(xué)院，北京 100044)

0 引言

多源遙感圖像已應(yīng)用于地震災(zāi)害監(jiān)測(cè)、國(guó)土資源勘查和運(yùn)動(dòng)目標(biāo)識(shí)別等各種領(lǐng)域，而圖像配準(zhǔn)是遙感應(yīng)用的重要基礎(chǔ)。單一傳感器圖像的幾何配準(zhǔn)就是通過(guò)尋找2 景圖像的相同或相似區(qū)域的偏移量，對(duì)待配準(zhǔn)圖像進(jìn)行重采樣，受噪聲和灰度變化的影響比較大。多傳感器獲取的圖像具有互補(bǔ)性和冗余性，包含了更加豐富的信息，可反映地物特點(diǎn)的全面性;因此，多源遙感圖像的幾何配準(zhǔn)是利用相似性測(cè)度計(jì)算圖像間的變換參數(shù)，將同一區(qū)域的多景圖像變換到同一坐標(biāo)系統(tǒng)下，實(shí)現(xiàn)像元的最佳匹配。將配準(zhǔn)后的圖像應(yīng)用于地震等級(jí)評(píng)估中，可解決地震救援和震害評(píng)估對(duì)基礎(chǔ)資料需求的問(wèn)題。地震災(zāi)害信息提取是指利用多源遙感圖像的光譜信息或紋理特征判斷地震災(zāi)害信息變化程度。

在特征定位、變化檢測(cè)和立體視覺(jué)等不同應(yīng)用中，需要進(jìn)行亞像元級(jí)的圖像配準(zhǔn)，而像元級(jí)的配準(zhǔn)算法(例如:空域互相關(guān)配準(zhǔn)［1］、傅里葉變換配準(zhǔn)［2］等)很難滿(mǎn)足亞像元級(jí)圖像配準(zhǔn)的高精度定位要求。目前，學(xué)術(shù)界傾向于研究亞像元級(jí)配準(zhǔn)方法，包括插值方法［3］、微分法［4］和最大互信息法［5］等。插值法取決于插值的質(zhì)量;微分法對(duì)光照的變化敏感;而最大互信息法則需要人工干預(yù)選取特征點(diǎn)。上述方法均不適合多源遙感圖像的自動(dòng)化配準(zhǔn)。

針對(duì)此問(wèn)題，本文提出了多源遙感圖像自動(dòng)配準(zhǔn)算法和分類(lèi)流程，由以下幾個(gè)步驟組成:①提取影像特征點(diǎn)以獲取同名點(diǎn)，采用隨機(jī)抽樣一致性(random sample consensus，RANSAC)算法剔除誤匹配點(diǎn);②通過(guò)對(duì)數(shù)極坐標(biāo)變換，求解旋轉(zhuǎn)和縮放因子;③采用頻域相位相關(guān)算法計(jì)算圖像間的偏移量，利用插值算法對(duì)待配準(zhǔn)圖像重采樣，實(shí)現(xiàn)待配準(zhǔn)圖像與基準(zhǔn)圖像的幾何配準(zhǔn);④利用地震前圖像和與之配準(zhǔn)的地震后圖像進(jìn)行重點(diǎn)區(qū)域特征提取;⑤面向特征信息提取結(jié)果，開(kāi)展基于特征的震害信息提取。

1 算法

1.1 Moravec 特征點(diǎn)提取

Moravec 是點(diǎn)特征提取的經(jīng)典算子［6］，在以每一個(gè)像元為中心的影像窗口(w，w)中，計(jì)算對(duì)角線(xiàn)方向4個(gè)相鄰像元灰度差的平方和作為感興趣值;給定一個(gè)經(jīng)驗(yàn)閾值t(一般取整景圖像的灰度均值)，將感興趣值大于閾值的點(diǎn)作為候選點(diǎn);選取一定大小的窗口，從候選點(diǎn)中選擇一個(gè)最大的感興趣值作為特征點(diǎn)。該算子具有計(jì)算量小、不丟失灰度信息等優(yōu)點(diǎn)，由(w，w)和t 決定特征點(diǎn)的數(shù)量。

RANSAC 算法是一個(gè)以迭代法來(lái)估計(jì)參數(shù)的數(shù)學(xué)模型［7］，應(yīng)用于基于特征的圖像配準(zhǔn)中，具有較高的剔除誤匹配點(diǎn)的能力。首先，它依據(jù)設(shè)定的容許誤差，將所有的初步匹配點(diǎn)對(duì)分為內(nèi)點(diǎn)和外點(diǎn);然后，采用比較準(zhǔn)確的內(nèi)點(diǎn)來(lái)進(jìn)行參數(shù)估計(jì)，剔除誤匹配點(diǎn)。該算法具有較高的可靠性、配準(zhǔn)精度和魯棒性等優(yōu)點(diǎn)。點(diǎn)匹配與誤匹配點(diǎn)剔除的技術(shù)流程如圖1 所示。

圖1 特征點(diǎn)匹配流程Fig.1 Flow chart of matching process for feature points

以地震前全色波段遙感圖像為基準(zhǔn)圖像、地震后高分辨率數(shù)碼圖像為待配準(zhǔn)圖像，采用Moravec特征提取算子進(jìn)行特征點(diǎn)提取(圖2)。

圖2 Moravec 特征點(diǎn)提取效果Fig.2 Extraction of feature points with Moravec

1.2 對(duì)數(shù)極坐標(biāo)變換

對(duì)數(shù)極坐標(biāo)變換是從笛卡爾坐標(biāo)系變換到對(duì)數(shù)極坐標(biāo)系，將圖像匹配中笛卡爾坐標(biāo)系下的旋轉(zhuǎn)和縮放變換轉(zhuǎn)化為對(duì)數(shù)極坐標(biāo)下的平移變換［8］。笛卡爾坐標(biāo)系可表示為

式中:z 為笛卡爾坐標(biāo)空間;(x，y)為點(diǎn)坐標(biāo);ρ 為點(diǎn)的極徑;θ 為點(diǎn)的極角。

對(duì)數(shù)極坐標(biāo)系可表示為

映射變換將圖像在軸向縮放和旋轉(zhuǎn)上的變化轉(zhuǎn)變?yōu)閷?duì)數(shù)極坐標(biāo)上的水平和垂直方向位移，從而實(shí)現(xiàn)從一個(gè)圓形區(qū)域的圖像映射成一個(gè)矩形區(qū)域的圖像，即為對(duì)數(shù)極坐標(biāo)的二維不變性。縮放和旋轉(zhuǎn)可表示為

式中:Δρ 為放大倍數(shù);Δθ 為旋轉(zhuǎn)角度。

水平和垂直方向上的縮放和旋轉(zhuǎn)表示為

圖3 示出對(duì)數(shù)極坐標(biāo)變換的結(jié)果。

圖3 對(duì)數(shù)極坐標(biāo)變換Fig.3 Logarithmic polar coordinate transformation

1.3 相位相關(guān)與插值

首先，利用相位相關(guān)和插值算法對(duì)2 景圖像進(jìn)行插值，以提高圖像分辨率［9］;然后，進(jìn)行傅里葉變換，利用頻域相位相關(guān)的多相位分解求出圖像之間的位移關(guān)系，則可得到圖像之間的亞像元級(jí)配準(zhǔn)結(jié)果。

2 景圖像存在的水平和豎直方向的平移矢量，通過(guò)插值放大一定的倍數(shù)，以解決頻域相位相關(guān)非整數(shù)不能求解的問(wèn)題。對(duì)2 景插值前圖像的歸一化互功率譜進(jìn)行傅里葉逆變換，得到單位脈沖響應(yīng)函數(shù)。它采用二維sinc 函數(shù)來(lái)近似，通過(guò)解求sinc 函數(shù)，求得相位相關(guān)峰值處坐標(biāo)，即可解出平移矢量。常用的插值算法有最鄰近插值法、雙線(xiàn)性插值法和三線(xiàn)性插值法。圖4 示出一維sinc 函數(shù)。

圖4 sinc 二維空間(左)和三維空間顯示(右)Fig.4 Display of sinc in 2D (left)and 3D(right)

1.4 基于特征的圖像分類(lèi)

應(yīng)用于地物分類(lèi)的圖像含有不同類(lèi)型的地物，需要選擇一定的特征代表不同類(lèi)型的地物;以一定的方法將特征所對(duì)應(yīng)的空間劃分為不同的子空間，使圖像中代表各類(lèi)地物的像元?dú)w為相對(duì)應(yīng)的子空間。因此，在進(jìn)行遙感圖像分類(lèi)時(shí)，既要考慮影像的灰度信息，也要顧及影像的結(jié)構(gòu)信息，故需要選擇不同的影像特征參數(shù)來(lái)反映這些要求。目前，常用的影像特征參數(shù)包括均值、方差、能量、對(duì)比度、熵、相關(guān)性和均質(zhì)度等［10］。采用光譜信息與結(jié)構(gòu)信息等多種特征相結(jié)合的方法，可以充分利用影像灰度及其分布的信息，改善圖像分類(lèi)的效果。

2 實(shí)驗(yàn)

2.1 數(shù)據(jù)源

實(shí)驗(yàn)數(shù)據(jù)為四川省玉樹(shù)地區(qū)2009年5月Quick-Bird 全色波段遙感圖像和2010年6月ADS40 機(jī)載高分辨率數(shù)字圖像。選取其中1個(gè)波段的圖像進(jìn)行幾何配準(zhǔn)實(shí)驗(yàn)，配準(zhǔn)精度為0.1個(gè)像元;并進(jìn)行灰度變化和抗噪性實(shí)驗(yàn)，測(cè)試兩者對(duì)配準(zhǔn)結(jié)果的影響。實(shí)驗(yàn)證明，基于傅里葉變換的配準(zhǔn)方法具有較高的魯棒性和較強(qiáng)的抗噪性。在旋轉(zhuǎn)和縮放的情況下，本文提出的自動(dòng)配準(zhǔn)方法具有較高的穩(wěn)定性。分類(lèi)后的變化檢測(cè)實(shí)驗(yàn)表明，穩(wěn)定的幾何配準(zhǔn)對(duì)信息變化檢測(cè)有較高的魯棒性。圖5 示出本次實(shí)驗(yàn)的流程。

圖5 實(shí)驗(yàn)流程Fig.5 Flow chart of experiment

2.2 配準(zhǔn)實(shí)驗(yàn)

以震前QuickBird 全色圖像為基準(zhǔn)圖像(圖6(a))，以震后ADS40 航攝圖像為待配準(zhǔn)圖像(圖6(b))，使用本文提出的自動(dòng)配準(zhǔn)方法進(jìn)行多源遙感圖像的幾何配準(zhǔn)實(shí)驗(yàn)，獲得與QuickBird 全色圖像配準(zhǔn)的ADS40 航攝圖像(圖7)，用于配準(zhǔn)的特征點(diǎn)均方根誤差RMSE=0.96 像元(表1)。

圖6 地震前QuickBird(左)和地震后ADS40(右)圖像Fig.6 QuickBird image before earthquake (left)and ADS40 image after earthquake (right)

圖7 配準(zhǔn)后的ADS40 圖像Fig.7 ADS40 image after registration

表1 部分特征點(diǎn)誤差Tab.1 Error of some feature points (像元)

不同時(shí)相的多源遙感圖像，因受大氣、太陽(yáng)入射角、視角及地形的影響，導(dǎo)致地震前后影像灰度變化大，影響到特征匹配的效果，所以在幾何配準(zhǔn)前要對(duì)地震后圖像進(jìn)行相對(duì)輻射校正(即輻射歸一化)。經(jīng)灰度調(diào)整后，采用Moravec 對(duì)2 景圖像提取匹配的特征點(diǎn);利用RANSAC 剔除誤匹配點(diǎn)對(duì)，對(duì)匹配的特征點(diǎn)對(duì)進(jìn)行對(duì)數(shù)極坐標(biāo)變換;解求變換點(diǎn)對(duì)的旋轉(zhuǎn)角和縮放系數(shù)，進(jìn)行相位相關(guān)求解偏移量;采用雙線(xiàn)性插值重采樣完成配準(zhǔn)過(guò)程。

2.3 灰度變化和噪聲對(duì)配準(zhǔn)的影響

通過(guò)灰度線(xiàn)性變化函數(shù)D=aT+b(其中a，b 為灰度變化因子)，使原來(lái)的灰度值T 變?yōu)镈。由于噪聲對(duì)圖像的配準(zhǔn)存在影響，故利用衡量噪聲強(qiáng)弱的信噪比的對(duì)數(shù)函數(shù)來(lái)加入噪聲。圖8 示出灰度變化后和加入椒鹽噪聲后的圖像。

圖8 灰度變化后(左)和加入椒鹽噪聲后(右)圖像Fig.8 Images after gray change(left)and adding noise(right)

從表2 可以看出，灰度變化和噪聲對(duì)配準(zhǔn)誤差的影響都比較小。

表2 灰度變化和噪聲對(duì)配準(zhǔn)誤差的影響Tab.2 Influence of gray change and noise to registration error

2.4 配準(zhǔn)誤差對(duì)分類(lèi)精度的影響

對(duì)自動(dòng)配準(zhǔn)后的2 期遙感圖像進(jìn)行了分類(lèi)比較。結(jié)果表明，當(dāng)配準(zhǔn)誤差較大時(shí)，分類(lèi)面積的平均變化率會(huì)較大，但也會(huì)穩(wěn)定在一定范圍;當(dāng)配準(zhǔn)誤差較小時(shí)，分類(lèi)面積的平均變化率較小，分類(lèi)精度較高。表3 示出不同配準(zhǔn)誤差下的地物分類(lèi)變化率。

表3 不同配準(zhǔn)誤差下的分類(lèi)變化率Tab.3 Change rate of classification under different registration errors

2.5 震害信息提取

本文采用分類(lèi)后處理的方式提取地震災(zāi)害信息。針對(duì)所提取的重點(diǎn)目標(biāo)區(qū)域，以多特征統(tǒng)計(jì)向量為分類(lèi)標(biāo)準(zhǔn)。首先對(duì)地震前后的重點(diǎn)區(qū)域分別進(jìn)行基于特征的分類(lèi)，然后根據(jù)2 期圖像的分類(lèi)結(jié)果進(jìn)行變化檢測(cè)(圖9)。

圖9 震害信息提取Fig.9 Extraction of earthquick damage informatiom

震害信息的提取結(jié)果詳見(jiàn)表4。

表4 震害信息提取結(jié)果Tab.4 Extraction results of earthquick damage information (像元)

從圖9 和表4 可以看出，地震前、后房屋、裸地和道路3 種類(lèi)別的區(qū)分明顯，分類(lèi)效果較好;地震后的總體分類(lèi)精度為68%，說(shuō)明地震后的分類(lèi)結(jié)果與地震前的實(shí)際類(lèi)型相差較大;地震后的裸地類(lèi)別變化相對(duì)較小，地震后道路和房屋變化相對(duì)較大。盡管研究區(qū)域中的房屋損毀情況比較大，但是對(duì)房屋的分類(lèi)較為準(zhǔn)確。震害信息變化圖(圖9(c))表明，對(duì)多源遙感圖像進(jìn)行幾何配準(zhǔn)、分類(lèi)后，能夠檢測(cè)出房屋、裸地、道路等不同類(lèi)別的損毀程度，可以為震害損毀分析與評(píng)價(jià)提供重要信息。

3 結(jié)論

本文提出了多源遙感圖像自動(dòng)配準(zhǔn)的算法，并通過(guò)震害信息檢測(cè)對(duì)配準(zhǔn)的效果進(jìn)行了驗(yàn)證。結(jié)果表明:

1)解決基準(zhǔn)圖像和待配準(zhǔn)圖像的旋轉(zhuǎn)、縮放以及配準(zhǔn)自動(dòng)化問(wèn)題是當(dāng)前多源遙感圖像幾何配準(zhǔn)的難點(diǎn)。幾何配準(zhǔn)精度的高低也決定了圖像分類(lèi)效果的好壞，并會(huì)對(duì)變化檢測(cè)的結(jié)果產(chǎn)生影響。

2)本文提出的面向地震災(zāi)害信息提取的多源遙感圖像自動(dòng)配準(zhǔn)算法，經(jīng)遙感圖像數(shù)據(jù)的幾何配準(zhǔn)實(shí)驗(yàn)證明，該算法綜合了Moravec 算子、RANSAC算子、對(duì)數(shù)極坐標(biāo)變換和相位相關(guān)等算法的優(yōu)點(diǎn)，實(shí)現(xiàn)了多源遙感像亞像元級(jí)的自動(dòng)配準(zhǔn)，較好地解決了灰度變化和噪聲對(duì)幾何配準(zhǔn)的影響以及基準(zhǔn)圖像和待配準(zhǔn)圖像的旋轉(zhuǎn)和縮放等問(wèn)題。

3)將用本文方法配準(zhǔn)后的震前和震后圖像應(yīng)用于地震災(zāi)害信息提取，震害類(lèi)型變化信息統(tǒng)計(jì)結(jié)果表明，較高的幾何配準(zhǔn)精度提高了震害的分類(lèi)精度，說(shuō)明本文提出的高精度自動(dòng)配準(zhǔn)方法和基于特征的分類(lèi)可以較好地應(yīng)用于地震震害評(píng)估、海洋水色反演等領(lǐng)域的實(shí)時(shí)變化檢測(cè)中。

本文的不足之處為:①因獲取遙感圖像的時(shí)相及傳感器等方面的不同，會(huì)出現(xiàn)實(shí)際未變化地物卻存在灰度變化的情況，所以本文的幾何配準(zhǔn)方法研究應(yīng)采用相對(duì)輻射校正處理后的遙感圖像;②本文涉及的算法比較復(fù)雜，下一步工作將研究如何運(yùn)用并行計(jì)算方法提高計(jì)算效率。

［1］吳金鹿，陳紹煒，孫 磊.一種基于互相關(guān)匹配的空域目標(biāo)提取方法［J］.科學(xué)技術(shù)與工程，2011，11(32):7948－7951.Wu J L，Chen S W，Sun L.An airspace object extraction method based on cross－correlation matching［J］.Science Technology and Engineering，2011，11(32):7948－7951.

［2］林 茂.基于改進(jìn)的曲線(xiàn)傅里葉變換圖像配準(zhǔn)研究［J］.計(jì)算機(jī)仿真，2011，28(10):265－268.Lin M.Improved digital image contour curve Fourier descriptor algorithm［J］.Computer Simulation，2011，28(10):265－268.

［3］詹 毅.基于泰勒展開(kāi)式的圖像插值方法［J］.計(jì)算機(jī)工程，2012，38(13):202－204.Zhan Y.Image interpolation method based on Tyler expansion［J］.Computer Engineering，2012，38(13):202－204.

［4］張曉芬.基于偏微分方程的圖像配準(zhǔn)技術(shù)研究［D］.河北:華北電力大學(xué)，2009:32－39.Zhang X F.Research on image registration technique based on PDE’s［D］.Hebei:North China Electric Power University，2009:32－39.

［5］程有娥.基于角點(diǎn)特征和最大互信息的圖像配準(zhǔn)［J］.計(jì)算機(jī)系統(tǒng)應(yīng)用，2012，21(6):186－190.Cheng Y E.Image registration based on mutual information and corner points［J］.Computer Systems and Applications，2012，21(6):186－190.

［6］劉曉龍，李英成.基于多源遙感影像融合的影像匹配技術(shù)［J］.測(cè)繪科學(xué)，2007，32(3):59－61.Liu X L，Li Y C.Image matching based on multi－source remotesensing image fusion［J］.Science of Surveying and Mapping，2007，32(3):59－61.

［7］周劍軍，歐陽(yáng)寧，張 彤，等.基于RANSAC 的圖像拼接方法［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2009，30(24):5692－5694.Zhou J J，Ouyang N，Zhang T，et al.Image mosaic method based on RANSAC［J］.Computer Engineering and Design，2009，30(24):5692－5694.

［8］雷 凱，劉艷瀅，王延杰，等.基于特征點(diǎn)的對(duì)數(shù)極坐標(biāo)變換圖像配準(zhǔn)算法［J］.光學(xué)技術(shù)，2007，32(3):435－437.Lei K，Liu Y Y，Wang Y J，et al.Log polar transformation based on feature points for image registration［J］.Optical Technique，2007，32(3):435－437.

［9］劉衛(wèi)光，崔江濤，周利華.插值和相位相關(guān)的圖像亞像素配準(zhǔn)方法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，17(6):1273－1277.Liu W G，Cui J T，Zhou L H.Subpixel registration based on interpolation and extension phase correlation［J］.Journal of Computer－aided Design and Computer Graphics，2005，17(6):1273－1277.

［10］賈永紅.數(shù)字圖像處理［M］.武漢:武漢大學(xué)出版社，2003:182－186.Jia Y H.Digital image processing［M］.Wuhan:Wuhan University Press，2003:182－186.