李冰川，潘紅播，田佳，崔澤華

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

0 引言

GLT數(shù)據(jù)是常見(jiàn)于擺掃式衛(wèi)星影像中的定位格網(wǎng)數(shù)據(jù)?；贕LT數(shù)據(jù)的坐標(biāo)反算是衛(wèi)星數(shù)據(jù)幾何校正中的重要環(huán)節(jié)，尋求高質(zhì)量坐標(biāo)反算方法對(duì)擺掃式衛(wèi)星數(shù)據(jù)處理至關(guān)重要。由于GLT模型是由影像坐標(biāo)獲取地理位置的正向索引，沒(méi)有確定的反算函數(shù)，因此基于GLT數(shù)據(jù)的坐標(biāo)反算較為困難。

國(guó)內(nèi)外學(xué)者圍繞坐標(biāo)反算進(jìn)行了一系列研究。丁一帆等[1]根據(jù)單幀影像建立基于視線(xiàn)的幾何模型，針對(duì)每一幀數(shù)據(jù)進(jìn)行嚴(yán)格計(jì)算，建立統(tǒng)一的局部平面模型，將各幀影像進(jìn)行連接，然后再構(gòu)建整體RPC模型進(jìn)行反算。郭廣猛等[2]將GLT作為地面控制點(diǎn)，計(jì)算多項(xiàng)式系數(shù)，建立物像轉(zhuǎn)換多項(xiàng)式模型，實(shí)現(xiàn)坐標(biāo)反算。這類(lèi)算法都是通過(guò)構(gòu)建數(shù)學(xué)模型替代擺掃成像嚴(yán)密模型，建立物像轉(zhuǎn)換關(guān)系實(shí)現(xiàn)反算。但由于擺掃成像的特殊性，使用多項(xiàng)式模型替代該過(guò)程會(huì)引入較大的擬合誤差，導(dǎo)致該類(lèi)方法精度較低。另一類(lèi)方法是先由影像坐標(biāo)計(jì)算到物方坐標(biāo)，再推斷該物方點(diǎn)周?chē)c(diǎn)位對(duì)應(yīng)的像點(diǎn)應(yīng)在該物方點(diǎn)對(duì)應(yīng)像點(diǎn)周?chē)?，然后開(kāi)辟窗口，逐個(gè)搜索，根據(jù)距離確權(quán)，內(nèi)插出像點(diǎn)坐標(biāo)[3-4]。這類(lèi)算法都未考慮在無(wú)參考像點(diǎn)情況下的坐標(biāo)反算，并且有較大計(jì)算冗余，精度較低。賈益等[5]提出一種基于路線(xiàn)追蹤的方法快速定位經(jīng)緯度對(duì)應(yīng)原始影像中位置。文獻(xiàn)[6-7]先去除蝴蝶結(jié)效應(yīng)，然后基于GLT數(shù)據(jù)使用一階線(xiàn)性迭代逼近將輸入物方點(diǎn)反算至像方系統(tǒng)。這些算法都是根據(jù)經(jīng)緯度數(shù)據(jù)變化的連續(xù)性，迭代計(jì)算改正坐標(biāo)向量，逼近像方坐標(biāo)，其精度較之前方法有所提升，但沒(méi)有考慮地形起伏影響。

針對(duì)上述問(wèn)題，本文提出了一種顧及地形效應(yīng)的線(xiàn)性逼近高精度坐標(biāo)反算方法。先用線(xiàn)性逼近法將輸入物方點(diǎn)反算至未經(jīng)地形改正像方點(diǎn)，然后改正由高差引起的像點(diǎn)位移實(shí)現(xiàn)高精度坐標(biāo)反算。該方法能夠克服目前存在的受地形起伏影響大、精度低的缺點(diǎn)。

1 算法原理

蝴蝶結(jié)(bowtie)效應(yīng)又稱(chēng)雙眼皮效應(yīng)，是擺掃式成像中隨著擺掃角度增大，兩邊影像空間分辨率降低，產(chǎn)生的相鄰掃描幀數(shù)據(jù)重復(fù)的現(xiàn)象[8]。在GLT數(shù)據(jù)中表現(xiàn)為從兩邊到中心的逐漸減少的數(shù)據(jù)重疊以及經(jīng)緯度數(shù)據(jù)在不同幀之間的不連續(xù)。蝴蝶結(jié)效應(yīng)的存在，使得線(xiàn)性迭代逼近時(shí)產(chǎn)生多解，降低精度甚至無(wú)法收斂。但單幀數(shù)據(jù)內(nèi)部經(jīng)緯度連續(xù)變化[9]，線(xiàn)性逼近法精度可達(dá)子像素，且收斂較快。因此，先根據(jù)輸入物方經(jīng)度(λ)、緯度(φ)和高程(h)，計(jì)算對(duì)應(yīng)像方點(diǎn)所在幀號(hào)；然后根據(jù)幀號(hào)抽取單幀數(shù)據(jù)；再線(xiàn)性逼近反算像方坐標(biāo)(xu，yu)；最后改正地形起伏引起的像點(diǎn)位移，得到最終的像方坐標(biāo)(xu+dx，yu+dy)。具體步驟如圖1所示。

圖1 坐標(biāo)反算流程

1.1 改正高差引起像點(diǎn)位移

星載擺掃式成像系統(tǒng)是通過(guò)掃描方向的擺掃及沿軌道方向上運(yùn)動(dòng)實(shí)現(xiàn)二維影像的獲取[10]，因此影像的成像視角在掃描向變化較大。為保證平臺(tái)穩(wěn)定性，衛(wèi)星姿態(tài)保持本體坐標(biāo)系與軌道坐標(biāo)系一致，因此可認(rèn)為影像在飛行方向上成像光線(xiàn)近似平行。地形起伏引起的像點(diǎn)位移主要在掃描向，飛行方向上位移相對(duì)較小。

改正高差引起像點(diǎn)位移基本思想是將真實(shí)高程和線(xiàn)性插值高程之間的差值轉(zhuǎn)化為掃描光軸的角度變化量，進(jìn)而計(jì)算像點(diǎn)位移改正。GLT具有一定分辨率的地理定位格網(wǎng)數(shù)據(jù)，基于該數(shù)據(jù)使用線(xiàn)性逼近法可反算得到未經(jīng)地形改正的像方坐標(biāo)。此過(guò)程假設(shè)高程是局部線(xiàn)性變化，但實(shí)際地形變化復(fù)雜，與假設(shè)不相符，因此二者的差值將引起像點(diǎn)位移。線(xiàn)性逼近反算像點(diǎn)坐標(biāo)(xu，yu)一般為浮點(diǎn)型，如圖2所示，其左右臨近像方點(diǎn)為m1、m2，對(duì)應(yīng)的物方點(diǎn)為G1和G2?；贕1、G2及(λ，φ，h)線(xiàn)性插值計(jì)算像點(diǎn)坐標(biāo)時(shí)，本質(zhì)是預(yù)設(shè)待定點(diǎn)的高程等于G1、G2線(xiàn)性插值高程，然而待定點(diǎn)G的實(shí)際高程為h，與G′的高差為Δh。高差Δh的存在直接導(dǎo)致掃描光軸角度變化，進(jìn)而產(chǎn)生像點(diǎn)位移Δd。

由于擺掃式成像系統(tǒng)為等角成像系統(tǒng)，因此像方的偏移Δd可轉(zhuǎn)換為成像視角α上的改變量Δα。擺掃式成像系統(tǒng)掃描向相鄰像元之間的角度差γ可由采樣頻率和掃描速度所確定(或由視場(chǎng)角與影像寬計(jì)算)，飛行向相鄰像元角度差為傳感器探元成像瞬時(shí)視場(chǎng)角(instantaneous field of view，IFoV)。綜上，像點(diǎn)偏移Δd，可由高差Δh、傳感器下視角α、相鄰像元角度差γ或瞬時(shí)視場(chǎng)角IFoV計(jì)算得到。

如圖2所示，S為衛(wèi)星位置，G為衛(wèi)星采樣點(diǎn)，則掃描向天頂角β、星地距離l可由α角與衛(wèi)星高、地球半徑算出。由三角形內(nèi)正弦定理，掃描向角度變化量 Δα求解如式(1)所示。

圖2 地形起伏引起在線(xiàn)性?xún)?nèi)插中產(chǎn)生像點(diǎn)位移示意圖

(1)

而天頂角β可由式(2)求得。

(2)

式中：RN為地球半徑；H為衛(wèi)星高。

衛(wèi)星傳感器下視角α可由下視角在掃描向上分量αx和飛行向上分量αy計(jì)算得到，表達(dá)如式(3)所示。

tan2α=tan2αx+tan2αy

(3)

由于擺掃式衛(wèi)星姿態(tài)角通常較小，αx約為衛(wèi)星姿態(tài)角滾動(dòng)角(ω)與掃描鏡掃描角αm之和，如式(4)所示。

αx≈αm+ω

(4)

掃描鏡掃描角αm由線(xiàn)性逼近所求未經(jīng)地形改正的影像坐標(biāo)(xu，yu)計(jì)算求得，表達(dá)如式(5)所示。

(5)

式中：影像寬為w；傳感器視場(chǎng)角為FOV；掃描角αm正負(fù)代表方向；飛行向傳感器下視角αy約等于俯仰角φ。

將上述變量代入公式，可解得高差引起的像點(diǎn)偏移Δd，表達(dá)如式(6)所示。

(6)

將Δd在飛行向和掃描向進(jìn)行分解，根據(jù)坐標(biāo)幾何關(guān)系，并按照各自的像元分辨率換算，表達(dá)如式(7)所示。

(7)

改正高差引起像點(diǎn)位移后影像定位結(jié)果為(xu+dx，yu+dy)。

1.2 線(xiàn)性逼近法坐標(biāo)反算

基于單幀影像地理信息數(shù)據(jù)進(jìn)行坐標(biāo)反算過(guò)程如圖3所示。線(xiàn)性逼近坐標(biāo)反算思想是根據(jù)GLT數(shù)據(jù)中像方點(diǎn)的二維梯度和輸入物方經(jīng)緯度，任意給定初始像點(diǎn)坐標(biāo)，計(jì)算像方改正量，迭代計(jì)算像點(diǎn)坐標(biāo)。當(dāng)二元函數(shù)模型λ(x，y)、φ(x，y)并未直接給出時(shí)，很難根據(jù)(λ，φ，h)求出像方坐標(biāo)(x，y)?？紤]到幀內(nèi)連續(xù)變化，可根據(jù)所給經(jīng)緯度二維數(shù)據(jù)進(jìn)行線(xiàn)性迭代逼近。將λ、φ泰勒展開(kāi)得式(8)。

(8)

式中：φ′y、φ′x、λ′y、λ′x為該點(diǎn)的一階偏導(dǎo)數(shù)?？梢愿鶕?jù)該像點(diǎn)鄰近像點(diǎn)經(jīng)緯度數(shù)據(jù)計(jì)算梯度，解得如式(9)所示。

(9)

式中：D=λ′xφ′y-φ′xλ′y。預(yù)設(shè)給定任意初始像點(diǎn)坐標(biāo)可由(x0，y0)，再由所給地理坐標(biāo)(λ，φ，h)，以及初始點(diǎn)二維梯度代入式(9)算出對(duì)應(yīng)像方改正量(Δx，Δy)。將(x0，y0)+(Δx，Δy)作為新計(jì)算的像點(diǎn)坐標(biāo)，并繼續(xù)迭代，直到兩次迭代計(jì)算出的坐標(biāo)在一個(gè)GLT網(wǎng)格內(nèi)結(jié)束迭代，再由該GLT網(wǎng)格經(jīng)緯度梯度和所給物方坐標(biāo)線(xiàn)性?xún)?nèi)插出未經(jīng)地形改正的像方坐標(biāo)(xu，yu)。

圖3 迭代逼近示意圖

1.3 反算像方點(diǎn)幀號(hào)

逐幀取每幀影像中間一行地理定位數(shù)據(jù)，組成似GLT二維數(shù)組，該數(shù)據(jù)在飛行方向不存在跳變。本文反算像方點(diǎn)幀號(hào)具體步驟如下。

步驟1：逐幀抽取掃描幀中間一行數(shù)據(jù)，或者掃描幀中間兩行數(shù)據(jù)均值組成似GLT數(shù)組，如MODIS抽取后為203×1 354維數(shù)組。

步驟2：將物方坐標(biāo)(λ，φ，h)與代入步驟1中數(shù)組，經(jīng)線(xiàn)性逼近反算得像方坐標(biāo)。

步驟3：抽取距離步驟2中像點(diǎn)坐標(biāo)最近的兩幀數(shù)據(jù)，分別由線(xiàn)性逼近反算影像坐標(biāo)，計(jì)算像點(diǎn)與各幀中心距離，取距離較短的幀號(hào)作為該物方點(diǎn)對(duì)應(yīng)像方點(diǎn)所在幀號(hào)。

1.4 特殊區(qū)域處理

當(dāng)影像跨越子午線(xiàn)時(shí)，影像覆蓋區(qū)域的經(jīng)度會(huì)從180°變成-180°，導(dǎo)致GLT數(shù)據(jù)出現(xiàn)跳變，若不特殊處理則會(huì)產(chǎn)生粗差，或迭代不收斂。本文算法先判斷該數(shù)據(jù)是否包含子午線(xiàn)部分，若包含則將經(jīng)度全部加上360°以保證經(jīng)度連續(xù)。經(jīng)過(guò)該調(diào)整后對(duì)反算效果不產(chǎn)生影響。

2 實(shí)驗(yàn)設(shè)計(jì)與分析

2.1 實(shí)驗(yàn)設(shè)計(jì)

MODIS的GLT數(shù)據(jù)是1 km分辨率地理信息格網(wǎng)，以數(shù)組的形式存儲(chǔ)在地理定位產(chǎn)品MOD03文件中，其定位標(biāo)稱(chēng)精度為51 m。為驗(yàn)證本文算法，利用中國(guó)高差較大區(qū)域青藏高原地區(qū)MODIS影像做實(shí)驗(yàn)分析。本文實(shí)驗(yàn)算例A2020321.0535.061.2020321105426，實(shí)驗(yàn)環(huán)境如表1所示。使用高精度DEM數(shù)據(jù)可驗(yàn)證本文算法，但考慮高程和經(jīng)緯度的一致性，本文選取MODIS數(shù)據(jù)中DEM作為輸入DEM。MODIS數(shù)據(jù)提供1 000 m、500 m、250 m分辨率影像和1KM-DEM格網(wǎng)。因此，本文選取MOD03-1KM-DEM格網(wǎng)抽稀后的MOD03-4KM-DEM格網(wǎng)、MOD03-2KM-DEM格網(wǎng)作為控制格網(wǎng)，將余下的DEM格網(wǎng)點(diǎn)作為檢查點(diǎn)驗(yàn)證精度。

表1 實(shí)驗(yàn)硬件環(huán)境

2.2 精度評(píng)定與結(jié)論

為驗(yàn)證反算幀號(hào)準(zhǔn)確度，對(duì)GLT數(shù)據(jù)逐像元反算幀號(hào)并與其真實(shí)所在幀號(hào)對(duì)比。對(duì)GLT數(shù)據(jù)逐點(diǎn)測(cè)試本文算法，然后將結(jié)果與真實(shí)幀號(hào)相減。結(jié)果顯示誤差為零，這表明本文算法可由物方點(diǎn)精確定位對(duì)應(yīng)像方點(diǎn)所在幀。為驗(yàn)證子像素級(jí)精度，先由線(xiàn)性逼近反算得到未經(jīng)地形改正影像坐標(biāo)，再改正高差引起像點(diǎn)位移，將反算結(jié)果與原像點(diǎn)真實(shí)坐標(biāo)對(duì)比分析?？紤]蝴蝶結(jié)效應(yīng)引起多解，本文將距離掃描幀中心像點(diǎn)較近的點(diǎn)位作為輸出像方坐標(biāo)。

將本文算法與經(jīng)典的線(xiàn)性逼近和三次多項(xiàng)式插值進(jìn)行對(duì)比。如表2所示，在4∶1抽取實(shí)驗(yàn)中，在掃描向上，本文算法最大誤差為0.021個(gè)像元。線(xiàn)性逼近法和三次多項(xiàng)式插值法最大誤差分別為1.002和0.996個(gè)像元，相對(duì)于線(xiàn)性逼近和三次多項(xiàng)式插值，采樣誤差影響降低到可以忽略。如表3所示，在2∶1抽取實(shí)驗(yàn)中，在掃描向上，本文算法最大誤差為0.009個(gè)像元。線(xiàn)性逼近法和三次多項(xiàng)式插值法最大誤差分別為0.647和0.557個(gè)像元，相對(duì)于線(xiàn)性逼近和三次多項(xiàng)式插值在精度，采樣誤差影響降低到可以忽略。

表2 GLT 4∶1抽稀各個(gè)方法誤差對(duì)比

由于MODIS傳感器本身俯仰角較小且穩(wěn)定，在飛行方向地形起伏對(duì)像點(diǎn)位移的影響較小，因此本文算法和線(xiàn)性逼近法的反算精度都很高。如表2、表3所示，反算一次平均耗時(shí)不超過(guò)2.474 μs，這表明本文算法效率與傳統(tǒng)線(xiàn)性逼近法無(wú)顯著差異。

表3 GLT 2∶1抽稀各個(gè)方法誤差對(duì)比

在顧及地形改正的GLT數(shù)據(jù)的反算中，由于衛(wèi)星在飛行方向光線(xiàn)近似平行且俯仰角較小，精度雖有所提升但不顯著。在掃描方向上，衛(wèi)星成像下視角與掃鏡角度相關(guān)，易受地形起伏影響，成像時(shí)會(huì)產(chǎn)生大小不等的像點(diǎn)位移。如圖4所示，因此會(huì)對(duì)線(xiàn)性逼近法和三次多項(xiàng)式插值法的定位精度產(chǎn)生很大影響，最大誤差超過(guò)一個(gè)像元。而本文方法通過(guò)改正地形起伏引起像點(diǎn)位移，將采樣誤差影響降低到可以忽略。由圖4(b)可得，在削弱地形起伏影響后，因在掃描方向影像分辨率變化為非線(xiàn)性，隨著掃描角增大，線(xiàn)性方式產(chǎn)生的舍入誤差逐漸增大，因此精度逐漸下降。

圖4 抽稀比例4∶1時(shí)三種方法殘差分布和影像DEM圖

3 結(jié)束語(yǔ)

本文就顧及地形效應(yīng)的GLT數(shù)據(jù)反算難度大、精度低等問(wèn)題，針對(duì)受地形起伏影響大，提出一種帶有改正高差引起像點(diǎn)位移的線(xiàn)性逼近法，該方法具有精度高、速度快、受地形起伏影響較小等特點(diǎn)，在對(duì)GLT數(shù)據(jù)的坐標(biāo)反算中具有較大優(yōu)勢(shì)。利用青藏高原地形起伏較大的地區(qū)進(jìn)行實(shí)驗(yàn)，結(jié)果表明本文算法精度很高，優(yōu)于現(xiàn)有反算算法。