基于參叉像元和非均勻B樣條曲面的遙感圖像超分辨率重建

王京萌， 張愛武， 孟憲剛， 劉詔

(1.首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京 100048；2.三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048；3.資源環(huán)境與地理信息系統(tǒng)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048；4.城市環(huán)境過程與數(shù)字模擬國家重點(diǎn)實(shí)驗(yàn)室培育基地，北京 100048)

基于參叉像元和非均勻B樣條曲面的遙感圖像超分辨率重建

王京萌1，2，3，4， 張愛武1，2，3，4， 孟憲剛1，2，3，4， 劉詔1，2，3，4

(1.首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京 100048；2.三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048；3.資源環(huán)境與地理信息系統(tǒng)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048；4.城市環(huán)境過程與數(shù)字模擬國家重點(diǎn)實(shí)驗(yàn)室培育基地，北京 100048)

遙感圖像之間的信息互補(bǔ)可以提高圖像分辨率，但插值方法易使圖像邊界模糊、部分細(xì)節(jié)信息丟失。針對(duì)這一問題，提出一種基于參叉像元與非均勻B樣條插值相結(jié)合的遙感圖像超分辨率重建方法。利用經(jīng)灰度匹配和亞像元級(jí)幾何配準(zhǔn)的2景低分辨率圖像，通過參叉交錯(cuò)像元采樣到原圖像網(wǎng)格2倍的網(wǎng)格中；對(duì)于沒有值的坐標(biāo)處用三次B樣條插值，插值時(shí)選用非均勻的節(jié)點(diǎn)參數(shù)化方法，曲面圖像網(wǎng)格點(diǎn)由鄰域36個(gè)已知像元組成；在求解待插值點(diǎn)參數(shù)值時(shí)引入平行線法和黃金分割法迭代尋找最優(yōu)值，使插值更準(zhǔn)確；最后對(duì)插值后的圖像進(jìn)行復(fù)原處理，重建可視效果更好的“高”分辨率圖像。對(duì)實(shí)驗(yàn)圖像的評(píng)價(jià)表明，用本文方法重建的圖像在清晰度、信息量、信噪比和分辨率等方面均有較大的提高。

圖像超分辨率重建；非均勻B樣條曲面；參叉像元

0 引言

近年來，采用亞像元參叉線陣提高圖像分辨率的成像技術(shù)已在很多領(lǐng)域應(yīng)用，并收到了很好的效果。德國宇航中心研制的熱點(diǎn)識(shí)別傳感器(hot spot recognition sensors，HSRS)以及Lange等[1]提出的新型CCD設(shè)計(jì)方案中，在焦平面上放置2片線陣方向上錯(cuò)開半個(gè)像元的CCD，將微掃描轉(zhuǎn)化為靜態(tài)錯(cuò)位，從而避免了移位精度的難題。此外，航空攝影系統(tǒng)ADS40的全色波段也采用了和SPOT5類似的設(shè)計(jì)，以實(shí)現(xiàn)超分辨率[2]。繼Andrews等[3]首先提出將三次樣條插值應(yīng)用于圖像的重采樣之后，Panda等[4]提出了廣義B樣條(Bezier樣條是樣條曲線的一種特殊表示形式，它是B樣條基曲線的線性組合)插值可用于圖像放大。Unser等[5]將數(shù)字圖像看成是一個(gè)以像元點(diǎn)組成的網(wǎng)格控制點(diǎn)的曲面，構(gòu)造B樣條曲面(用多個(gè) B 樣條曲線畫出的曲面)進(jìn)行插值。馮杰飛[6]采用基于圖像亮度變化的準(zhǔn)均勻參數(shù)化方法來構(gòu)造B樣條曲面。安博文等[7]則對(duì)圖像進(jìn)行亞像元級(jí)幾何配準(zhǔn)后，利用待插值點(diǎn)周圍的4條B樣條曲線實(shí)現(xiàn)圖像插值融合。在圖像應(yīng)用中,均勻或準(zhǔn)均勻B樣條參數(shù)化的方法居多，由于圖像像元的亮度值不一定是等距分布的，故節(jié)點(diǎn)參數(shù)化時(shí)選擇均勻和準(zhǔn)均勻的方法欠妥；B樣條曲線的方法只考慮了線性一維的因素，沒有考慮周圍各像元對(duì)中心像元的影響；插值或者融合的方法給圖像帶來的模糊與噪聲很大，這些算法并沒有考慮插值后的復(fù)原問題。

本文參考亞像元參叉線陣的思想，對(duì)同一研究區(qū)經(jīng)過灰度匹配和亞像元級(jí)配準(zhǔn)的2景圖像，參叉采樣到原圖像網(wǎng)格的2倍的網(wǎng)格中，使用非均勻B樣條插值。由于控制點(diǎn)之間并不一定呈等距分布，因此本文節(jié)點(diǎn)參數(shù)化選擇非均勻法。將待插值點(diǎn)鄰域36個(gè)已知像元張成插值曲面，反算控制頂點(diǎn)時(shí)引入“非節(jié)點(diǎn)”的方法；計(jì)算待插值點(diǎn)參數(shù)值時(shí)利用平行線法和黃金分割法，迭代尋找到最優(yōu)的參數(shù)值。最終經(jīng)過圖像復(fù)原得到1景高分辨率的重建圖像。

1 研究方法

1.1 參叉線陣思想

保持CCD探測(cè)器的像元尺寸和相機(jī)的光學(xué)系統(tǒng)不變，但在相機(jī)焦面上除原來的CCD線陣外，增加一條相同的CCD線陣。2線陣間在X方向錯(cuò)開1/2像元，Y方向錯(cuò)開N+1/2像元(N為盡可能小的整數(shù))[8]是亞像元?jiǎng)討B(tài)成像的基本思想。本文參考這種方法，對(duì)線陣CCD在不同時(shí)刻、不同角度分別獲取的同一區(qū)域的2景圖像，經(jīng)灰度匹配和亞像元級(jí)精確幾何配準(zhǔn)后，分別填入到原始格網(wǎng)2倍采樣密度的新格網(wǎng)中；對(duì)沒有值的坐標(biāo)處用零補(bǔ)足(圖1)，插值時(shí)要在零值的位置處插值；最后進(jìn)行圖像復(fù)原，得到超分辨率圖像。

圖1 亞像元參叉數(shù)組采樣網(wǎng)格Fig.1 Sampling grid by sub-pixel staggered array

盡管在水平、垂直方向采樣間隔分別增加了1倍，但根據(jù)參叉像元的B樣條超分辨率復(fù)原理論，提高圖像分辨率倍數(shù)為接近1.6倍，而非2倍[9]，在本質(zhì)上不是真正將分辨率增加1倍，而是通過插值方法獲得2倍采樣網(wǎng)格中沒有的像元值。

1.2 三次B樣條插值

1.2.1 B樣條曲線的節(jié)點(diǎn)參數(shù)化

B樣條方法具有表示和設(shè)計(jì)自由型曲線、曲面的強(qiáng)大功能，是最廣為流行的形狀數(shù)學(xué)描述的主流方法之一[10-12]。B樣條的方法應(yīng)用到圖像插值中，不同的節(jié)點(diǎn)參數(shù)化方法對(duì)圖像處理效果有差異，本文使用非均勻的節(jié)點(diǎn)參數(shù)化方法。

令U={u0,u1,…,um}是一個(gè)單調(diào)不減的實(shí)數(shù)序列，即ui≤ui+1，i=0,±1,…,m-1。其中ui稱為節(jié)點(diǎn)，U稱為節(jié)點(diǎn)矢量，區(qū)間[ui,ui+1]為第i個(gè)節(jié)點(diǎn)區(qū)間。用Ni,p(u)表示第i個(gè)p次(p+1階)B樣條基函數(shù)，即

(1)

(2)

同時(shí)約定0/0=0，B樣條曲面是由2個(gè)方向的控制點(diǎn)網(wǎng)絡(luò)、2個(gè)節(jié)點(diǎn)矢量和單變量B樣條基函數(shù)的乘積來定義，即

(3)

式中：Pi,j為空間中給定的(n+1)(m+1)個(gè)網(wǎng)點(diǎn)，通常稱為S(u,v)的控制頂點(diǎn)；Ni,p(u)和Nj,q(v)分別為關(guān)于節(jié)點(diǎn)向量U和V的p+1階的B樣條基函數(shù)。

1.2.2 圖像的B樣條插值

以往的基于B樣條的圖像插值算法認(rèn)為圖像像元點(diǎn)在二維空間呈等距分布，故采用均勻參數(shù)化方法。若將圖像置于三維空間中，圖像空間中像元點(diǎn)的位置坐標(biāo)作為三維空間中的x，y軸坐標(biāo)值，圖像像元值作為p軸坐標(biāo)值[6]，則各控制點(diǎn)之間并不一定呈等距分布。因此，本文對(duì)節(jié)點(diǎn)參數(shù)化選擇非均勻法。綜合考慮插值的效果和效率，選擇了三次B樣條曲面插值。在進(jìn)行插值之前，使用反射填充法補(bǔ)齊邊界點(diǎn)。

1.2.2.1 反算控制頂點(diǎn)

通常的算法是將數(shù)據(jù)點(diǎn)Q0與Qn分別作為三次B樣條插值曲線的首、末端點(diǎn)，將內(nèi)部點(diǎn)Q1,Q2,…,Qn-1依次作為分段連接點(diǎn)，則將曲線分為n段，而所求的控制頂點(diǎn)數(shù)為n+3個(gè)，方程為n+1個(gè)，需要補(bǔ)充2個(gè)邊界條件確定的附加方程才能求出全部控制頂點(diǎn)。依據(jù)文獻(xiàn)[13]可以采用“非節(jié)點(diǎn)”條件，引入t1與t2這2個(gè)參數(shù)，使得待求的控制點(diǎn)減少為n+1個(gè)，從而反算出控制頂點(diǎn)。

給定(n+1)(m+1)個(gè)數(shù)據(jù)點(diǎn)，創(chuàng)建一個(gè)非有理的(p，q)次B樣條插值曲面使其通過這些點(diǎn)，即

(4)

(5)

1.2.2.2 計(jì)算待插值點(diǎn)參數(shù)值

要迅速找到二元函數(shù)z=f(x,y)的最小值及其對(duì)應(yīng)的(x，y)點(diǎn)。在2個(gè)因素中如果一個(gè)因素易于調(diào)整，另一個(gè)不易調(diào)整，則使用“平行線法”，先將y固定在范圍(c，d)的0.618處，即取

y=c+ω(d-c),

(6)

式中ω為系數(shù)，ω=0.618 033 988。用單因素法(黃金分割法)找最小值，假定在P點(diǎn)取得這一值，再把y固定在范圍(c，d)的0.382處，即取

y=c+(1-ω)(d-c)。

(7)

假定在Q點(diǎn)取得最小值，則比較P和Q的結(jié)果。如果P為所求，則去掉Q點(diǎn)下面的部分，剩余的部分同上。黃金分割法(0.618法)[14]的具體算法如下：

設(shè)φ(k):[a,b]→R，若存在點(diǎn)k*(a,b)，使得φ(k)在[a,k*]上單調(diào)減小、在[k*,b]上單調(diào)增加，則稱[a,b]為φ(k)的單峰區(qū)間，φ(k)為[a.b]上的單峰函數(shù)。

1)計(jì)算k1=a+(1-ω)(b-a),k2=a+ω(b-a),φ(k1)=φ1,φ(k2)=φ2；

3)若φ1<φ2，則b←k2,k2←k1,φ2←φ1，轉(zhuǎn)向4)；若φ1>φ2，則a←k1,k1←k2，φ1←φ2，轉(zhuǎn)向5)；若φ1=φ2，則a←k1,b←k2，轉(zhuǎn)向1)；

4)計(jì)算k1=a+(1-ω)(b-a)，φ1←φ(k1)，轉(zhuǎn)向2)；

5)計(jì)算k2=a+ω(b-a)，φ2←φ(k2)，轉(zhuǎn)向2)。

將每個(gè)控制點(diǎn)的屬性值p值獨(dú)立出來構(gòu)成矩陣Pz，然后將待插值點(diǎn)的參數(shù)值帶入曲面公式，從而得到待插值點(diǎn)的像元值。

1.3 實(shí)驗(yàn)處理流程

本文建立了基于參叉像元和非均勻B樣條曲面的遙感圖像超分辨率重建技術(shù)流程(圖2)。

圖2 遙感圖像超分辨率重建技術(shù)流程Fig.2 Technical process of super-resolution reconstruction for remote sensing image

對(duì)技術(shù)流程的詳細(xì)說明如下：

1)輸入： 2景低分辨率圖像A和B。

①預(yù)處理。先對(duì)獲取的2景圖像A和B分別進(jìn)行輻射校正和幾何糾正，得到圖像A1和B1；再對(duì)A1和B1進(jìn)行像元級(jí)幾何配準(zhǔn)，得到圖像A1和B2；②對(duì)圖像A1和B2分別進(jìn)行灰度匹配，得到圖像A2和B3； ③對(duì)灰度匹配后的A2和B3進(jìn)行亞像元幾何配準(zhǔn)，得到圖像A2和B4；④圖像重建。先將圖像A2和B4用參叉像元采樣到1景是原來分辨率2倍的圖像網(wǎng)格P中，無像元的位置用零補(bǔ)充；再對(duì)零值位置進(jìn)行非均勻B樣條插值，得到圖像P1；⑤圖像復(fù)原。對(duì)插值后圖像P1進(jìn)行三階波特沃斯帶阻濾波器去除噪聲復(fù)原，得到圖像P2。

2)輸出： 1景高分辨率圖像P2。

為了說明本文方法的效果，對(duì)原低分辨率圖像進(jìn)行了雙三次插值和1組均勻B樣條插值，將這2組插值處理結(jié)果與本文算法進(jìn)行了對(duì)比實(shí)驗(yàn)。

2 遙感圖像超分重建實(shí)驗(yàn)

實(shí)驗(yàn)數(shù)據(jù)為北京某濕地線陣CCD相機(jī)航空攝影數(shù)據(jù)，飛行高度約300 m，圖像大小為800像元×800像元，使用灰度圖像進(jìn)行處理。原始圖像數(shù)據(jù)采用CaliGeo軟件進(jìn)行輻射校正和幾何初步糾正。對(duì)預(yù)處理后的2景待超分重建的圖像先進(jìn)行像元級(jí)的幾何配準(zhǔn)。

2.1 灰度匹配

原始圖像經(jīng)過預(yù)處理和初步幾何配準(zhǔn)后，由于待復(fù)原的2景圖像是在不同時(shí)刻獲取的，在獲取過程中因受大氣、光照等多種因素的影響，2景圖像的灰度值會(huì)出現(xiàn)差異；所以必須要用直方圖匹配的方法對(duì)待復(fù)原圖像進(jìn)行灰度匹配，即以其中1景圖像的直方圖作為基準(zhǔn)直方圖，對(duì)另1景圖像進(jìn)行直方圖匹配。假設(shè)r被歸一到區(qū)間[0,1]，1景圖像的灰度級(jí)可被視為該區(qū)間的隨機(jī)變量，則對(duì)于任意r存在變換

s=T(r)，

(8)

式中：s為對(duì)于每個(gè)像元值r產(chǎn)生的1個(gè)對(duì)應(yīng)灰度值，0≤r≤1；T(r)為變換函數(shù)。

令Pr(r)為隨機(jī)變量r的概率密度函數(shù)。對(duì)于1景圖像，L為可能的灰度級(jí)總數(shù)，則有

(9)

令pz(z)為基準(zhǔn)圖像的直方圖，pr(r)為待匹配圖像直方圖，n為圖像中像元總數(shù)，nj為灰度級(jí)為rj的像元數(shù)，設(shè)存在變換函數(shù)G(s)，則變量vk和zk分別為

(10)

zk=G-1[T(rk)]=G-1[T(sk)]。

(11)

式(9)-(11)中：k為元素在陣列中的位置，k=0,1,2,…,L-1；s為該位置的值；L=256，且為整數(shù)。初步幾何配準(zhǔn)前后圖像如圖3所示; 灰度匹配前后的直方圖如圖4所示。

(a) 低分辨率圖像1 (b) 低分辨率圖像2(c) 以圖像1為基準(zhǔn)配準(zhǔn)的圖像2

圖3 初步幾何配準(zhǔn)前后圖像

Fig.3 Images before and after preliminary geometric registration

(a) 低分辨率圖像1 (b) 低分辨率圖像2 (c) 灰度匹配后

圖4 灰度匹配前后直方圖

Fig.4 Histograms before and after gray level matching

2.2 亞像元級(jí)配準(zhǔn)

當(dāng)待超分重建圖像的位移量為整數(shù)像元時(shí)，超分辨率的效果并不明顯；但存在亞像元位移時(shí)，通過信息融合即可獲得更多信息，得到高質(zhì)量圖像。

由于在圖像預(yù)處理中進(jìn)行了像元級(jí)的初步幾何配準(zhǔn)，使縮放尺度和幾何變化基本確定，故亞像元的配準(zhǔn)主要估算平移和旋轉(zhuǎn)參數(shù)。利用傅立葉域的平移不變性可以很好地估計(jì)平移參數(shù)和旋轉(zhuǎn)參數(shù)，對(duì)旋轉(zhuǎn)角度θ的估算使用Vanewalle[15]等提出的頻域配準(zhǔn)算法。設(shè)參考圖像和待配準(zhǔn)圖像分別為s(x)和f(x)，x=(x1,x2)T，則

s(x)=f[R(x+△x)],

(12)

(13)

求式(13)中S(u)的模，可得

|S(u)|=|F(R(u))|。

(14)

經(jīng)過用式(14)對(duì)幾何配準(zhǔn)后的圖像模型求模說明，圖像旋轉(zhuǎn)后的振幅只與旋轉(zhuǎn)角度有關(guān)，與位移無關(guān)，這樣就可以先求出旋轉(zhuǎn)角度。經(jīng)過計(jì)算，本實(shí)驗(yàn)中圖像2相對(duì)于圖像1的旋轉(zhuǎn)角度為0.02°。

亞像元的位移量估計(jì)使用多峰擬合的相位相關(guān)配準(zhǔn)方法。文獻(xiàn)[16]指出，當(dāng)圖像出現(xiàn)亞像元位移時(shí)，位移量落在相位相關(guān)頻譜的最高峰與其他幾個(gè)峰值之間。設(shè)最高峰為M1，其他峰值為M2，M3，M4，對(duì)應(yīng)位置坐標(biāo)分別為(x1,y1)，(x2,y2),(x3,y3),(x4,y4),實(shí)際位移量為(x0,y0)[7]，使用雙線性擬合即可求出圖像2相對(duì)于圖像1的亞像元位移量為(0.27,-0.13)，糾正后的亞像元配準(zhǔn)結(jié)果為(x+0.27,y-0.13)(右、上為正，左、下為負(fù))。

對(duì)于精度的測(cè)定，選取圖像1作為參考圖像，對(duì)其先進(jìn)行垂直方向和水平方向的亞像元精度的位移，再進(jìn)行小角度旋轉(zhuǎn)，得到待配準(zhǔn)的圖像；使用本文方法對(duì)待配準(zhǔn)圖像進(jìn)行亞像元的參數(shù)估計(jì)，將其與已知參數(shù)比較，用均方誤差來衡量配準(zhǔn)精度，即

(15)

式中：p(i)為已知配準(zhǔn)參數(shù)；p′(i)為估算出的配準(zhǔn)參數(shù)；n為圖像的像元總數(shù)。

除使用本文圖像進(jìn)行配準(zhǔn)精度計(jì)算外，還選取了另一組圖像作為衡量配準(zhǔn)精度的參考，以說明此配準(zhǔn)精度不僅限于本文的圖像，其他圖像配準(zhǔn)精度也都與本文圖像配準(zhǔn)精度相近。本文亞像元配準(zhǔn)算法的精度如表1，平均MSE=0.281 8<1。

表1 配準(zhǔn)參數(shù)估計(jì)值及配準(zhǔn)精度Tab.1 Estimated values of registration coefficient and registration accuracy

亞像元配準(zhǔn)的精度直接影響圖像復(fù)原的效果，應(yīng)該使其精度的MSE控制在小于1的范圍內(nèi)。本文亞像元配準(zhǔn)精度的MSE大約在0.281 8左右。即可使每個(gè)像元的信息較原單個(gè)像元信息增加，從而提高單位像元的信息量，達(dá)到超分辨率的效果。

2.3 非均勻B樣條曲面圖像插值

本文選用6像元×6像元的數(shù)據(jù)點(diǎn)，p=q=7(p，q分別為水平和垂直方向B樣條基函數(shù)的次數(shù))，這36個(gè)點(diǎn)為距離待插值點(diǎn)最近的點(diǎn)，形狀如菱形(圖5)。

圖5 曲面B樣條插值鄰域示意圖Fig.5 Sketch map of curved surface B-spline interpolation neighborhood

為了便于計(jì)算，將標(biāo)號(hào)為1和6的交點(diǎn)作為原點(diǎn)，這樣可以使得36個(gè)已知數(shù)據(jù)點(diǎn)的坐標(biāo)均為非負(fù)。按照標(biāo)號(hào)順序，36個(gè)數(shù)據(jù)點(diǎn)的坐標(biāo)依次為： (5，0)，(4，1)，(3，2)，(2，3)，(1，4)，(0，5)；(6，1)，(5，2)，(4，3)，(3，4)，(2，5)，(1，6)；(7，2)，(6，3)，(5，4)，(4，5)，(3，6)，(2，7)；(8，3)，(7，4)，(6，5)，(5，6)，(4，7)，(3，8)；(9，4)，(8，5)，(7，6)，(6，7)，(5，8)，(4，9)；(10，5)，(9，6)，(8，7)，(7，8)，(6，9)，(5，10)。待插值點(diǎn)的坐標(biāo)為(5，5)。此點(diǎn)在u和v軸上的參數(shù)值的估計(jì)采用中雙因素優(yōu)選法中的平行線法。

使用本文1.2.1中的式(3)完成參數(shù)化。為了避免另外給出邊界條件，本文采用文獻(xiàn)[8]提出的“非節(jié)點(diǎn)”條件，這樣可以使得控制頂點(diǎn)減少到n+1個(gè)；對(duì)于構(gòu)造三次B樣條插值曲線而言，假設(shè)給定n+1個(gè)數(shù)據(jù)點(diǎn)Qi,其中i=0,1,…,n。具體做法為： 首、末數(shù)據(jù)點(diǎn)Q0與Qn依然為曲線的首、末端點(diǎn)，把Q2,…,Qn-2依次作為曲線的分段連接點(diǎn)，而Q1與Qn-1分別是首段上和末段上對(duì)應(yīng)于參數(shù)t1，t2的點(diǎn)；t1，t2滿足t1[u3,u4]，t2[un,un+1]，即Q1與Qn-1不再作為分段連接點(diǎn)。

2.4 圖像復(fù)原

插值方法都會(huì)引入噪聲，故本文復(fù)原的方法僅考慮噪聲存在情況下的復(fù)原，插值后圖像表現(xiàn)為類似沖擊的串，這種串在傅立葉頻譜中的表現(xiàn)如圖6所示。在高頻部分有4個(gè)明顯的亮點(diǎn)，經(jīng)分析發(fā)現(xiàn)這種噪聲為周期噪聲，通過陷波濾波可濾除這些成分。本文使用三階巴特沃斯帶阻濾波器進(jìn)行噪聲去除[17]。

圖6 非均勻B樣條曲面插值后頻域圖像(黑色圓圈為帶阻濾波器形狀)Fig.6 Image in frequency domain after non-uniform B-spline curved surface interpolation

對(duì)不同噪聲去除方法的試驗(yàn)對(duì)比表明，本文復(fù)原的方法要比均值濾波器、中值濾波器等其他方法的噪聲去除效果要好，帶阻濾波器更多地保存了高頻細(xì)節(jié)信息。圖7示出本文超分辨率復(fù)原結(jié)果與其他插值方法的結(jié)果對(duì)比。

(a) 原圖像(b) 雙三次插值(c) 均勻B樣條曲面插值(d) 本文超分辨率復(fù)原

圖7 超分辨率復(fù)原結(jié)果與其他插值方法對(duì)比

Fig.7 Comparison between non-uniform surface B-spline interpolation and other interpolation methods

從圖7可以看出，直觀評(píng)價(jià)本文方法的復(fù)原效果更清晰。再對(duì)比細(xì)節(jié)的提升情況。分別截取原圖像、雙三次內(nèi)插、均勻B樣條曲面插值和非均勻B樣條曲面插值結(jié)果的局部(圖7(a)中的黑色方框區(qū)域)放大圖進(jìn)行對(duì)比(圖8)，可以看出本文方法對(duì)圖像細(xì)節(jié)的提升效果較好；而雙三次插值的方法降低了圖像的清晰度；均勻B樣條曲面插值的方法則給圖像帶來了許多嚴(yán)重噪聲。

(a) 原圖像 (b) 雙三次插值(c) 均勻B樣條曲面插值 (d) 本文超分辨率重建

圖8 細(xì)節(jié)對(duì)比圖

(圖7(a)中的黑色方框區(qū)域放大圖)

Fig.8 Comparison in details

3 超分重建效果評(píng)價(jià)

用平均梯度來表征圖像的清晰度，反映圖像中的微小細(xì)節(jié)反差和紋理變換特征；平均梯度越大，則圖像的清晰度越高，微小細(xì)節(jié)及紋理反映越好。圖像峰值信噪比(peak signal to noise ratio，PSNR)[18]反映了圖像信噪比變化情況的統(tǒng)計(jì)平均，PSNR值越大，代表失真越少。圖像熵反映了圖像的信息量大小，熵值越大，表明信息量越大?；谌祟愐曈X系統(tǒng)(human vision system，HVS)適于提取目標(biāo)結(jié)構(gòu)信息的特點(diǎn)，Wang等[19]提出了圖像結(jié)構(gòu)相似度(structural similarity image model，SSIM)的概念。該指標(biāo)反映了高分辨率圖像與超分辨率復(fù)原圖像的結(jié)構(gòu)差異，SSIM的值越大，表明超分辨率重構(gòu)效果越好。但是，由于遙感圖像的超分辨率復(fù)原沒有供參考用的高分辨率圖像，本文構(gòu)造了無參考圖像結(jié)構(gòu)相似度(no-reference SSIM，NR-SSIM)的評(píng)價(jià)方法： 對(duì)待評(píng)價(jià)圖像進(jìn)行低通濾波得到1景參考圖像，然后基于2次模糊的清晰度算法(Re-Blur)[20]來比較參考退化圖像與待評(píng)價(jià)圖像的SSIM值——如果1景圖像已經(jīng)模糊了，那么對(duì)它再進(jìn)行1次模糊處理后，高頻分量變化不大，SSIM值較大；但如果原圖像是清晰的，對(duì)它再進(jìn)行1次模糊處理后，則高頻分量變化會(huì)非常大，SSIM值反而較小。因此NR-SSIM的計(jì)算結(jié)果越小，表明圖像越清晰，反之則越模糊。

將原圖像、基于雙三次插值和均勻節(jié)點(diǎn)“邊界”條件的B樣條曲面插值結(jié)果圖像與本文方法重建的圖像進(jìn)行對(duì)比，各評(píng)價(jià)指標(biāo)的對(duì)比結(jié)果見表2。

表2 不同插值方法評(píng)價(jià)結(jié)果對(duì)比Tab.2 Comparison of evaluation results for different interpolation methods

從表2可以看出，用本文方法重建的圖像，不僅清晰度、信噪比和信息量都有所提高，而且NR-SSIM評(píng)價(jià)也是最優(yōu)的?？梢姡诒疚姆椒ǖ膱D像超分辨率重建使圖像更加清晰、信息量更大、信噪比更高，可以應(yīng)用于實(shí)際遙感圖像復(fù)原中。

本文使用截止頻率[21]來評(píng)價(jià)圖像空間分辨率提高的程度。用刀刃法(edge method)[22]計(jì)算調(diào)制傳遞函數(shù)(modulation transfer function，MTF)，當(dāng)MTF曲線不再明顯下降而開始進(jìn)入震蕩階段時(shí)，所對(duì)應(yīng)的頻率被稱為截止頻率；而截止頻率的倒數(shù)反映了能分辨出圖像細(xì)節(jié)的最小距離，即空間分辨率[23]，二者的對(duì)比結(jié)果見表3。

表3 截止頻率與空間分辨率對(duì)比Tab.3 Comparison between cut-off frequency and spatial resolution

從截止頻率的評(píng)價(jià)結(jié)果可以看出，用本文方法重建的遙感圖像使圖像的空間分辨率比原圖像提高了1.43倍，比用雙三次插值和均勻節(jié)點(diǎn)“邊界”條件B樣條曲面插值復(fù)原的遙感圖像的空間分辨率都有較大的提高。

4 結(jié)論

本文從低成本“軟”的圖像處理技術(shù)出發(fā)，研究提高遙感圖像空間分辨率的方法，得到如下結(jié)論：

1)參考亞像元參叉線陣的思想，提出一種參叉采樣的超分辨率復(fù)原方法。將2景不同時(shí)刻獲取的同一地區(qū)的圖像經(jīng)過灰度匹配和亞像元級(jí)的精確幾何配準(zhǔn)后，采用參叉交錯(cuò)像元采樣到原來圖像網(wǎng)格2倍的網(wǎng)格中，使2景圖像的信息互補(bǔ)，更有助于圖像細(xì)節(jié)的提升。

2)結(jié)合非均勻B樣條曲面插值算法，將重采樣的圖像置于三維空間中，對(duì)沒有值的坐標(biāo)處用三次B樣條插值計(jì)算像元值。構(gòu)造6像元×6像元的插值曲面網(wǎng)絡(luò)點(diǎn)鄰域，每一個(gè)待求像元值都是由鄰域擬合的曲面產(chǎn)生。真實(shí)圖像的待求像元的鄰域并不一定是均勻分布，因此本文在節(jié)點(diǎn)參數(shù)化時(shí)選用非均勻方法更接近真實(shí)情況。非均勻的參數(shù)化方法比均勻的參數(shù)化方法在實(shí)現(xiàn)圖像插值上略微復(fù)雜，在處理方法上本文選擇了自適應(yīng)參數(shù)求解方法。在反算雙三次插值控制曲面時(shí)使用“非節(jié)點(diǎn)”條件求解的方法，形成自適應(yīng)的控制曲面，無需人為賦值，也可以避免賦值的不科學(xué)性。

3)在求解待插值點(diǎn)參數(shù)值時(shí)使用了平行線法和黃金分割法，迭代尋找最優(yōu)值。通過與雙三次插值方法和均勻節(jié)點(diǎn)“邊界”條件B樣條曲面插值法對(duì)比，本文方法效果均為最優(yōu)。通過細(xì)節(jié)放大圖像的視覺評(píng)價(jià)亦可看出，本文方法對(duì)提升圖像的細(xì)節(jié)信息有很大幫助。另外，本文算法的魯棒性良好，適用于各種遙感圖像。

[1] Lange D A,Vu P,Wang S C H,et al.6000-element infrared focal plane array for reconnaissance applications[C]//The International Society for Optical Engineering.Denver:SPIE,1999.

[2] 劉軍,張永生,范永弘.ADS40機(jī)載數(shù)字傳感器的攝影測(cè)量處理與應(yīng)用[J].測(cè)繪學(xué)院學(xué)報(bào),2002,19(3):186-194. Liu J,Zhang Y S,Fan Y H.The photogrammetric processing and applications of the ADS40 airborne digital sensor[J].Journal of Institute of Surveying and Mapping,2002,19(3):186-194.

[3] Andrews H C,Hunt B R.Digital image restoration[C]//Prentice-Hall Signal Processing Series.Englewood Cliffs,NJ:Prentice-Hall,1977:136-140.

[4] Panda R,Chatterji B N.Least squares generalized B-spline signal and image processing[J].Signal Processing,2001,81(10):2005-2017.

[5] Unser M,Aldroubi A,Eden M.Fast B-spline transforms for continuous image representation and interpolation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1991,13(3):277-285.

[6] 馮杰飛.基于非均勻B樣條曲面的圖像插值算法[D].濟(jì)南:山東經(jīng)濟(jì)學(xué)院,2010. Feng J F.The Image Interpolation Algorithm Based on the Non-uniform B-spline Surface[D].Jinan:Shandong Economic University,2010.

[7] 安博文,薛冰玢.光纖耦合的亞像元圖像超分辨率重建[J].計(jì)算機(jī)工程,2012,38(17):222-225. An B W,Xue B B.Super-resolution reconstruction for sub-pixel image with fiber bundle coupling[J].Computer Engineering,2012,38(17):222-225.

[8] 胡萌萌.基于亞像素級(jí)位移的超分辨率成像系統(tǒng)設(shè)計(jì)[D].西安:西安工業(yè)大學(xué),2012. Hu M M.Image Super-resolution Reconstruction System Design Based on Subpixel Displacement[D].Xi’an:Xi’an Technology University,2012.

[9] 周峰,王世濤,王懷義.關(guān)于亞像元成像技術(shù)幾個(gè)問題的探討[J].航天返回與遙感,2002,23(4):26-33. Zhou F,Wang S T,Wang H Y.Study of several points about subpixel imaging technology[J].Spacecraft Recovery and Remote Sensing,2002,23(4):26-33.

[10]Xiong Y H,Li G Q,Mao A H.Convergence analysis for B-spline geometric interpolation[J].Computers and Graphics,2012,36(7):884-891.

[11]Piegl L,Tiller W.The NURBS Book[M].Berlin:Springer-Verlag,1997.

[12]施法中.計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條[M].北京:北京航空航天大學(xué)出版社,1994. Shi F Z.Computer Aided Geometric Design and Non-nuiform Rational B-Spline[M].Beijing:Beijing University of Aeronautics and Astronautics Press,1994.

[13]王飛.三次B樣條反算的一種簡(jiǎn)便算法[J].北京郵電大學(xué)學(xué)報(bào),1996,19(3):82-88. Wang F.An algorithm for the inverse calculation of the triple B-Spline[J].Journal of Beijing University of Posts and Telecommunications,1996,19(3):82-88.

[14]王開榮.最優(yōu)化方法[M].北京:科學(xué)出版社,2012:94-95. Wang K R.Optimization Methods[M].Beijing:Science Press,2012:94-95.

[15]Vandewalle P,Süsstrunk S,Vetterli M.A frequency domain approach to registration of aliased images with application to super-resolution[J].EURASIP Journal on Advances in Signal Processing,2006:071459.

[16]Stone H S,Wolpov R.Blind cross-spectral image registration using prefiltering and Fourier-based translation detection[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(3):637-650.

[17]Gonzalez R C,Woods R E.數(shù)字圖像處理[M].阮秋琦,阮宇智等譯.北京:電子工業(yè)出版社,2007:194-195. Gonzalez R C,Woods R E.Digital Image Processing[M].Beijing:Publishing House of Electronics Industry,2007:194-195.

[18]Kotwal K,Chaudhuri S.A novel approach to quantitative evaluation of hyperspectral image fusion techniques[J].Information Fusion,2013,14(1):5-18.

[19]Wang Z,Bovik A C,Sheikh H R,et al.Image quality assessment:From error visibility to structural similarity[J].IEEE Transactions on Image Processing,2004,13(4):600-612.

[20]李祚林,李曉輝,馬靈玲,等.面向無參考圖像的清晰度評(píng)價(jià)方法研究[J].遙感技術(shù)與應(yīng)用,2011,26(2):239-246. Li Z L,Li X H,Ma L L,et al.Research of definition assessment based on no-reference digital image quality[J].Remote Sensing Technology and Application,2011,26(2):239-246.

[21]李亞鵬,何斌.采用MTF定量評(píng)估CCD錯(cuò)位成像的成像質(zhì)量[J].紅外與激光工程,2013,42(2):443-448. Li Y P,He B.Quantitative evaluation of image quality of CCD subpixel imaging using MTF[J].Infrared and Laser Engineering,2013,42(2):443-448.

[22]Greer P B,Van D T.Evaluation of an algorithm for the assessment of the MTF using an edge method[J].Medical Physics,2000,27(9):2048-2059.

[23]Pan Z X,Huang H J,Yu J,et al.Super-resolution method based on CS and structural self-similarity for remote sensing images[J].Signal Processing,2012,28(6):859-872.

(責(zé)任編輯： 劉心季)

Super-resolution reconstruction of remote sensing image based on staggered pixels and non-uniform B-spline curved surface

WANG Jingmeng1，2，3，4， ZHANG Aiwu1，2，3，4， MENG Xiangang1，2，3，4， LIU Zhao1，2，3，4

(1.CollegeofResourcesEnvironmentandTourism，CapitalNormalUniversity，Beijing100048，China；2.Laboratoryof3DInformationAcquisitionandApplication，MOST，Beijing100048，China；3.BeijingMunicipalKeyLaboratoryofResourcesEnvironmentandGIS，Beijing100048，China；4.StateKeyLaboratoryIncubationBaseofUrbanEnvironmentalProcessesandDigitalSimulation，Beijing100048，China)

The adoption of information from two images can improve resolution of images; nevertheless, the method based on interpolation is apt to make fussy boundary or cause partial loss of detailed information. To solve this problem, the authors employed an approach to super-resolution reconstruction based on staggered pixels and non-uniform B-spline surface. Two frames of images with low resolution，which were matched and geometrically registered at sub-pixel level，were adopted，and then these two images were staggered and re-sampled to double times of the origin grid；the positions with non-value were interpolated and filled with tri-B-spline surface interpolation，and non-uniform node parametric method was adopted. And the curved surface was composed of 36 known neighborhood pixels. For solving the value of interpolated points，the authors introduced parallel and golden section methods to iterate searching for the optimal value，which made interpolation more accurate; the last but not unimportant, the interpolated images were restored to reconstruct better visualization high resolution images. The assessment of the results demonstrates that this approach can considerably improve definition, quantity of information, signal-to-noise ratio and resolution.

image super-resolution restoration；non-uniform B-spline curved surface；staggered pixels

2013-11-06；

2013-12-24

國家“十二五”科技支撐課題“超高空飛艇載荷集成與定量處理技術(shù)”(編號(hào)： 2012BAH31B01)和北京市自然科學(xué)基金重點(diǎn)項(xiàng)目(B類)“基于斜模成像的高光譜圖像超分重建方法研究”(編號(hào)： KZ201310028035)共同資助。

10.6046/gtzyyg.2015.01.06

王京萌，張愛武，孟憲剛，等.基于參叉像元和非均勻B樣條曲面的遙感圖像超分辨率重建[J].國土資源遙感，2015，27(1)：35-43.(Wang J M，Zhang A W，Meng X G，et al.Super-resolution reconstruction of remote sensing image based on staggered pixels and non-uniform B-spline curved surface[J].Remote Sensing for Land and Resources，2015，27(1)：35-43.)

TP 751.1

A

1001-070X(2015)01-0035-09

王京萌(1986-)，女，博士研究生，主要從事圖像超分辨率重建、高光譜圖像應(yīng)用等方面研究。Email： bj_wjm1018@163.com。