濾光片陣列多光譜圖像拼接與灰度調(diào)整算法

李銅哨， 孫文邦*， 趙漢東， 白新偉， 王志磊

(1. 空軍航空大學(xué)航空作戰(zhàn)勤務(wù)學(xué)院， 長(zhǎng)春 130022； 2. 陸軍勤務(wù)學(xué)院， 重慶 400050； 3. 93525部隊(duì)， 日喀則 857000)

航空濾光片陣列多光譜相機(jī)同一時(shí)刻獲取的數(shù)據(jù)為包含多個(gè)波段的窄帶圖像，需要對(duì)不同時(shí)刻獲取的圖像數(shù)據(jù)進(jìn)行裁剪、拼接處理，才能得到單波段大區(qū)域圖像。受機(jī)械快門影響，相機(jī)在不同時(shí)刻成像時(shí)的曝光量存在一定差異，導(dǎo)致單波段拼接圖像容易出現(xiàn)具有一定寬度的明暗條紋，破壞了地物灰度一致性，影響了光譜圖像后期處理及應(yīng)用。因此，為滿足濾光片陣列多光譜圖像應(yīng)用需求，必須解決濾光片陣列多光譜圖像條帶間存在的地物灰度差異。

目前，圖像灰度校正算法大致可分為兩種類型：一是在預(yù)處理階段對(duì)圖像進(jìn)行校正，如Retinex算法[1-2]、Wallis算法[3]、直方圖匹配算法[4]等；二是在圖像拼接為全景圖時(shí)，對(duì)重疊區(qū)域像素灰度進(jìn)行融合，如漸入漸出融合[5]、加權(quán)平均融合[6]、弧函數(shù)權(quán)重模型[7]等。上述算法主要從視覺(jué)效果上消除圖像灰度差異，而光譜圖像灰度校正需要考慮校正前后地物光譜的保真度，盡可能避免光譜信息的損失。此外，Tian等[8]提出一種改進(jìn)Wallis濾波方法，通過(guò)逐行/列計(jì)算圖像行/列重疊區(qū)域像素灰度均值比，進(jìn)一步將圖像灰度調(diào)整為一致，但灰度校正過(guò)程屬于強(qiáng)制改正圖像灰度，并不適用于包含地物光譜信息的多光譜圖像。為解決分視場(chǎng)多光譜圖像的灰度差異，方秀秀等[9]提出灰度線性變換算法，通過(guò)同名點(diǎn)灰度計(jì)算圖像灰度變換關(guān)系。但該方法的灰度處理效果依賴于同名點(diǎn)數(shù)量及分布，同名點(diǎn)數(shù)量少且分布不均時(shí)，灰度變換關(guān)系不一定能夠準(zhǔn)確反映條帶圖像間灰度變換關(guān)系。

因此，現(xiàn)提出基于各波段重疊區(qū)域灰度均值比的均值統(tǒng)調(diào)算法，利用加速穩(wěn)健特征(speeded up robust features，SURF)算法對(duì)投影后圖像進(jìn)行配準(zhǔn)，單波段圖像重疊區(qū)域灰度平均值比的均值作為圖像灰度調(diào)整系數(shù)，并以參考圖像灰度為基準(zhǔn)，依次對(duì)配準(zhǔn)后圖像灰度進(jìn)行調(diào)整，以解決單波段多光譜拼接圖像中的灰度差異問(wèn)題，且最大限度減少地物光譜信息的損失。

1 濾光片陣列多光譜成像特點(diǎn)

濾光片陣列多光譜相機(jī)采用畫幅式成像，其在電荷耦合元件(charge-coupled device，CCD)探測(cè)器前沿著垂直飛行方向放置一個(gè)鍍有數(shù)個(gè)窄帶帶通濾光片的濾光板，獲得多光譜圖像。

1.1 圖像幾何特點(diǎn)

由于每個(gè)帶通濾光片只能通過(guò)一個(gè)指定波長(zhǎng)范圍的圖像，探測(cè)器被分割成若干個(gè)光譜帶，各光譜帶之間存在一定寬度的隔離帶，以避免波段圖像之間的相互干擾。通過(guò)平臺(tái)的運(yùn)動(dòng)，可以獲得某一區(qū)域內(nèi)的多光譜圖像數(shù)據(jù)，而大區(qū)域單波段多光譜圖像則需要通過(guò)圖像拼接的方式才能夠得到，如圖1所示。

圖1中所示的是以波段2和波段6為例顯示拼接過(guò)程。4個(gè)時(shí)刻t1、t2、t3、t4分別獲得4幅含有8個(gè)波段的多條帶圖像，然后從4幅圖像中裁切各波段條帶圖像，再拼接在一起形成單波段圖像。

圖1 濾光片陣列多光譜數(shù)據(jù)面陣處理Fig.1 Area array processing of filter array multi-spectral dataset

1.2 圖像灰度特點(diǎn)

一般情況下，多光譜相機(jī)一次曝光就可以獲得某一區(qū)域多個(gè)波段大區(qū)域灰度一致圖像，如線性漸變?yōu)V光片式多光譜相機(jī)[10]、可調(diào)諧濾光片式多光譜相機(jī)[11]，但濾光片陣列多光譜相機(jī)一次曝光只能獲得某一區(qū)域不同波段的窄條帶圖像，再通過(guò)裁剪拼接獲得單波段大區(qū)域圖像。

由于不同時(shí)刻成像時(shí)曝光量可能存在差異，獲取的部分圖像會(huì)偏亮(暗)。因此，經(jīng)裁剪、拼接得到的單波段大區(qū)域圖像可能存在條帶狀的灰度差異。例如，圖1所示的t3時(shí)刻與t1、t2、t4時(shí)刻存在明顯的曝光量差異，t3時(shí)刻成像總體上偏暗，則裁切拼接后的每個(gè)單波段圖像中均出現(xiàn)較暗的條帶。正是由于多光譜圖像中可能存在的條帶灰度差異，破壞了地物灰度一致性，破壞了地物目標(biāo)的光譜特征信息準(zhǔn)確度。

2 灰度差異調(diào)整

濾光片陣列多光譜圖像灰度差異調(diào)整主要包括圖像幾何關(guān)系確定、圖像灰度調(diào)整和圖像拼接3個(gè)步驟。

2.1 圖像幾何關(guān)系確定

圖像幾何關(guān)系確定主要包括條帶圖像有效區(qū)域模板確定、圖像投影變換和圖像配準(zhǔn)3個(gè)部分。

2.1.1 條帶圖像有效區(qū)域模板確定

以8譜段濾光片陣列多光譜相機(jī)為例，其多光譜圖像數(shù)據(jù)如圖2(a)所示。條帶圖像間存在過(guò)渡區(qū)域，如圖2(a)中黃色框線區(qū)域和圖2(b)中黑色區(qū)域。因此，需要計(jì)算各波段條帶圖像有效區(qū)域的4個(gè)頂點(diǎn)坐標(biāo)，以構(gòu)建多條帶圖像有效區(qū)域模板，如圖2(b)所示。

圖2 各條帶圖像有效區(qū)域Fig.2 Effective area of strip image

2.1.2 圖像投影變換

濾光片陣列多光譜相機(jī)成像時(shí)除了記錄圖像數(shù)據(jù)外，還記錄了成像時(shí)平臺(tái)的姿態(tài)信息(俯仰角φ、翻滾角ω、航向角κ等)。因此，選擇t1、t2、…、tn時(shí)刻的n幅序列多條帶圖像，并根據(jù)式(1)和式(2)對(duì)各時(shí)刻圖像和對(duì)應(yīng)的圖像模板分別進(jìn)行投影[12]。

(1)

(2)

式中：X、Y為投影后像點(diǎn)坐標(biāo)，像素；H為航高，m；D為投影后圖像地面空間分辨率，m；x、y為條帶圖像上像點(diǎn)坐標(biāo)，mm；f為相機(jī)焦距，mm；a1、a2、a3、b1、b2、b3、c1、c2、c3由式(2)計(jì)算得到。

2.1.3 圖像配準(zhǔn)

當(dāng)多光譜相機(jī)成像時(shí)記錄的平臺(tái)姿態(tài)信息精度足夠高時(shí)，投影變換后的圖像在空間上只有平移關(guān)系。但是，由于平臺(tái)姿態(tài)信息均存在一定的誤差，若利用平臺(tái)姿態(tài)信息進(jìn)行投影后圖像直接進(jìn)行拼接處理，則拼接圖像中很容易出現(xiàn)地物未對(duì)齊的現(xiàn)象，如圖3所示。

從圖3中可以看出，拼接縫處的地物存在明顯錯(cuò)位，如圖中紅色虛線圓圈區(qū)域。因此，通過(guò)提取投影后圖像間的匹配點(diǎn)，并且通過(guò)匹配點(diǎn)坐標(biāo)計(jì)算單應(yīng)性轉(zhuǎn)換矩陣H，對(duì)圖像進(jìn)行配準(zhǔn)處理。主要步驟如下。

圖3 偏移量拼接圖像Fig.3 Offset stitching image

步驟1匹配點(diǎn)提取。目前文獻(xiàn)中，匹配點(diǎn)的提取算法較多，例如，SIFT算法[13]、SURF算法[14]等。

由于同一區(qū)域的各單波段條帶圖像需要由多張濾光片陣列式多光譜序列圖像經(jīng)裁切、拼接得到，處理的數(shù)據(jù)量較大。利用SIFT算法雖然能夠提取亞像素級(jí)的特征點(diǎn)，但特征提取與匹配的耗時(shí)過(guò)長(zhǎng)[15]。因此，現(xiàn)采用運(yùn)行速度較快的SURF算法完成特征提取和匹配，通過(guò)RANSAC算法篩選出高精度匹配點(diǎn)[16]。

步驟2轉(zhuǎn)換矩陣計(jì)算。設(shè)圖像I0匹配點(diǎn)坐標(biāo)為(X1,Y1)、(X2,Y2)、…、(Xn,Yn)，圖像I1的匹配點(diǎn)坐標(biāo)為(x1,y1)、(x2,y2)、…、(xn,yn)，則圖像I1向圖像I0的坐標(biāo)轉(zhuǎn)換為

(3)

令

則式(3)可寫為

P0=H1P1

(4)

故圖像I1向圖像I0坐標(biāo)轉(zhuǎn)換矩陣H1為

(5)

同理可得，圖像I0向圖像I1坐標(biāo)轉(zhuǎn)換為

(6)

步驟3空間位置確定。圖像拼接過(guò)程如圖4所示。首先確定基準(zhǔn)圖像I0，利用圖像I1匹配點(diǎn)坐標(biāo)P1和圖像I0匹配點(diǎn)坐標(biāo)P0，通過(guò)式(5)計(jì)算圖像I1向圖像I0的單應(yīng)性轉(zhuǎn)換矩陣H1，再通過(guò)式(4)計(jì)算圖像I1在圖像I0坐標(biāo)系下的坐標(biāo)，確定二者之間的位置關(guān)系。

圖4 全局矩陣配準(zhǔn)圖像Fig.4 Global matrix alignment image

再通過(guò)類似方法，運(yùn)用式(7)計(jì)算圖像I2在圖像I1坐標(biāo)系中坐標(biāo)，公式為

P1=H2P2

(7)

再將式(7)代入式(4)中可以計(jì)算得到待拼接圖像I2在基準(zhǔn)圖像I0坐標(biāo)系中坐標(biāo)，如式(8)所示，可以確定圖像I2與圖像I0間的位置關(guān)系，表達(dá)式為

P0=H1H2P2

(8)

類似地，若有n張待拼接圖像，則第n張待拼接圖像In與基準(zhǔn)圖像I0間的位置關(guān)系為

(9)

由于采用單應(yīng)性轉(zhuǎn)換矩陣進(jìn)行配準(zhǔn)處理，存在誤差積累效應(yīng)，并且隨著拼接數(shù)量增加，拼接圖像畸變會(huì)越來(lái)越明顯[17]。

故為減小圖像拼接過(guò)程中畸變，選擇序列中間圖像為基準(zhǔn)，采用兩邊向中間拼接的策略，如圖5所示。

設(shè)圖5中基準(zhǔn)圖像I0右側(cè)有n張待拼接圖像(I1,I2, …,In)，第n張待拼接圖像In與基準(zhǔn)圖像I0間的位置關(guān)系如式(9)所示。

設(shè)基準(zhǔn)圖像左側(cè)有m張待拼接圖像(I′1,I′2, …,I′m)，左側(cè)圖像坐標(biāo)矩陣設(shè)為P′i，坐標(biāo)轉(zhuǎn)換矩陣為Ti，與右側(cè)轉(zhuǎn)換矩陣計(jì)算式(4)類似，可以得到左側(cè)待拼接圖像I′1向基準(zhǔn)圖像I0坐標(biāo)轉(zhuǎn)換為

P0=T1P′1

(10)

則左側(cè)第m張待拼接圖像I′m與基準(zhǔn)圖像I0間的位置關(guān)系為

(11)

從圖5可以看出，采用兩側(cè)向中間基準(zhǔn)圖像進(jìn)行配準(zhǔn)處理，則可以明顯減少配準(zhǔn)誤差的積累效應(yīng)。

圖5 圖像轉(zhuǎn)換矩陣計(jì)算Fig.5 Image conversion matrix calculation

2.2 圖像灰度調(diào)整

圖像灰度調(diào)整主要依據(jù)條帶圖像重疊區(qū)域的灰度進(jìn)行調(diào)整，涉及關(guān)鍵技術(shù)主要包括條帶圖像重疊區(qū)域確定、灰度調(diào)整系數(shù)計(jì)算、調(diào)整圖像灰度3個(gè)方面。

2.2.1 條帶圖像重疊區(qū)域確定

以相鄰兩幅相鄰的多光譜圖像為例。首先，利用成像時(shí)平臺(tái)姿態(tài)信息對(duì)多光譜圖像和對(duì)應(yīng)的模板均進(jìn)行投影變換；然后，利用計(jì)算得到的坐標(biāo)轉(zhuǎn)換矩陣Hi(或Ti)對(duì)投影后多光譜圖像和模板與基準(zhǔn)圖像的進(jìn)行配準(zhǔn)，模板與基準(zhǔn)圖像配準(zhǔn)后的位置關(guān)系如圖6(a)、圖6(b)所示。提取單波段圖像重疊區(qū)域時(shí)，只保留模板圖像中對(duì)應(yīng)波段區(qū)域，如圖6(c)、圖6(d)所示的為第一波段模板區(qū)域圖像。將相鄰的轉(zhuǎn)換后的模板圖像相乘，得到相應(yīng)波段重疊區(qū)域模板，如圖6(e)所示。

圖6 圖像重疊區(qū)域模板Fig.6 Image overlay area template

然后，將重疊區(qū)域模板分別與兩幅配準(zhǔn)后多光譜圖像相乘得到單波段的重疊區(qū)域圖像。以第一條帶圖像為例，如圖7所示。圖7(a)、圖7(b)分別為配準(zhǔn)后的第1、2幅多光譜圖像，圖7(c)為第1、2幅中第1波段的重疊區(qū)域圖像。

圖7 重疊區(qū)域圖像Fig.7 Image of overlapping area

2.2.2 灰度調(diào)整系數(shù)計(jì)算

灰度調(diào)整系數(shù)計(jì)算方式較多，例如文獻(xiàn)[8]針對(duì)濾光片陣列多光譜圖像提出灰度線性變換算法，核心思想是分別選擇各波段重疊區(qū)域匹配點(diǎn)的灰度均值比，作為待拼接圖像各波段的灰度調(diào)系數(shù)。但該方法灰度處理效果依賴于同名點(diǎn)的數(shù)量及分布，同名點(diǎn)數(shù)量少且分布不均時(shí)，灰度變換關(guān)系不一定能夠準(zhǔn)確反映條帶圖像間灰度變換關(guān)系。

故采用重疊區(qū)域灰度均值比的均值作為各波段統(tǒng)一的調(diào)整系數(shù)。設(shè)圖7(c)、圖7(d)分別是相鄰第i和第j兩幅圖像中第1個(gè)波段的重疊區(qū)域圖像，可以分別計(jì)算該重疊區(qū)域圖像的灰度均值ki1、kj1。同樣，可以計(jì)算得到其他7個(gè)波段的灰度均值，每幅圖像分別可以得到8個(gè)均值，即ki1、ki2、ki3、ki4、ki5、ki6、ki7、ki8和kj1、kj2、kj3、kj4、kj5、kj6、kj7、kj8。則可以通過(guò)式(12)計(jì)算得到第j幅圖像的灰度調(diào)整為第i幅圖像灰度的調(diào)整系數(shù)gij，即

(12)

式(12)中：m為波段數(shù)。

2.2.3 調(diào)整圖像灰度

求出圖像灰度調(diào)整系數(shù)gij后，就可以對(duì)投影后圖像進(jìn)行灰度調(diào)整，以9幅多光譜圖像為例，選擇中間的第5幅圖像為基準(zhǔn)，依次將各灰度調(diào)整系數(shù)與對(duì)應(yīng)序列圖像中的每個(gè)像素灰度相乘。各灰度調(diào)整系數(shù)計(jì)算公式為

(13)

2.3 圖像拼接

濾光片陣列式多光譜圖像拼接不同于其他圖像，圖像拼接前需要提取各單波段條帶圖像，如圖8所示。首先提取各單波段圖像在配準(zhǔn)后的圖像模板上對(duì)應(yīng)區(qū)域，即構(gòu)建配準(zhǔn)后圖像對(duì)應(yīng)的單波段圖像模板，并將配準(zhǔn)后的圖像與各單波段圖像模板相乘，得到單波段配準(zhǔn)后圖像。最后，對(duì)各單波段條帶圖像進(jìn)行拼接，得到灰度一致的濾光片陣列多光譜單波段圖像。

圖8 單波段圖像拼接Fig.8 Single-band image stitching

3 實(shí)驗(yàn)與分析

為了驗(yàn)證灰度調(diào)整算法的可行性，采用某航拍無(wú)人機(jī)濾光片陣列多光譜相機(jī)拍攝某地區(qū)的多光譜圖像數(shù)據(jù)進(jìn)行驗(yàn)證實(shí)驗(yàn)。其中，部分多光譜圖像如圖9所示。

圖9 實(shí)驗(yàn)圖像Fig.9 Experimental images

未進(jìn)行灰度處理和本文算法處理后的單波段多光譜圖像如圖10所示。

圖10 灰度調(diào)整對(duì)比圖Fig.10 Grayscale adjustment comparison chart

由圖10可知，未進(jìn)行灰度處理的單波段拼接圖像整體灰度不均勻，同一區(qū)域中的地物存在著較為明顯的灰度差異，如圖10中箭頭所指位置。而本文算法處理后的單波段圖像灰度整體比較均勻，不存在明顯的灰度差異，說(shuō)明本文方法較好地解決了條帶圖像灰度不均的問(wèn)題。

為進(jìn)一步說(shuō)明本文算法的有效性，使用文獻(xiàn)[9]中灰度一致性校正算法與本文算法進(jìn)行對(duì)比實(shí)驗(yàn)。文獻(xiàn)[9]與本文算法處理后的各單波段圖像如圖11所示，對(duì)比圖像分別為3波段、4波段、5波段、7波段多光譜圖像。

從目視效果上看，圖11(g)和圖11(h)整體灰度均勻，無(wú)明顯的灰度差異，說(shuō)明對(duì)比算法能夠解決圖像灰度不均勻問(wèn)題。但圖11(f)中地物灰度不均勻問(wèn)題有所改善，并未完全解決灰度不均勻問(wèn)題。而本文算法處理后圖像無(wú)明顯的灰度差異，而且與灰度未調(diào)整圖像差異較小。

圖11 灰度調(diào)整后效果對(duì)比Fig.11 Comparison of the effect after grayscale adjustment

對(duì)于第2波段圖像，由于無(wú)法提取匹配點(diǎn)，文獻(xiàn)[9]無(wú)法進(jìn)行灰度調(diào)整；對(duì)于第4波段圖像，由于部分相鄰圖像匹配點(diǎn)數(shù)量較少，灰度變換模型無(wú)法正確反映匹配點(diǎn)對(duì)應(yīng)圖像灰度的變化關(guān)系，調(diào)整后圖像仍存在灰度差異，部分第4波段圖像的匹配點(diǎn)分布如圖12所示。

圖12 第4波段匹配點(diǎn)分布Fig.12 Distribution of matching points for 4th band

此外，對(duì)比算法調(diào)整后的第7波段圖像中圖像上方較圖像下方灰度明顯偏亮，經(jīng)分析發(fā)現(xiàn)，使用RANSAC算法篩選出的精匹配點(diǎn)具有隨機(jī)性[18]，基于匹配點(diǎn)灰度值計(jì)算出的灰度變換模型系數(shù)可能存在一定差異，影響了條帶圖像灰度的正確調(diào)整，如圖13所示。

圖13 灰度調(diào)整后局部放大對(duì)比圖Fig.13 Grayscale adjusted partial comparison

因此，當(dāng)重疊區(qū)域圖像存在提取匹配點(diǎn)數(shù)量少、匹配地分布不均勻時(shí)，對(duì)比算法無(wú)法解決單波段圖像灰度不一致問(wèn)題。此外，匹配點(diǎn)處的灰度值直接影響了圖像灰度變換關(guān)系參數(shù)，不一定準(zhǔn)確反映相鄰圖像間的灰度變換關(guān)系。而本文算法處理后單波段圖像灰度分布均勻，圖像灰度不一致問(wèn)題得到較好的解決。

為進(jìn)一步說(shuō)明本文算法的性能，分別計(jì)算調(diào)整后各單波段條帶圖像重疊區(qū)域灰度平均值的差值來(lái)說(shuō)明，如圖14所示(僅顯示3、4、5波段圖像)。

圖14橫坐標(biāo)是圖像幅序號(hào)，縱坐標(biāo)是重疊區(qū)域平均灰度差，灰度級(jí)為0～255。由圖14可知，文獻(xiàn)[9]處理后部分單波段條帶圖像重疊區(qū)域灰度均值的變化劇烈，說(shuō)明調(diào)整后單波段圖像部分區(qū)域仍存在較明顯灰度差異。

圖14 重疊區(qū)域灰度均值差值對(duì)比Fig.14 Difference in mean grey value of overlapping areas

總體上看，本文算法處理后的重疊區(qū)域灰度差值要小于文獻(xiàn)[9]。經(jīng)分析可知，對(duì)單波段圖像進(jìn)行灰度調(diào)整時(shí)，當(dāng)重疊區(qū)域匹配點(diǎn)數(shù)量較多時(shí)，文獻(xiàn)[9]可以取得較好的灰度調(diào)整效果。但由于各單波段圖像調(diào)整時(shí)選擇的參考圖像不是同一時(shí)刻獲取的圖像數(shù)據(jù)，灰度調(diào)整后的圖像灰度存在偏暗(亮)的現(xiàn)象。

而本文算法處理后的圖像重疊區(qū)域平均灰度差未處理圖像整體較小，各波段重疊區(qū)域灰度差整體近似一條直線，不存在部分條帶圖像重疊區(qū)域灰度差突變的情況。本文算法處理后的各單波段圖像與參考圖像的灰度近似一致，較好地保持了地物在不同波段下的光譜反射特性。

由此可知，提出的基于重疊區(qū)域圖像的灰度調(diào)整算法，充分利用了重疊區(qū)域圖像的灰度信息，調(diào)整后的各單波段圖像灰度均勻，有效解決了單波段圖像中灰度不一致問(wèn)題。

4 結(jié)論

針對(duì)濾光片陣列式多光譜圖像灰度不一致的問(wèn)題，提出了一種基于重疊區(qū)域的圖像灰度調(diào)整算法，即根據(jù)各個(gè)波段圖像重疊區(qū)域灰度均值比的平均值，對(duì)圖像數(shù)據(jù)中各波段統(tǒng)一進(jìn)行灰度調(diào)整，調(diào)整后的拼接圖像對(duì)齊效果較好、同一地物灰度一致，降低了對(duì)圖像處理中特征提取與匹配的影響。通過(guò)實(shí)驗(yàn)對(duì)比，表明針對(duì)濾光片陣列多光譜圖像提出的灰度調(diào)整算法較目前算法有較強(qiáng)的適用性和較好的灰度處理效果，有效解決了濾光片陣列多光譜圖像中的地物灰度不一致問(wèn)題，能夠滿足多光譜圖像的后期處理及應(yīng)用需求。