顧及地形改正的地理定位查找表模型反算

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

(中南大學 地球科學與信息物理學院，長沙 410083)

0 引言

GLT數據是常見于擺掃式衛星影像中的定位格網數據。基于GLT數據的坐標反算是衛星數據幾何校正中的重要環節，尋求高質量坐標反算方法對擺掃式衛星數據處理至關重要。由于GLT模型是由影像坐標獲取地理位置的正向索引，沒有確定的反算函數，因此基于GLT數據的坐標反算較為困難。

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

針對上述問題，本文提出了一種顧及地形效應的線性逼近高精度坐標反算方法。先用線性逼近法將輸入物方點反算至未經地形改正像方點，然后改正由高差引起的像點位移實現高精度坐標反算。該方法能夠克服目前存在的受地形起伏影響大、精度低的缺點。

1 算法原理

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

圖1 坐標反算流程

1.1 改正高差引起像點位移

星載擺掃式成像系統是通過掃描方向的擺掃及沿軌道方向上運動實現二維影像的獲取[10]，因此影像的成像視角在掃描向變化較大。為保證平臺穩定性，衛星姿態保持本體坐標系與軌道坐標系一致，因此可認為影像在飛行方向上成像光線近似平行。地形起伏引起的像點位移主要在掃描向，飛行方向上位移相對較小。

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

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

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

圖2 地形起伏引起在線性內插中產生像點位移示意圖

(1)

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

(2)

式中：RN為地球半徑；H為衛星高。

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

tan2α=tan2αx+tan2αy

(3)

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

αx≈αm+ω

(4)

掃描鏡掃描角αm由線性逼近所求未經地形改正的影像坐標(xu，yu)計算求得，表達如式(5)所示。

(5)

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

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

(6)

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

(7)

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

1.2 線性逼近法坐標反算

基于單幀影像地理信息數據進行坐標反算過程如圖3所示。線性逼近坐標反算思想是根據GLT數據中像方點的二維梯度和輸入物方經緯度，任意給定初始像點坐標，計算像方改正量，迭代計算像點坐標。當二元函數模型λ(x，y)、φ(x，y)并未直接給出時，很難根據(λ，φ，h)求出像方坐標(x，y)。考慮到幀內連續變化，可根據所給經緯度二維數據進行線性迭代逼近。將λ、φ泰勒展開得式(8)。

(8)

式中：φ′y、φ′x、λ′y、λ′x為該點的一階偏導數。可以根據該像點鄰近像點經緯度數據計算梯度，解得如式(9)所示。

(9)

式中：D=λ′xφ′y-φ′xλ′y。預設給定任意初始像點坐標可由(x0，y0)，再由所給地理坐標(λ，φ，h)，以及初始點二維梯度代入式(9)算出對應像方改正量(Δx，Δy)。將(x0，y0)+(Δx，Δy)作為新計算的像點坐標，并繼續迭代，直到兩次迭代計算出的坐標在一個GLT網格內結束迭代，再由該GLT網格經緯度梯度和所給物方坐標線性內插出未經地形改正的像方坐標(xu，yu)。

圖3 迭代逼近示意圖

1.3 反算像方點幀號

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

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

步驟2：將物方坐標(λ，φ，h)與代入步驟1中數組，經線性逼近反算得像方坐標。

步驟3：抽取距離步驟2中像點坐標最近的兩幀數據，分別由線性逼近反算影像坐標，計算像點與各幀中心距離，取距離較短的幀號作為該物方點對應像方點所在幀號。

1.4 特殊區域處理

當影像跨越子午線時，影像覆蓋區域的經度會從180°變成-180°，導致GLT數據出現跳變，若不特殊處理則會產生粗差，或迭代不收斂。本文算法先判斷該數據是否包含子午線部分，若包含則將經度全部加上360°以保證經度連續。經過該調整后對反算效果不產生影響。

2 實驗設計與分析

2.1 實驗設計

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

表1 實驗硬件環境

2.2 精度評定與結論

為驗證反算幀號準確度，對GLT數據逐像元反算幀號并與其真實所在幀號對比。對GLT數據逐點測試本文算法，然后將結果與真實幀號相減。結果顯示誤差為零，這表明本文算法可由物方點精確定位對應像方點所在幀。為驗證子像素級精度，先由線性逼近反算得到未經地形改正影像坐標，再改正高差引起像點位移，將反算結果與原像點真實坐標對比分析。考慮蝴蝶結效應引起多解，本文將距離掃描幀中心像點較近的點位作為輸出像方坐標。

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

表2 GLT 4∶1抽稀各個方法誤差對比

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

表3 GLT 2∶1抽稀各個方法誤差對比

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

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

3 結束語

本文就顧及地形效應的GLT數據反算難度大、精度低等問題，針對受地形起伏影響大，提出一種帶有改正高差引起像點位移的線性逼近法，該方法具有精度高、速度快、受地形起伏影響較小等特點，在對GLT數據的坐標反算中具有較大優勢。利用青藏高原地形起伏較大的地區進行實驗，結果表明本文算法精度很高，優于現有反算算法。