基于觀測角信息的HJ-1A/B衛(wèi)星光學(xué)影像幾何精糾正

楊 亮， 賈 益， 江萬壽， 張 過

(1.武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室,武漢 430079； 2.長光衛(wèi)星技術(shù)有限公司,長春 130000)

0 引言

雖然環(huán)境一號(hào)(HJ-1)A/B星已發(fā)射近9 a，但其存檔的海量CCD影像并未得到充分利用，幾何定位不準(zhǔn)確及山區(qū)地形畸變誤差是限制其應(yīng)用的主要瓶頸之一。針對(duì)這一大規(guī)模數(shù)據(jù)的應(yīng)用需求，中科院遙感所、武漢大學(xué)等單位已經(jīng)研制了環(huán)境衛(wèi)星影像幾何精糾正系統(tǒng)，采用全球ETM+影像作為參考數(shù)據(jù)進(jìn)行自動(dòng)處理[1-4]。

目前，環(huán)境衛(wèi)星對(duì)外提供的標(biāo)準(zhǔn)產(chǎn)品是二級(jí)影像產(chǎn)品，因?yàn)槲聪匦纹鸱挠绊?，給后續(xù)的變化檢測、定量遙感等應(yīng)用帶來了很大的麻煩，需要進(jìn)一步做幾何精糾正處理[1-6]?，F(xiàn)有的二級(jí)影像數(shù)據(jù)幾何精糾正方法可分為2種： ①假設(shè)可以通過影像匹配獲取大量的可靠控制點(diǎn)，采用多項(xiàng)式擬合[2]或小面元[3-6]進(jìn)行影像幾何配準(zhǔn)糾正； ②沿襲一級(jí)影像糾正的思路，先利用分布均勻且具有較好高程范圍的控制點(diǎn)建立一個(gè)近似的成像模型(如直接線性變換模型[6]或仿射變換模型[7-8])，再進(jìn)行正射糾正。其中，第1種方法的多項(xiàng)式模型無法表達(dá)地形起伏引起的畸變，三角面元的方法則需要大量分布均勻的控制點(diǎn)以反映地形的變化[7]，因而無法適應(yīng)控制點(diǎn)難以大量獲取和地形起伏較大地區(qū)的影像精糾正要求； 第2種方法在理論上比較合理，但也需要均勻分布的三維控制點(diǎn)，難以普遍使用。因此，對(duì)于全球規(guī)模的定量遙感而言，急需尋找一種能適合于全球不同地區(qū)的影像精糾正方法。

在“863”項(xiàng)目“星機(jī)地綜合定量遙感系統(tǒng)與應(yīng)用示范”幾何歸一化子課題研究中，筆者注意到中低空間分辨率的遙感影像數(shù)據(jù)產(chǎn)品一般都附帶有與影像對(duì)應(yīng)的衛(wèi)星觀測角文件，該文件實(shí)際隱含了原始影像的成像模型信息，卻沒有被充分利用。針對(duì)這一情況，本文提出一種利用觀測角信息重建CCD二級(jí)影像“成像模型”的方法，并利用90 m數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù)進(jìn)行正射糾正。

1 常用幾何糾正模型分析與對(duì)比

在衛(wèi)星光學(xué)影像的幾何糾正過程中，需要建立影像的行列號(hào)與對(duì)應(yīng)地面點(diǎn)之間的數(shù)學(xué)關(guān)系，然后利用該數(shù)學(xué)關(guān)系確定糾正影像的像元在原始影像中的位置，最后進(jìn)行灰度重采樣，計(jì)算該點(diǎn)的像元值賦給糾正影像。除采用共線方程的嚴(yán)格模型外，常用的幾何糾正模型包括仿射變換模型、多項(xiàng)式變換模型、直接線性變換(direct liner transformation，DLT)模型和有理函數(shù)模型(rational polynomial coefficient, RPC)[8-14]。其中，DLT模型和RPC模型可用于對(duì)多線陣共線方程進(jìn)行近似，而仿射變換模型和多項(xiàng)式變換模型可用于影像配準(zhǔn)或?qū)LT模型和RPC模型進(jìn)行像方改正。

1.1 像方改正/影像配準(zhǔn)模型

1.1.1 二維仿射變換模型

二維仿射變換模型將二維控制點(diǎn)與二維影像點(diǎn)以一種最簡單的數(shù)學(xué)關(guān)系表示出來，由于需要的控制點(diǎn)少，所以在影像幾何糾正中經(jīng)常被采用。但該模型沒有考慮高程引起的變形，因此不能用于衛(wèi)星影像傾角較大或地形起伏較明顯的情況。目前，該模型在有嚴(yán)格模型或替代模型的情況下，常用于像方改正，即

(1)

式中： (x,y)為影像像元坐標(biāo)； (X,Y)為控制點(diǎn)坐標(biāo)(一般為參考影像坐標(biāo))；ki(i=1,2,…，6)為系數(shù)。

1.1.2 二維多項(xiàng)式變換模型

與仿射變換模型相比，多項(xiàng)式變換模型能改正影像的平移、旋轉(zhuǎn)、縮放和扭曲等全局變形，在基于影像配準(zhǔn)的幾何糾正中經(jīng)常被采用； 其缺點(diǎn)是無法表達(dá)地形起伏變形等局部突變。其數(shù)學(xué)表達(dá)式為

(2)

1.2 近似成像模型

1.2.1 擴(kuò)展仿射變換模型

擴(kuò)展仿射變換模型考慮了地形的影響[8-9]，適合于視場角較小、近似平行投影的情況。由于需要擬合高程起伏的影響，對(duì)控制點(diǎn)的高程取值范圍有一定的要求。其數(shù)學(xué)表達(dá)式為

(3)

1.2.2 擴(kuò)展DLT模型

與擴(kuò)展仿射變換模型類似，DLT模型將影像的像元坐標(biāo)與地面點(diǎn)的坐標(biāo)進(jìn)行直接對(duì)應(yīng)[6]。為了更好地適應(yīng)大角度多線陣投影，可將分母相同的DLT模型擴(kuò)展為分母不同的DLT模型，其數(shù)學(xué)表達(dá)式可表示為

(4)

式中Li(i=1,2,…，14)為系數(shù)。

該模型對(duì)控制點(diǎn)的高程變化有一定的要求，但與擴(kuò)展仿射變換模型相比，DLT模型可以適應(yīng)視場角更大的情況。

1.2.3 RPC模型

RPC模型[10]將地面點(diǎn)的三維坐標(biāo)(P,L,H)與影像的行列號(hào)(l,s)以比值多項(xiàng)式的形式聯(lián)系起來。RPC模型的表達(dá)式為

(5)

式中Nl,Dl,Ns,Ds分別為4個(gè)三次多項(xiàng)式，即

(6)

(7)

(8)

(9)

式中： (l,s)為去中心歸一化的控制點(diǎn)影像行列號(hào)； (P,L,H)為去中心歸一化的控制點(diǎn)地理坐標(biāo)；ai,bi,ci和di(i=1,2,...,20)為模型的系數(shù)。

因?yàn)橄禂?shù)過多，RPC模型主要作為嚴(yán)格模型的替代模型使用，一般不宜利用地面控制點(diǎn)進(jìn)行擬合。在用于影像糾正時(shí)，RPC模型常與像方改正模型配合使用，可利用控制點(diǎn)進(jìn)一步擬合像方光學(xué)畸變和姿態(tài)變化等誤差。

2 基于觀測角信息的成像模型擬合

2.1 影像幾何誤差來源分析

環(huán)境衛(wèi)星二級(jí)影像已經(jīng)過系統(tǒng)幾何糾正和相對(duì)輻射校正，鏡頭畸變、傳感器平臺(tái)位置與姿態(tài)的變化以及地球曲率對(duì)影像的影響可以忽略不計(jì)[2]，但地形起伏在CCD影像中引起的畸變卻不能忽略不計(jì)，由于采用雙CCD相機(jī)拼接對(duì)地成像，較大的地形起伏在雙拼相機(jī)的影像外側(cè)會(huì)產(chǎn)生很大的畸變。影像成像時(shí)的傳感器天頂角和方位角及與地形起伏的關(guān)系如圖1所示。

圖1 影像的觀測角與高程起伏引起的畸變

圖1中，S為影像的投影中心，β為影像的傳感器天頂角，α為傳感器方位角。其中，β最大可達(dá)35°。根據(jù)幾何關(guān)系可知，高程起伏Δh在影像上引起的畸變?chǔ)為

(10)

在影像最外側(cè)，100 m地形起伏引起的影像畸變可達(dá)57 m，換算到像方達(dá)到1.92個(gè)像元。因此，高差幾km的地區(qū)影像畸變可達(dá)近20個(gè)像元。

2.2 基于觀測角信息的近似成像模型

由于RPC模型解算簡單，能夠很好地逼近嚴(yán)格成像模型，故本文選擇RPC模型來擬合成像模型，并用于環(huán)境衛(wèi)星二級(jí)光學(xué)影像的進(jìn)一步幾何精糾正，以消除地形起伏引起的畸變。

該影像二級(jí)產(chǎn)品的衛(wèi)星角度文件SatAngle.txt包含了傳感器天頂角和方位角的信息?？紤]到數(shù)據(jù)量的問題，SatAngle.txt中的角度信息是33像元×33像元間隔的采樣數(shù)據(jù)。由圖1可知，從每個(gè)采樣點(diǎn)的角度信息可以恢復(fù)該點(diǎn)的原始光束。因此，對(duì)每個(gè)像元的光線在多個(gè)高程面進(jìn)行采樣，就可以獲得RPC模型擬合需要的控制網(wǎng)(圖2)。

圖2 根據(jù)角度信息建立的RPC擬合空間控制格網(wǎng)

利用環(huán)境衛(wèi)星二級(jí)產(chǎn)品自帶的觀測角信息對(duì)RPC系數(shù)進(jìn)行求解的過程如下：

1)讀取角度信息文件并計(jì)算每個(gè)抽樣像元的投影坐標(biāo)。設(shè)有m×n個(gè)抽樣像元的角度記錄，每一行記錄了L1及L2影像中的行列坐標(biāo)、橢球高為0時(shí)的像元經(jīng)緯度坐標(biāo)、傳感器方位角和傳感器高度角。選擇通用橫軸墨卡托投影作為投影坐標(biāo)系，然后由經(jīng)緯度計(jì)算高程為0時(shí)的投影坐標(biāo)(X0,Y0)。

2)根據(jù)經(jīng)緯度坐標(biāo)，從全球1 km格網(wǎng)DEM估計(jì)影像范圍內(nèi)的最小高程Zmin和最大高程Zmax，并分為5個(gè)高程面，每層之間的間隔Δh為(Zmax-Zmin)/4，第i個(gè)高程面的高程為Zmin+iΔh，i=0,...,4。

3)對(duì)于每個(gè)高程面Zi，計(jì)算每個(gè)采樣像元光線對(duì)應(yīng)的投影坐標(biāo)。由(X0,Y0)和式(11)可計(jì)算Zi對(duì)應(yīng)的投影坐標(biāo)(Xi,Yi)，即

(11)

4)通過m×n×5個(gè)地面點(diǎn)坐標(biāo)與m×n個(gè)二級(jí)影像像素坐標(biāo)，解算模型參數(shù)[11]。

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

3.1 實(shí)驗(yàn)數(shù)據(jù)

為驗(yàn)證本文提出的基于觀測角信息的中等空間分辨率衛(wèi)星影像幾何精糾正算法的精度，本文選取2景高程起伏較大地區(qū)的HJ-1A/B CCD影像進(jìn)行實(shí)驗(yàn)。影像1為西藏西南部HJ-1A CCD1影像，拍攝時(shí)間為2015年5月28日，影像成像范圍為E82.7°～87.8°，N28.3°～32.2°，影像大小為16 717像元×14 407像元，共4個(gè)波段，地面海拔范圍為88～6 597 m，平均海拔為3 774 m； 影像2為河北省HJ-1B CCD1影像，拍攝時(shí)間為2016年1月25日，成像范圍為E111.6°～ 117.2°，N36.4°～40.4°，影像大小為16 516像元×14 771像元，共4個(gè)波段，地面海拔范圍為-2～2 836 m，平均海拔為995 m。

3.2 實(shí)驗(yàn)流程

為了對(duì)比本文方法的效果，本文設(shè)計(jì)了3個(gè)實(shí)驗(yàn)方案，對(duì)第2節(jié)介紹的像方改正方法和近似成像模型方法進(jìn)行組合(圖3)。

(a) 方案1和方案2 (b) 方案3

圖33種實(shí)驗(yàn)方案流程

Fig.3Threeprocessingflowsofexperiments

方案1和方案2如圖3(a)所示，采用全球ETM+作為參考影像，通過影像匹配獲取控制點(diǎn)，然后解算變換參數(shù)。方案1是影像配準(zhǔn)的思路，通過多項(xiàng)式擬合影像之間的變形進(jìn)行糾正； 方案2在方案1的基礎(chǔ)上增加高程信息，通過擴(kuò)展DLT模型擬合高程的影響，然后采用正射糾正的方式進(jìn)行糾正。如圖3(b)所示，方案3直接利用衛(wèi)星觀測角信息擬合成像模型，然后利用擬合得到的成像模型和90 m空間分辨率的DEM將參考影像模擬到影像成像時(shí)的大致狀態(tài)，以消除姿態(tài)和高程起伏等引起的畸變，最后與影像進(jìn)行匹配。

實(shí)驗(yàn)中先采用自動(dòng)匹配的方法在30像元×30像元的規(guī)則格網(wǎng)上自動(dòng)匹配控制點(diǎn)，然后通過粗差剔除和人工檢查得到最終的控制點(diǎn)。為便于對(duì)比分析，不同方案采用統(tǒng)一的控制點(diǎn)，參數(shù)解算中不再自動(dòng)剔除控制點(diǎn)。

3.3 結(jié)果與分析

像方改正采用仿射變換、二次多項(xiàng)式和三次多項(xiàng)式進(jìn)行實(shí)驗(yàn)。分別對(duì)一次多項(xiàng)式、二次多項(xiàng)式和三次多項(xiàng)式進(jìn)行了RPC參數(shù)擬合實(shí)驗(yàn)，發(fā)現(xiàn)一次多項(xiàng)式擬合誤差較大，但不影響最終糾正結(jié)果； 二次多項(xiàng)式擬合精度最好，糾正精度最佳； 三次多項(xiàng)式可能存在過擬合，導(dǎo)致最終糾正誤差很大。限于篇幅，像方改正只保留二次多項(xiàng)式和三次多項(xiàng)式的結(jié)果，RPC擬合只保留二次多項(xiàng)式結(jié)果。

影像1經(jīng)匹配得到589個(gè)控制點(diǎn)，未經(jīng)人工編輯，全部參與計(jì)算，控制點(diǎn)的行方向誤差如圖4所示，不同糾正方案的精度結(jié)果如表1所示。

(a) 三次多項(xiàng)式擬合誤差 (b) 擴(kuò)展DLT+二次多項(xiàng)式擬合誤差

(c) 擴(kuò)展DLT+三次多項(xiàng)式擬合誤差 (d) 二次RPC+三次多項(xiàng)式擬合誤差

圖4 不同糾正方法結(jié)果的行方向誤差曲線Fig.4 Error curves in line direction yielded by different rectification methods表1 影像1糾正誤差Tab.1 Rectification error of image 1 (像元)

①mx和my分別為糾正影像在行方向和列方向的中誤差； ② maxVx和maxVy為糾正影像在行方向和列方向的最大絕對(duì)誤差。

從圖4和表1可以看出，未考慮高程影響的多項(xiàng)式糾正影像在行方向存在多達(dá)20個(gè)像元的誤差，在列方向誤差相對(duì)較小，整個(gè)誤差曲線不均衡； 考慮高程影響的正射糾正結(jié)果在行方向的誤差減小到和列方向相當(dāng)，列方向的誤差也較未考慮高程影響的方案有較大的減小。采用二次多項(xiàng)式的誤差曲線在列方向呈現(xiàn)明顯的系統(tǒng)性； 采用三次多項(xiàng)式的誤差比采用二次多項(xiàng)式的誤差有較大的減小，且基本消除了列方向的系統(tǒng)性。對(duì)比由觀測角文件重建RPC模型和控制點(diǎn)擬合DLT的糾正結(jié)果可以看出，兩者的中誤差和最大絕對(duì)誤差比較接近，但是由觀測角文件重建的RPC模型在行方向精度更高，大多數(shù)點(diǎn)的誤差均在2個(gè)像素以內(nèi)。

影像2經(jīng)匹配共得到373個(gè)控制點(diǎn)，因參考影像和HJ-1B影像差異很大，對(duì)匹配點(diǎn)進(jìn)行了人工檢查和編輯。與影像1的實(shí)驗(yàn)結(jié)果類似，考慮高程的糾正結(jié)果能夠很好地消除地形起伏的影響?？赡苡捎谠黾恿巳斯た刂泣c(diǎn)編輯，減少了自動(dòng)匹配控制點(diǎn)存在的誤差，本景影像幾何糾正的中誤差和最大絕對(duì)誤差都有較大幅度的減小(表2)。

表2 影像2糾正誤差Tab.2 Rectification error of image 2 (像元)

4 結(jié)論

針對(duì)常規(guī)中等空間分辨衛(wèi)星影像的幾何糾正方法無法適應(yīng)地形起伏較大區(qū)域的問題，提出了一種基于觀測角信息的環(huán)境衛(wèi)星CCD影像幾何精糾正方法。該方法利用觀測角信息重建影像的“成像模型”并用RPC模型表示，從而可以采用正射糾正的方法對(duì)二級(jí)影像進(jìn)行正射糾正。實(shí)驗(yàn)證明，該方法無需控制點(diǎn)信息即可重建環(huán)境衛(wèi)星光學(xué)二級(jí)影像坐標(biāo)與地面坐標(biāo)的關(guān)系，用于正射糾正的精度明顯高于多項(xiàng)式法，并略優(yōu)于DLT擬合法。該方法的最大優(yōu)勢(shì)在于重建成像模型過程不需參考影像及其控制點(diǎn)，不受控制點(diǎn)平面分布及高程分布的影響，具有較大的普適性。

參考文獻(xiàn)(References)：

[1] 唐 娉,張宏偉,趙永超,等.全球30 m分辨率多光譜影像數(shù)據(jù)自動(dòng)化處理的實(shí)踐與思考[J].遙感學(xué)報(bào),2014,18(2):231-253.

Tang P,Zhang H W,Zhao Y C,et al.Practice and thoughts of the automatic processing of multispectral images with 30 m spatial resolution on the global scale[J].Journal of Remote Sensing,2014,18(2):231-253.

[2] 熊文成,申文明,王 橋,等.環(huán)境一號(hào)衛(wèi)星A/B星數(shù)據(jù)自動(dòng)幾何精校正設(shè)計(jì)研究[J].中國工程科學(xué),2014,16(3):37-42.

Xiong W C,Shen W M,Wang Q,et al.Research on HJ-1A/B satellites data automatic geometric precision correction design[J].Engineering Science,2014,16(3):37-42.

[3] 唐 亮,唐 娉,閻福禮,等.HJ-1 CCD圖像自動(dòng)幾何精糾正系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J].計(jì)算機(jī)應(yīng)用,2012,32(s2):237-241.

Tang L,Tang P,Yan F L,et al.Design and implementation of automatic accurate geometric correction system for HJ-1 CCD images[J].Journal of Computer Applications,2012,32(s2):237-241.

[4] 單小軍,唐 娉,胡昌苗,等.圖像分層匹配的HJ-1A/B CCD影像自動(dòng)幾何精校正技術(shù)與系統(tǒng)實(shí)現(xiàn)[J].遙感學(xué)報(bào),2014,18(2):254-266.

Shan X J,Tang P,Hu C M,et al.Automatic geometric precise correction technology and system based on hierarchical image matching for HJ-1A/B CCD images[J].Journal of Remote Sensing,2014,18(2):254-266.

[5] 胡九超,周忠發(fā).基于快速匹配與線性橡皮拉伸的HJ-1CCD影像幾何精校正[J].湖北農(nóng)業(yè)科學(xué),2015,54(3):703-704,708.

Hu J C,Zhou Z F.HJ-1CCD image geometric accurate rectification based on fast matching and linear rubber sheeting[J].Hubei Agricultural Sciences,2015,54(3):703-704,708.

[6] 李全文,趙衛(wèi)林,楊曉梅,等.基于FAST測點(diǎn)的大幅寬HJ星圖像幾何精糾正方法[J].國土資源遙感,2014,26(4):14-22.doi:10.6046/gtzyyg.2014.04.03.

Li Q W,Zhao W L,Yang X M,et al. AST algorithm-based geometric accurate rectification of large HJ satellite image[J].Remote Sensing for Land and Resources,2014,26(4):14-22.doi:10.6046/gtzyyg.2014.04.03.

[7] 潘建平,郝建明,趙繼萍.基于SURF的圖像配準(zhǔn)改進(jìn)算法[J].國土資源遙感,2017,29(1):110-115.doi:10.6046/gtzyyg.2017.01.17.

Pang J P,Hao J M,Zhao J P.Improved algorithm based on SURF for image registration[J].Remote Sensing for Land and Resources,2017,29(1):110-115.doi:10.6046/gtzyyg.2017.01.17.

[8] Savopol F,Armenakis C.Modelling of the IRS-1C satellite pan stereo-imagery using the DLT approach[C]//Proceedings of 1998 IAPRS Commission IV Symposium on GIS-Between Visions and Applications,Stuttgart:[s.n.],1998:511-514.

[9] Hattori S,Ono T,Raser C F,et al.Orientation of high-resolution satellite images based on affine projection[C]//Proceedings of International Archives of Photogrammetry and Remote Sensing.Amsterdam: ISPRS,2000:359-366.

[10] 張劍清,張祖勛.高分辨率遙感影像基于仿射變換的嚴(yán)格幾何模型[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2002,27(6):555-559.

Zhang J Q,Zhang Z X.Strict geometric model based on affine transformation for remote sensing image with high resolution[J].Geomatics and Information Science of Wuhan University,2002,27(6):555-559.

[11] Tao C V,Hu Y.A comprehensive study of the rational function model for photogrammetric processing[J].Photogrammetric Engineering and Remote Sensing,2001,67(12):1347-1357.

[12] 張 過,李 揚(yáng),祝小勇,等.有理函數(shù)模型在光學(xué)衛(wèi)星影像幾何糾正中的應(yīng)用[J].航天返回與遙感,2010,31(4):51-57.

Zhang G,Li Y,Zhu X Y,et al.Application of RFM in geometric rectification of optical satellite image[J].Spacecraft Recovery and Remote Sensing,2010,31(4):51-57.

[13] 鞏丹超,張永生.有理函數(shù)模型的解算與應(yīng)用[J].測繪學(xué)院學(xué)報(bào),2003,20(1):39-42,46.

Gong D C,Zhang Y S.The solving and application of rational function model[J].Journal of Institute of Surveying and Mapping,2003,20(1):39-42,46.

[14] 楊寶林,呂婷婷,王少軍,等.Pleiades衛(wèi)星圖像正射糾正方法[J].國土資源遙感,2015,27(3):25-29.doi:10.6046/gtzyyg.2015.03.05.

Yang B L,Lyu T T,Wang S J,et al. Ortho-rectification method for Pleiades satellite images[J].Remote Sensing for Land and Resources,2015,27(3):25-29.doi:10.6046/gtzyyg.2015.03.05.