濾光片陣列多光譜圖像拼接與灰度調整算法

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

(1. 空軍航空大學航空作戰勤務學院， 長春 130022； 2. 陸軍勤務學院， 重慶 400050； 3. 93525部隊， 日喀則 857000)

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

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

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

1 濾光片陣列多光譜成像特點

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

1.1 圖像幾何特點

由于每個帶通濾光片只能通過一個指定波長范圍的圖像，探測器被分割成若干個光譜帶，各光譜帶之間存在一定寬度的隔離帶，以避免波段圖像之間的相互干擾。通過平臺的運動，可以獲得某一區域內的多光譜圖像數據，而大區域單波段多光譜圖像則需要通過圖像拼接的方式才能夠得到，如圖1所示。

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

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

1.2 圖像灰度特點

一般情況下，多光譜相機一次曝光就可以獲得某一區域多個波段大區域灰度一致圖像，如線性漸變濾光片式多光譜相機[10]、可調諧濾光片式多光譜相機[11]，但濾光片陣列多光譜相機一次曝光只能獲得某一區域不同波段的窄條帶圖像，再通過裁剪拼接獲得單波段大區域圖像。

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

2 灰度差異調整

濾光片陣列多光譜圖像灰度差異調整主要包括圖像幾何關系確定、圖像灰度調整和圖像拼接3個步驟。

2.1 圖像幾何關系確定

圖像幾何關系確定主要包括條帶圖像有效區域模板確定、圖像投影變換和圖像配準3個部分。

2.1.1 條帶圖像有效區域模板確定

以8譜段濾光片陣列多光譜相機為例，其多光譜圖像數據如圖2(a)所示。條帶圖像間存在過渡區域，如圖2(a)中黃色框線區域和圖2(b)中黑色區域。因此，需要計算各波段條帶圖像有效區域的4個頂點坐標，以構建多條帶圖像有效區域模板，如圖2(b)所示。

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

2.1.2 圖像投影變換

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

(1)

(2)

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

2.1.3 圖像配準

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

從圖3中可以看出，拼接縫處的地物存在明顯錯位，如圖中紅色虛線圓圈區域。因此，通過提取投影后圖像間的匹配點，并且通過匹配點坐標計算單應性轉換矩陣H，對圖像進行配準處理。主要步驟如下。

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

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

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

步驟2轉換矩陣計算。設圖像I0匹配點坐標為(X1,Y1)、(X2,Y2)、…、(Xn,Yn)，圖像I1的匹配點坐標為(x1,y1)、(x2,y2)、…、(xn,yn)，則圖像I1向圖像I0的坐標轉換為

(3)

令

則式(3)可寫為

P0=H1P1

(4)

故圖像I1向圖像I0坐標轉換矩陣H1為

(5)

同理可得，圖像I0向圖像I1坐標轉換為

(6)

步驟3空間位置確定。圖像拼接過程如圖4所示。首先確定基準圖像I0，利用圖像I1匹配點坐標P1和圖像I0匹配點坐標P0，通過式(5)計算圖像I1向圖像I0的單應性轉換矩陣H1，再通過式(4)計算圖像I1在圖像I0坐標系下的坐標，確定二者之間的位置關系。

圖4 全局矩陣配準圖像Fig.4 Global matrix alignment image

再通過類似方法，運用式(7)計算圖像I2在圖像I1坐標系中坐標，公式為

P1=H2P2

(7)

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

P0=H1H2P2

(8)

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

(9)

由于采用單應性轉換矩陣進行配準處理，存在誤差積累效應，并且隨著拼接數量增加，拼接圖像畸變會越來越明顯[17]。

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

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

設基準圖像左側有m張待拼接圖像(I′1,I′2, …,I′m)，左側圖像坐標矩陣設為P′i，坐標轉換矩陣為Ti，與右側轉換矩陣計算式(4)類似，可以得到左側待拼接圖像I′1向基準圖像I0坐標轉換為

P0=T1P′1

(10)

則左側第m張待拼接圖像I′m與基準圖像I0間的位置關系為

(11)

從圖5可以看出，采用兩側向中間基準圖像進行配準處理，則可以明顯減少配準誤差的積累效應。

圖5 圖像轉換矩陣計算Fig.5 Image conversion matrix calculation

2.2 圖像灰度調整

圖像灰度調整主要依據條帶圖像重疊區域的灰度進行調整，涉及關鍵技術主要包括條帶圖像重疊區域確定、灰度調整系數計算、調整圖像灰度3個方面。

2.2.1 條帶圖像重疊區域確定

以相鄰兩幅相鄰的多光譜圖像為例。首先，利用成像時平臺姿態信息對多光譜圖像和對應的模板均進行投影變換；然后，利用計算得到的坐標轉換矩陣Hi(或Ti)對投影后多光譜圖像和模板與基準圖像的進行配準，模板與基準圖像配準后的位置關系如圖6(a)、圖6(b)所示。提取單波段圖像重疊區域時，只保留模板圖像中對應波段區域，如圖6(c)、圖6(d)所示的為第一波段模板區域圖像。將相鄰的轉換后的模板圖像相乘，得到相應波段重疊區域模板，如圖6(e)所示。

圖6 圖像重疊區域模板Fig.6 Image overlay area template

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

圖7 重疊區域圖像Fig.7 Image of overlapping area

2.2.2 灰度調整系數計算

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

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

(12)

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

2.2.3 調整圖像灰度

求出圖像灰度調整系數gij后，就可以對投影后圖像進行灰度調整，以9幅多光譜圖像為例，選擇中間的第5幅圖像為基準，依次將各灰度調整系數與對應序列圖像中的每個像素灰度相乘。各灰度調整系數計算公式為

(13)

2.3 圖像拼接

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

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

3 實驗與分析

為了驗證灰度調整算法的可行性，采用某航拍無人機濾光片陣列多光譜相機拍攝某地區的多光譜圖像數據進行驗證實驗。其中，部分多光譜圖像如圖9所示。

圖9 實驗圖像Fig.9 Experimental images

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

圖10 灰度調整對比圖Fig.10 Grayscale adjustment comparison chart

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

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

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

圖11 灰度調整后效果對比Fig.11 Comparison of the effect after grayscale adjustment

對于第2波段圖像，由于無法提取匹配點，文獻[9]無法進行灰度調整；對于第4波段圖像，由于部分相鄰圖像匹配點數量較少，灰度變換模型無法正確反映匹配點對應圖像灰度的變化關系，調整后圖像仍存在灰度差異，部分第4波段圖像的匹配點分布如圖12所示。

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

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

圖13 灰度調整后局部放大對比圖Fig.13 Grayscale adjusted partial comparison

因此，當重疊區域圖像存在提取匹配點數量少、匹配地分布不均勻時，對比算法無法解決單波段圖像灰度不一致問題。此外，匹配點處的灰度值直接影響了圖像灰度變換關系參數，不一定準確反映相鄰圖像間的灰度變換關系。而本文算法處理后單波段圖像灰度分布均勻，圖像灰度不一致問題得到較好的解決。

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

圖14橫坐標是圖像幅序號，縱坐標是重疊區域平均灰度差，灰度級為0～255。由圖14可知，文獻[9]處理后部分單波段條帶圖像重疊區域灰度均值的變化劇烈，說明調整后單波段圖像部分區域仍存在較明顯灰度差異。

圖14 重疊區域灰度均值差值對比Fig.14 Difference in mean grey value of overlapping areas

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

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

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

4 結論

針對濾光片陣列式多光譜圖像灰度不一致的問題，提出了一種基于重疊區域的圖像灰度調整算法，即根據各個波段圖像重疊區域灰度均值比的平均值，對圖像數據中各波段統一進行灰度調整，調整后的拼接圖像對齊效果較好、同一地物灰度一致，降低了對圖像處理中特征提取與匹配的影響。通過實驗對比，表明針對濾光片陣列多光譜圖像提出的灰度調整算法較目前算法有較強的適用性和較好的灰度處理效果，有效解決了濾光片陣列多光譜圖像中的地物灰度不一致問題，能夠滿足多光譜圖像的后期處理及應用需求。