基于ENVI二次開發的高光譜推掃圖像拼接技術

蓋穎穎，蓋志剛，禹定峰，劉恩曉，李輝，秦勝光

(1.齊魯工業大學(山東省科學院)，山東省科學院海洋儀器儀表研究所，山東省海洋監測儀器裝備技術重點實驗室，國家海洋監測設備工程技術研究中心，山東 青島 266061；2.青島鐳測創芯科技有限公司，山東 青島 266102)

隨著大視場光學和定標技術的發展，推掃型成像光譜儀能逐漸克服視場小、波段受限、定標困難等缺點，充分體現結構簡單、體積小、靈敏度高、光譜分辨率高和地面分辨率高的優勢，在遙感應用中成為主流發展方向[1]。推掃式光譜儀中的面陣CCD探測器單次探測到的是一條狹帶物體的光譜，其特點是平行于狹縫方向為狹帶物體的灰度分布，而垂直于狹縫方向為像在光譜上的展開，如果要完成對二維物體的光譜成像，需要將各個狹帶依次成像而后拼接為二維圖像。由于搭載光譜儀的飛行平臺姿態和飛行速度的不斷變化，獲取的狹帶圖像并非理想飛行狀態下的正射影像，因此，拼接之前需要對狹帶圖像進行幾何校正，為狹帶匹配精確的地理信息。

推掃狹帶圖像幾何校正通常有兩種方式，一種是基于地面控制點，結合插值方法對圖像進行地理化，如Strakhov等[2]基于簡化的相機運動物理模型，提出了一種新的基于控制點的推掃圖像幾何校正算法。另外一種是采用GPS/INS組合導航系統同步獲取每條狹帶的姿態和位置信息，通過后處理算法生成具有地理信息的圖像[3]，如馬偉波等[4]基于POS導航數據解算推掃高光譜影像成像瞬間外方位元素，通過改進的共線方程實現了幾何校正；Xu等[5]基于對推掃成像過程的分解，建立了機載長線性多元素掃描成像系統的幾何校正模型，提出了線性方程校正算法和一種新的切線校正算法。本文基于第二種處理方式，在處理過程中，需要進行高光譜圖像讀取、投影轉換、圖像插值和圖像拼接等多個任務，因此采用遙感圖像處理平臺(the environment for visualizing images, ENVI)進行二次開發，結合交互式數據語言(interactive data language, IDL)，實現校正拼接的自動化和快速可視化。ENVI是基于IDL語言開發的，對IDL進行了一定程度的封裝，一個通用的函數就可以完成圖像處理的復雜任務，還支持分塊讀取和處理，在遙感應用中使用廣泛[6]。如Jeong[7]采用ENVI和IDL對Landsat衛星影像進行了海表面溫度的反演；Che等[8]在ENVI平臺中采用IDL語言實現了高光譜礦物信息識別模塊的開發；李耀輝[9]利用ENVI二次開發平臺實現了風云遙感農情監測原型系統的開發。但是目前還沒有針對推掃狹帶影像拼接的二次開發模塊，無法滿足遙感推掃影像自動化處理業務需求。

為了解決上述問題，本文根據單應映射原理，建立了光譜儀傾斜和正射狀態下像點的映射關系，利用GPS/INS組合導航數據校正狹帶影像中的畸變，然后拼接成一幅完整的影像，并在ENVI二次開發平臺上實現了推掃狹帶影像的自動校正和拼接。

1 推掃影像畸變原因

推掃式成像是利用飛行平臺的向前運動，借助于與飛行方向垂直的掃描線記錄而構成二維圖像。推掃型成像光譜儀通常采用一個垂直于運動方向的面陣CCD來感應地面響應，在飛行平臺向前運動中完成二維空間掃描，平行于平臺運動方向，通過光柵和棱鏡分光，完成光譜維掃描[10]，如圖1所示。因此，CCD上一個點對應一個譜段，一條線對應一個譜面，CCD探測器每次成像是空間一條線上的光譜信息。為了獲得空間二維圖像，再通過機械推掃，完成整個平面的圖像和光譜數據采集。

推掃成像時，CCD探測器所記錄的高光譜圖像數據是沿著飛行方向的條幅。由于搭載光譜儀的飛行平臺在飛行過程中，不能一直保證理想的姿態正射獲取影像，速度和姿態的不穩定導致飛行平臺的位置、航偏角、俯仰角和橫滾角不斷隨機變化，引起光譜儀拍攝時外方位元素也不斷隨機變化。因此，CCD曝光時每條掃描線對應的光譜儀外方位元素不一致引起了圖像的幾何畸變。

圖1 推掃成像示意圖Fig.1 Schematic diagram of push-broom imaging

飛行平臺姿態不穩定造成地面掃描行之間相互交錯，圖像扭曲變形，影響后期地物目標的解析和判別[11]。為了描述飛行平臺的姿態以便校正圖像，引入了機體坐標系oxbybzb和地理坐標系oxnynzn，如圖2所示。由于姿態變化引起兩坐標系之間產生的角位移由航偏角ψ、俯仰角θ和橫滾角φ表示，也稱為歐拉角[12]。

圖2 歐拉角和三維坐標系旋轉Fig.2 Euler angles and three-dimensional coordinate system rotation

飛行平臺速度不穩定易造成掃描行之間的行間距忽大忽小，出現重疊或間隙[11]，為了獲得地面的完整影像，通常推掃成像需保證一定的采樣率，因此在圖像拼接時就需要借助GPS位置信息，對重疊的掃描行進行幾何糾正和圖像融合處理。

2 單應映射

。

(1)

則可以將單應性表示為：

(2)

其中，s是任意尺度的比例，H包含物體平面的物理變換和成像儀內參數矩陣兩部分。

圖3 平面單應性Fig.3 Homography of the planes

物理變換部分是觀測到的圖像平面相關的旋轉R和平移t的影響之和，即成像儀的外參數，表示為：

(3)

成像儀內參數矩陣M表示為：

(4)

其中，fx,fy為焦距，cx,cy為主點坐標，則單應性就可以表示為：

(5)

因此，已知成像儀內外參數，就可以獲得單應矩陣H。內參數由成像儀標定獲得，外參數由GPS/INS測量系統獲取的歐拉角和位移量計算得到。

qb=sbMWbQ。

(6)

同理，正射時地面點Q與成像平面點qn之間的單應關系滿足：

qn=snMWnQ。

(7)

qb=MRbRn-1M-1qn，

(8)

其中單應矩陣H=MRbRn-1M-1。

。

(9)

那么，同一個成像中心獲取的兩幅狹帶圖像之間的平面坐標轉換關系可以表示為：

(10)

3 IDL實現

IDL是美國ITT VIS公司推出的第四代交互式、跨平臺、面向矩陣處理的編程語言，具有快速的數據分析、圖像處理和強大的可視化功能[16-17]。采用IDL語言調用ENVI平臺中的圖像處理函數，可以很方便地進行二次開發，實現遙感數據的快速分析和可視化。推掃圖像的自動拼接主要包括如下3個基本步驟：

(1)影像和GPS/INS數據讀取。遙感影像數據包含圖像本身和頭文件，ENVI二次開發提供了函數讀取遙感影像及其屬性，如ENVI_OPEN_FILE、ENVI_FILE_QUERY、ENVI_GET_SLICE等。GPS/INS數據存儲于文本文件中，按照文本文件讀取方式即可獲得狹帶影像獲取時光譜儀的姿態和位置信息。

(2)單應矩陣計算和單應映射。以北東地坐標系為地理坐標系，依據公式(10)計算得到單應矩陣H，由于需要將機體坐標系下的平面點坐標計算至北東地坐標系下，因此旋轉矩陣Cbn需求逆，主要代碼命令如下：

H=M_inv # MATRIX_POWER(C,-1) # M;計算單應矩陣。

所以，對于狹帶影像上的每一個二維點(xb,yb)，都可以獲得校正后的對應點(xn,yn)，點(xn,yn)的灰度值即為點(xb,yb)的灰度值。

(3)圖像拼接。校正后的每條狹帶圖像中心點的二維地理坐標即光譜儀成像中心的GPS二維坐標，根據光譜儀的成像地面分辨率，選定影像投影方式，可以為每條狹帶設置地理信息，主要代碼命令如下：

map_info = ENVI_MAP_INFO_CREATE(/geographic, mc=mc, ps=ps);為狹帶添加地理信息。

帶有地理信息的狹帶通過鑲嵌批處理函數MOSAIC_DOIT可以拼接成為一幅完整的影像，主要代碼如下：

ENVI_DOIT, 'MOSAIC_DOIT', fid=mosaic_fids, pos=corr_pos, dims=corr_dims, out_name=outname_mosaic, r_fid=out_mosaic_fid, xsize=xsize, ysize=ysize, x0=x0, y0=y0, georef=1, map_info=mosaic_map_info, out_dt=2, pixel_size=ps, see_through_val=see_through_val, use_see_through=use_see_through;圖像鑲嵌。

拼接后的影像認為是光譜儀理想姿態下獲取的正射影像，具有與GPS獲取的一致的位置信息，拼接影像點的高光譜曲線與原始掃描行對應點的一致，能夠真實地反映地面的空間特征和光譜特征。

4 實驗結果

本文選擇河北省保定市郊區的高光譜影像進行校正拼接實驗，影像由搭載于無人機的Pika L高光譜儀拍攝獲取。Pika L高光譜儀由美國Resonon公司設計生產，光譜范圍為400～1000 nm，光譜分辨率為2.1 nm，CCD掃描行寬度為900像素。飛行過程中同時搭載法國SBG Ellipse2-D慣導系統，實時獲取光譜儀的姿態位置信息。高光譜儀將推掃獲取的原始狹帶影像先簡單拼接起來，存儲于固態硬盤中，此時的地理信息并未經過糾正，圖像存在幾何畸變，圖4a所示為原始圖像的假彩色圖像。狹帶經過幾何校正和拼接后才能正確顯示地面目標的特征，如圖4b所示。

為了能夠定量檢驗該幾何校正方法的效果，同時采用Pika L光譜儀自帶的Georectify Airborne Datacube軟件對原始影像進行幾何校正，將兩種方法得到的校正影像進行比較。兩種校正方法均采用UTM投影，以WGS-84為基準面。首先在軟件校正影像中隨機選取10個均勻分布的明顯地物點，讀取其坐標值，作為采樣點用于評定校正精度，然后從本文方法校正后影像中讀取其相應坐標值，經過對10個采樣點殘差的計算得到如表1所示的精度檢驗結果。

圖4 原始影像和校正影像Fig.4 Original image and corrected image

表1 幾何校正隨機采樣精度分析Table 1 Accuracy analysis of random sampling in geometric correction

表1中北向距離均方根誤差為0.632 7 m,東向距離均方根誤差為0.381 7 m，總均方根誤差為0.738 9 m。由表1可以看出，采樣點在x和y方向上的坐標偏移均不超過1 m，兩種方法得到的校正圖像地理信息較為接近；y方向坐標均方根誤差大于x方向坐標均方根誤差，即像點坐標的經度值準確性高于緯度值。對于某些地理精度要求不高的航空高光譜遙感應用來說，本方法取得的校正效果已滿足需求，如果需要進一步提高精度，可以通過增加地面控制點或與高精度地圖進行圖像配準實現幾何精校正。

5 結語

本文根據推掃成像和單應映射原理，結合GPS/INS組合導航系統實時獲取光譜儀姿態角度和位置信息，在ENVI二次開發平臺上，采用IDL語言實現了高光譜儀推掃狹帶影像的自動校正和拼接。驗證實驗表明，本方法與自帶軟件校正拼接效果接近，均方根誤差基本滿足一般的高光譜遙感應用。雖然本文方法能夠取得較為理想的校正拼接效果，但是單掃描行的校正過程耗時較長，無法實時獲取校正影像，下一步將就提高校正拼接效率展開更加深入的研究。另外，拼接過程中不同成像條件下的勻色處理同樣是后續需要研究的內容。