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

(1.齊魯工業(yè)大學(xué)(山東省科學(xué)院)，山東省科學(xué)院海洋儀器儀表研究所，山東省海洋監(jiān)測(cè)儀器裝備技術(shù)重點(diǎn)實(shí)驗(yàn)室，國(guó)家海洋監(jiān)測(cè)設(shè)備工程技術(shù)研究中心，山東 青島 266061；2.青島鐳測(cè)創(chuàng)芯科技有限公司，山東 青島 266102)

隨著大視場(chǎng)光學(xué)和定標(biāo)技術(shù)的發(fā)展，推掃型成像光譜儀能逐漸克服視場(chǎng)小、波段受限、定標(biāo)困難等缺點(diǎn)，充分體現(xiàn)結(jié)構(gòu)簡(jiǎn)單、體積小、靈敏度高、光譜分辨率高和地面分辨率高的優(yōu)勢(shì)，在遙感應(yīng)用中成為主流發(fā)展方向[1]。推掃式光譜儀中的面陣CCD探測(cè)器單次探測(cè)到的是一條狹帶物體的光譜，其特點(diǎn)是平行于狹縫方向?yàn)楠M帶物體的灰度分布，而垂直于狹縫方向?yàn)橄裨诠庾V上的展開，如果要完成對(duì)二維物體的光譜成像，需要將各個(gè)狹帶依次成像而后拼接為二維圖像。由于搭載光譜儀的飛行平臺(tái)姿態(tài)和飛行速度的不斷變化，獲取的狹帶圖像并非理想飛行狀態(tài)下的正射影像，因此，拼接之前需要對(duì)狹帶圖像進(jìn)行幾何校正，為狹帶匹配精確的地理信息。

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

為了解決上述問題，本文根據(jù)單應(yīng)映射原理，建立了光譜儀傾斜和正射狀態(tài)下像點(diǎn)的映射關(guān)系，利用GPS/INS組合導(dǎo)航數(shù)據(jù)校正狹帶影像中的畸變，然后拼接成一幅完整的影像，并在ENVI二次開發(fā)平臺(tái)上實(shí)現(xiàn)了推掃狹帶影像的自動(dòng)校正和拼接。

1 推掃影像畸變?cè)?/h2>

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

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

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

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

圖2 歐拉角和三維坐標(biāo)系旋轉(zhuǎn)Fig.2 Euler angles and three-dimensional coordinate system rotation

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

2 單應(yīng)映射

。

(1)

則可以將單應(yīng)性表示為：

(2)

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

圖3 平面單應(yīng)性Fig.3 Homography of the planes

物理變換部分是觀測(cè)到的圖像平面相關(guān)的旋轉(zhuǎn)R和平移t的影響之和，即成像儀的外參數(shù)，表示為：

(3)

成像儀內(nèi)參數(shù)矩陣M表示為：

(4)

其中，fx,fy為焦距，cx,cy為主點(diǎn)坐標(biāo)，則單應(yīng)性就可以表示為：

(5)

因此，已知成像儀內(nèi)外參數(shù)，就可以獲得單應(yīng)矩陣H。內(nèi)參數(shù)由成像儀標(biāo)定獲得，外參數(shù)由GPS/INS測(cè)量系統(tǒng)獲取的歐拉角和位移量計(jì)算得到。

qb=sbMWbQ。

(6)

同理，正射時(shí)地面點(diǎn)Q與成像平面點(diǎn)qn之間的單應(yīng)關(guān)系滿足：

qn=snMWnQ。

(7)

qb=MRbRn-1M-1qn，

(8)

其中單應(yīng)矩陣H=MRbRn-1M-1。

。

(9)

那么，同一個(gè)成像中心獲取的兩幅狹帶圖像之間的平面坐標(biāo)轉(zhuǎn)換關(guān)系可以表示為：

(10)

3 IDL實(shí)現(xiàn)

IDL是美國(guó)ITT VIS公司推出的第四代交互式、跨平臺(tái)、面向矩陣處理的編程語言，具有快速的數(shù)據(jù)分析、圖像處理和強(qiáng)大的可視化功能[16-17]。采用IDL語言調(diào)用ENVI平臺(tái)中的圖像處理函數(shù)，可以很方便地進(jìn)行二次開發(fā)，實(shí)現(xiàn)遙感數(shù)據(jù)的快速分析和可視化。推掃圖像的自動(dòng)拼接主要包括如下3個(gè)基本步驟：

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

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

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

所以，對(duì)于狹帶影像上的每一個(gè)二維點(diǎn)(xb,yb)，都可以獲得校正后的對(duì)應(yīng)點(diǎn)(xn,yn)，點(diǎn)(xn,yn)的灰度值即為點(diǎn)(xb,yb)的灰度值。

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

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

帶有地理信息的狹帶通過鑲嵌批處理函數(shù)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;圖像鑲嵌。

拼接后的影像認(rèn)為是光譜儀理想姿態(tài)下獲取的正射影像，具有與GPS獲取的一致的位置信息，拼接影像點(diǎn)的高光譜曲線與原始掃描行對(duì)應(yīng)點(diǎn)的一致，能夠真實(shí)地反映地面的空間特征和光譜特征。

4 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

5 結(jié)語

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