張晗 謝妮慧 于飛 楊進(jìn)

（北京空間機(jī)電研究所，北京 100094）

0 引言

隨著星載遙感相機(jī)的分辨率由1～2m向著0.1m量級(jí)的不斷提高，衛(wèi)星抖動(dòng)已經(jīng)成為制約高分辨率遙感成像的關(guān)鍵瓶頸[1]。國(guó)內(nèi)目前像移估計(jì)在實(shí)際應(yīng)用的場(chǎng)景精度均不高，且多應(yīng)用于車載成像系統(tǒng)，估計(jì)精度達(dá)到0.5個(gè)像素。應(yīng)用于航天遙感的，如：“高分五號(hào)”衛(wèi)星對(duì)日成像跟蹤控制器，估計(jì)精度僅1個(gè)像素[2-3]。而針對(duì)高分辨率遙感圖像的清晰化處理，像移估計(jì)精度要達(dá)到0.1個(gè)像素，如圖1所示。本文將結(jié)合現(xiàn)有像移估計(jì)算法，分析影響算法估計(jì)精度的若干條件，獲得高精度像移估計(jì)實(shí)現(xiàn)方案。

圖1 不同精度下的圖像清晰度Fig.1 Image sharpness with different precision values

比較經(jīng)典的圖像像移估計(jì)算法是塊匹配法[4-5]，它具有精度高、實(shí)現(xiàn)簡(jiǎn)單的優(yōu)點(diǎn)，但其計(jì)算量龐大，所以無法用于實(shí)時(shí)性要求較高的場(chǎng)合[6]?；叶韧队胺╗7]是對(duì)塊相關(guān)法的改進(jìn)，雖然實(shí)現(xiàn)效率高，但其測(cè)量誤差會(huì)隨著圖像的位移增加而增大[8-9]。相位相關(guān)法[10]利用圖像的頻域信息，其相關(guān)峰具有單峰值特性，很容易檢測(cè)，且其計(jì)算精度與運(yùn)算量和圖像的位移大小無關(guān)[11]，本文即采用此種方式來實(shí)現(xiàn)穩(wěn)像。

相位相關(guān)算法必須作用于高速CCD相機(jī)拍攝的連續(xù)幀圖像。主要原理如圖2所示。

實(shí)際CCD相機(jī)與主成像相機(jī)同軸同步對(duì)地成像，但主成像相機(jī)所成圖像為遙感圖像，成像尺寸大、曝光時(shí)間長(zhǎng)，衛(wèi)星在成像期間的抖動(dòng)無法被準(zhǔn)確記錄[12-13]。衛(wèi)星平臺(tái)抖動(dòng)功率譜密度匯集在圖3中，包括美國(guó)“陸地衛(wèi)星四號(hào)”（Landsat-4）測(cè)量值、歐空局“奧林普斯”衛(wèi)星（OLYMPUS）測(cè)量值、日本“工程試驗(yàn)”衛(wèi)星ETS系列（ETS-VI）實(shí)測(cè)值、國(guó)際空間站的有限元分析值等[14]。

圖2 利用像移估計(jì)實(shí)現(xiàn)遙感圖像清晰化原理框圖Fig.2 The schematic of image motion estimation to make the remote sensing image clear

圖3 各種平臺(tái)的功率譜密度曲線Fig.3 The spectral density at different platforms

經(jīng)功率譜對(duì)比分析可知，衛(wèi)星振動(dòng)頻譜較寬，最寬可以達(dá)到1kHz，但是一般在80Hz處功率譜會(huì)有明顯下降，如圖3紅色框所示。因此，對(duì)成像載荷影響最大的振動(dòng)主要集中在80Hz以內(nèi)。要實(shí)現(xiàn)高精度的運(yùn)動(dòng)估計(jì)，根據(jù)工程經(jīng)驗(yàn)，運(yùn)動(dòng)測(cè)量的頻率不能夠低于運(yùn)動(dòng)頻率的5倍，因此對(duì)80Hz的抖動(dòng)進(jìn)行運(yùn)動(dòng)估計(jì)，最低需要為400幀/s拍攝速度。高速CCD相機(jī)成像尺寸小，曝光時(shí)間短，可以滿足400幀/s拍攝速度要求，所以可以詳細(xì)記錄主成像相機(jī)成像期間衛(wèi)星的抖動(dòng)。抖動(dòng)的軌跡既是連續(xù)幀圖像間像移估計(jì)結(jié)果隨時(shí)間的曲線[15]。最后，利用像移估計(jì)結(jié)果反作用于遙感圖像，即可將其清晰化。

由于國(guó)內(nèi)外現(xiàn)有的星載光學(xué)成像遙感器沒有 CCD相機(jī)與主成像相機(jī)同軸同步成像的應(yīng)用實(shí)現(xiàn)，故而無法獲得實(shí)際的連續(xù)幀高速 CCD相機(jī)成像圖片。本文通過對(duì)已有的高分辨率遙感圖像進(jìn)行裁剪，從而獲得等效的高速 CCD相機(jī)成像圖片。通過裁剪同一區(qū)域遙感圖像不同的偏移來模擬因衛(wèi)星抖動(dòng)而造成的像移偏差，最后利用相位相關(guān)算法對(duì)這些像移偏差進(jìn)行估計(jì)，并評(píng)價(jià)估計(jì)結(jié)果準(zhǔn)確度。

1 相位相關(guān)法像移估計(jì)原理

相位相關(guān)的實(shí)現(xiàn)原理[16]是基于傅里葉變換的平移性質(zhì)。假設(shè)一幅圖像f1(x,y)平移后形成另一幅圖像f2(x,y)，相對(duì)平移量為 (x0,y0)像素，則兩幅圖像分別進(jìn)行傅里葉變換后的關(guān)系式如式（1）所示

式中F1(u,v)是圖像f1(x,y)的傅里葉變換；F2(u,v)是圖像f2(x,y)的傅里葉變換。

式（1）說明兩幅圖像變換到頻域中有相同的幅值，而且相對(duì)關(guān)系僅與相對(duì)平移量有關(guān)。兩個(gè)圖像的相位差等于圖像歸一化互功率譜的相位，如式（2）

式中為F2(u,v)的共軛。

為求出相對(duì)平移量 (x0,y0)的具體值，只需在對(duì)式（2）進(jìn)行傅里葉反變換即可，如式（3）

式中 IFFT表示傅里葉反變換；δ表示沖擊函數(shù)，沖激函數(shù)只有在零點(diǎn)才有值，其他點(diǎn)均為零。即當(dāng)(x,y)的值為 (x0,y0)時(shí)該函數(shù)才有值，也是沖擊函數(shù)的最大值。

由式（3）可以看出，互功率譜的傅里葉反變換恰好為平移點(diǎn) (x0,y0)的沖擊函數(shù)，由此可以得到圖片運(yùn)動(dòng)的估計(jì)值。

2 相位相關(guān)法測(cè)試精度分析

以局部城市遙感圖片為例，設(shè)定偏移大小為（5，4.6）像素。截取原始實(shí)際地面城市區(qū)域局部圖片（256像素×256像素）、CCD成像圖片、CCD成像偏移圖片分別如圖4所示。圖4中（a）表示原始圖像，（b）表示（a）分離出的圖像，（c）表示（b）平移后的圖像，（d）是（b）與（c）做相位相關(guān)算法后的結(jié)果。（d）圖Z軸為沖擊函數(shù)值（無單位）。

為了測(cè)試相位相關(guān)法的穩(wěn)像精度，設(shè)置如下測(cè)試條件。

1）分別選取大小為5 000像素×5 000像素，分辨率為0.5m的城市區(qū)域圖片做為實(shí)際地面參考圖像。

2）相位相關(guān)算法理論能夠準(zhǔn)確識(shí)別的相對(duì)運(yùn)動(dòng)分辨率為 0.1個(gè)像素[17]。為了使CCD圖片擁有0.1像素分辨率的圖像信息，截取地面參考圖片的大小為2 560像素×2 560像素，在進(jìn)行10倍的縮小處理，處理方式是每10個(gè)×10個(gè)參考圖片的像素值的和取平均作為CCD圖片的像素點(diǎn)值[18]。

3）由于相位相關(guān)算法求出兩個(gè)圖像的相對(duì)平移量為 (x0,y0)是整數(shù)像素值[19]。在圖片平移并非整數(shù)像素時(shí)，雖然 (x0,y0)坐標(biāo)處的值仍然為沖擊函數(shù)的最大值，但已經(jīng)不是唯一有值的點(diǎn)，它周圍的坐標(biāo)也有值，并且形成比較突出的圓錐形，如圖4（d）。這表明兩個(gè)圖片的互功率譜能量主要集中在這個(gè)錐形上面，只要截取代表這個(gè)錐形的坐標(biāo)值并做平均加權(quán)處理，即可得到一個(gè)非整數(shù)的相對(duì)平移量。所以這里截取以 (x0,y0)為中心的相鄰區(qū)域的一些像素值做平均處理，分別采取 5像素×5像素的矩形圖像灰度值擬合，3像素×3像素的矩形圖像灰度值擬合，質(zhì)心擬合，過 (x0,y0)的3點(diǎn)連線像素灰度值擬合方式[20]，求得實(shí)際圖片的亞像素級(jí)偏移量，并評(píng)估不同擬合方式的估計(jì)精度。

4）因CCD成像速度與圖片大小為反比，為達(dá)到400幀/s的成像速度，用于穩(wěn)像的CCD圖片的大小最大為256像素×256像素[21]，同時(shí)測(cè)試圖片大小為128像素×128像素、64像素×64像素、32像素×32像素。

5）CCD成像相機(jī)X軸方向和Y軸方向的平移范圍均為-5～+5個(gè)像素，即需要實(shí)際地面參考圖片平移-50～+50個(gè)像素。

圖4 使用相位相關(guān)法進(jìn)行像移估計(jì)流程舉例Fig.4 An example using the phase correlation method in image motion estimation

由圖4可以看出，最高脈沖峰的坐標(biāo)在（50, 90），求得5像素×5像素圖像灰度值擬合后得到的平移坐標(biāo)為（55.0024, 94.6476）。因?yàn)榇怂惴ㄓ肕atlab軟件實(shí)現(xiàn)，已經(jīng)將零頻率得到的坐標(biāo)移到圖4中心，需要減去平移量，得到的最終平移坐標(biāo)為（5.0024，4.6476）像素，相對(duì)誤差為（0.002，0.0476）個(gè)像素。

由于CCD相機(jī)圖片大小的限制，僅截取了參考遙感圖片的一小部分作為相位相關(guān)算法的原始圖像。為了方便得到不同的 CCD圖像大小、不同的圖像偏移量、不同的圖像灰度值擬合方式對(duì)算法精度的影響，使用Matlab制作了測(cè)試軟件，如圖5所示。

圖5 像移估計(jì)算法精度分析軟件Fig.5 The precision analysis software of the image motion estimation algorithm

當(dāng)載入?yún)⒖歼b感圖片后，使用測(cè)試軟件截取參考圖片的某一部分，指定該部分的大小與需要偏移的像素值，生成原始CCD圖片與偏移后的CCD圖像，就可以進(jìn)行像移估計(jì)了。

首先使用測(cè)試軟件測(cè)試不同沖擊函數(shù)擬合方式對(duì)像移估計(jì)的影響，將城市區(qū)域圖片進(jìn)行分塊處理，按2 560像素×2 560像素分割成256個(gè)不同區(qū)域的圖片，并將這些圖片分別作為原始圖片。按上述5條測(cè)試條件進(jìn)行相位相關(guān)算法的像移估計(jì)的精度分析。將分析結(jié)果進(jìn)行統(tǒng)計(jì)，如表1所示。

表1 不同擬合算法精度分布Tab.1 The accuracy distribution of different fitting algorithms 單位：像素

原始圖像與平移后圖像，在水平與垂直方向的平移設(shè)定為最大±5個(gè)像素的前提下，水平與垂直方向的平移大小以 0.1像素分辨率任意組合的情況下，以圖像明暗的方式分別列出每種組合對(duì)應(yīng)的估計(jì)誤差大小。以表1中第一個(gè)小圖為例，圖中任意點(diǎn)表示的是該點(diǎn)對(duì)應(yīng)的水平和垂直坐標(biāo)值下，水平方向的誤差估計(jì)值。誤差越小，圖像點(diǎn)越暗；誤差越大，圖像點(diǎn)越亮。對(duì)誤差值大小進(jìn)行歸納總結(jié)，如表2所示。

表2 相位相關(guān)算法誤差綜合分析Tab.2 Comprehensive error analysis of the phase correlation algorithm 單位：像素

由表2可以看出，最小誤差近似為0，出現(xiàn)在CCD連續(xù)幀圖像相對(duì)平移量為（0，0）坐標(biāo)。由表1得出，最大誤差均出現(xiàn)在0.5個(gè)像素整數(shù)倍平移位置，具體坐標(biāo)與圖像本身有關(guān)。

表3 最優(yōu)運(yùn)動(dòng)誤差估計(jì)Tab.3 Maximum motion error estimation 單位：像素

表3進(jìn)一步對(duì)表2的精度分布進(jìn)行了總結(jié)。得出使用相位相關(guān)算法，結(jié)合質(zhì)心擬合計(jì)算做亞像素的精確擬合，能獲得最優(yōu)像移估計(jì)結(jié)果。

在固定的成像幀頻，圖像越小，計(jì)算量越小，所以如果能夠用更小的圖片達(dá)到相同的效果，則將大大地減少硬件壓力。

為了充分驗(yàn)證相關(guān)算法對(duì)不同圖像大小進(jìn)行像移估計(jì)的準(zhǔn)確度，將城市區(qū)域圖片進(jìn)行分塊處理，分別按320像素×320像素、640像素×640像素、1 280像素×1 280像素和2 560像素×2 560像素分割成256個(gè)不同區(qū)域的圖片，并將這些圖片分別作為原始圖像，進(jìn)行像移估計(jì)的精度分析。使用質(zhì)心擬合方式作為像移估計(jì)擬合方式，按成像大小將分析結(jié)果進(jìn)行統(tǒng)計(jì)，如表4～6所示。

表4 CCD成像大小64像素×64像素的相位相關(guān)算法精度分析Tab.4 Accuracy analysis of the phase correlation algorithm for a CCD image with 64 pixels×64 pixels 單位：像素

表5 CCD成像大小128像素×128像素的相位相關(guān)算法精度分析Tab.5 Accuracy analysis of the phase correlation algorithm for a CCD image with 128 pixels×128 pixels 單位：像素

表6 CCD成像大小256像素×256像素的相位相關(guān)算法精度分析Tab.6 Accuracy analysis of the phase correlation algorithm for a CCD image with 256 pixels×256 pixels 單位：像素

由表4、5、6可以看出，最小誤差均為0，其出現(xiàn)點(diǎn)為成像圖片和偏移圖片重合點(diǎn)位置。最大誤差均出現(xiàn)在 0.5個(gè)像素整數(shù)倍平移位置，具體坐標(biāo)與圖像本身有關(guān)。成像算法的最大誤差值與成像圖像大小成反比，大約按 1.5倍關(guān)系遞減，這說明如果需要較高的像移估計(jì)精度，則需采用更大的成像圖像，但耗時(shí)更多，需要折中處理。

在使用CCD成像32像素×32像素圖片進(jìn)行像移估計(jì)時(shí)，有可能導(dǎo)致相關(guān)算法失效，即當(dāng)平移圖像平移到某個(gè)點(diǎn)或某片區(qū)域時(shí)，利用相關(guān)算法計(jì)算出的誤差已經(jīng)為多個(gè)像素。具體失效區(qū)域與圖片的具體內(nèi)容相關(guān)。對(duì)所有的分割圖片進(jìn)行統(tǒng)計(jì)，假定像素偏移誤差大于0.5個(gè)像素則為失效，失效概率為92%；若假定像素偏移誤差大于1個(gè)像素則為失效，失效概率為86%。這說明針對(duì)32×32像素圖片進(jìn)行像移估計(jì)，失效概率過高，該方法無法使用。失效形式如圖6所示兩種情況。

圖6 32像素×32像素CCD成像圖片像移估計(jì)算法失效舉例Fig.6 A failure example of the image motion estimation algorithm for a CCD image with 32 pixels×32 pixels

圖6與表1表達(dá)的物理意義是一樣的，不同點(diǎn)在于表 1每個(gè)小圖的值是隨X和Y的值連續(xù)變化，而圖6的值出現(xiàn)了突變，且突變點(diǎn)是沒有規(guī)律的，分別為圖6（a）的大面積突變失效和圖6（b）的某幾個(gè)點(diǎn)失效，如果實(shí)際像移點(diǎn)出現(xiàn)在失效點(diǎn)，像移估計(jì)將出現(xiàn)極大誤差，甚至超過了真實(shí)像移值本身。這說明當(dāng)CCD成像大小為32像素×32像素時(shí)，由于圖像信息采集量不足，在某些情況下使用相位相關(guān)算法采集到的互功率譜的傅里葉反變換峰值不再特別獨(dú)立[22]。互功率譜的傅里葉反變換圖像如圖7所示。

圖7 32像素×32像素CCD成像圖像的互功率譜的傅里葉反變換圖片F(xiàn)ig.7 The inverse fourier transform image of cross power spectrum of the CCD image with 32 pixels×32 pixels

由圖7可以看出，極大值不唯一，雖然最大值唯一，但算法可能無法提取圖像的最大值，而錯(cuò)把臨近的極大值作為最大值對(duì)應(yīng)的坐標(biāo)輸出。這使得相位相關(guān)算法失效，獲得的值也是錯(cuò)值。

3 結(jié)束語

本文主要分析了對(duì)高速 CCD相機(jī)對(duì)地成像的圖片使用相位相關(guān)算法提取像移信息，復(fù)原遙感圖像所能達(dá)到的精度與算法使用條件。分析結(jié)果表明，在滿足 0.1像素分辨率的前提下，更大的圖像、整數(shù)倍像素偏移、使用質(zhì)心擬合方式可獲得更好的估計(jì)精度。當(dāng)圖片小到32像素×32像素時(shí)，有可能引起該方法失效。

雖然是使用仿真的方式對(duì)相位相關(guān)算法的像移估計(jì)精度進(jìn)行的分析，但使用的遙感圖像等數(shù)據(jù)源均來自于實(shí)際遙感成像圖像，因此仿真算法也可直接使用于真實(shí)系統(tǒng)中。本研究也為未來實(shí)際的在軌應(yīng)用提供了實(shí)驗(yàn)依據(jù)。

本研究沒有考慮到實(shí)際測(cè)量噪聲對(duì)相位相關(guān)算法的影響[23]，后續(xù)將開展不同信噪比對(duì)算法估計(jì)精度影響的分析。