一種線陣推掃式行星遙感影像攝影測量快速處理方法

耿 迅，徐 青

（1.河南大學 地理與環(huán)境學院，開封 475004；2.河南大學 河南省時空大數(shù)據(jù)產業(yè)技術研究院，鄭州 450000；3.信息工程大學 地理空間信息學院，鄭州 450052）

引 言

地外天體表面高精度的正射影像圖（Digital Orthophoto Map，DOM）與數(shù)字高程模型（Digital Elevation Model，DEM）是深空探測任務實施與行星科學研究的重要基礎數(shù)據(jù)[1-4]。使用攝影測量方法對地外天體遙感影像進行幾何處理是制作行星制圖產品的主要技術手段[5-7]?，F(xiàn)有行星遙感制圖產品的精度與分辨率仍然較低，不能很好地滿足深空探測工程任務以及深層次行星科學研究的需要。當前，利用不同航天機構獲取的多探測任務數(shù)據(jù)制作高分辨率、高精度的行星制圖產品是深空探測遙感測繪的研究重點，但是融合多探測任務數(shù)據(jù)進行攝影測量處理的技術難度很大，而且行星攝影測量處理過程復雜，可用軟件資源有限。

結合地外天體遙感影像的特點，有針對性地設計、優(yōu)化行星攝影測量處理方法，是深空探測任務科學數(shù)據(jù)處理與應用的重要工作。深空探測器通常搭載線陣傳感器獲取地外天體的遙感影像用于制作行星制圖產品[8-11]。然而，現(xiàn)有行星攝影測量處理方法在用于線陣推式遙感影像時的效果不佳，通常表現(xiàn)出較低的處理效率[12]，美國地質調查局（United States Geological Survey，USGS）開發(fā)的ISIS（Integrated System for Imagers and Spectrometers）軟件在處理生成大數(shù)據(jù)量（例如GB級）線陣推掃式行星遙感影像的正射影像圖時，通常需耗時幾十分鐘甚至若干小時[13]。另外，在構建連接點控制網(wǎng)時，現(xiàn)有行星攝影測量處理方法針對線陣推掃式遙感影像的適用性有限，經常出現(xiàn)匹配失效的情況，導致后序攝影測量處理步驟以及行星制圖產品精度受到影響[14]。針對現(xiàn)有行星攝影測量處理方法在線陣推掃式遙感影像上的不足，本文采用快速反投影算法進行正射影像糾正，并且在近似正射影像空間匹配獲取連接點，提升長條帶、大數(shù)據(jù)量線陣推掃式行星遙感影像的攝影測量處理效率。

1 行星攝影測量處理技術

經過幾十年的積累，美國在深空探測遙感測繪領域已經形成了一套比較成熟的行星攝影測量處理技術體系，可以將這個技術體系概括為“PDS+SPICE+ISIS”，其中PDS（Planetary Data System）是行星遙感數(shù)據(jù)系統(tǒng)，SPICE是深空探測任務輔助信息系統(tǒng)，ISIS是開源行星遙感影像處理軟件。另外，近年來美國國家航空航天局（National Aeronautics and Space Administration，NASA）研發(fā)的ASP（Ames Stereo Pipeline）開源軟件也逐漸在行星制圖領域中得到應用。

PDS是NASA研發(fā)的行星遙感數(shù)據(jù)系統(tǒng)，目前作為深空探測任務科學數(shù)據(jù)的標準文件格式被世界各國廣泛采用。中國“嫦娥”工程與“天問一號”火星探測工程的科學數(shù)據(jù)也采用了PDS文件格式。PDS文件除了記錄影像信息，也包含有影像獲取時間、分辨率、地理坐標范圍等輔助信息。可以使用開源圖像處理引擎GDAL對PDS文件進行讀寫、格式轉換等操作。

行星遙感影像幾何處理涉及到復雜的行星坐標系轉換、各種時間系統(tǒng)轉換、探測器位置與姿態(tài)計算、相機模型構建、光束法平差等內容。NASA使用SPICE庫文件對這些深空探測任務輔助數(shù)據(jù)進行歸檔[15]。SPICE庫文件也是行星攝影測量處理時重要的輔助數(shù)據(jù)，其中SPK文件對應于探測器位置數(shù)據(jù)，CK文件對應于探測器姿態(tài)數(shù)據(jù)，利用這兩個文件可以生成影像的外方位元素信息，用于構建相機模型以及光束法平差處理。另外，為便于利用這些SPICE庫文件進行各種計算，SPICE輔助信息系統(tǒng)也提供了相應的庫函數(shù)，可以使用C語言進行調用。

美國自實施“阿波羅”登月計劃以來，一直十分重視深空探測遙感測繪工作，基于USGS ISIS軟件制作了月球、火星、小行星等多個星體的制圖產品，在深空探測著陸區(qū)選址、任務路徑規(guī)劃以及行星科學研究等方面發(fā)揮了重要作用。USGS開發(fā)的ISIS行星攝影測量軟件系統(tǒng)能夠支持NASA的多個深空探測任務數(shù)據(jù)。ISIS提供了很多工具軟件用于行星遙感影像的攝影測量處理，其它航天機構，如歐洲航天局（European Space Agency，ESA）、日本宇宙航空研究開發(fā)機構（Japan Aerospace Exploration Agency，JAXA）也使用ISIS進行深空探測科學數(shù)據(jù)的預處理，ISIS在深空探測遙感測繪領域得到了行業(yè)的廣泛認可。

ASP是由NASA的艾姆斯研究中心（Ames Research Center）開發(fā)的立體攝影測量軟件，可以支持行星遙感影像與對地觀測遙感影像[16]。ASP軟件在生成DEM方面的效果較好，并且自動化程度較高，提供有stereo、point2dem等多個工具軟件用于自動生成DEM，stereo軟件主要是匹配立體影像生成點云文件，point2dem可以將點云文件轉換為常用格式的DEM文件（如GeoTiff格式），并且生成DEM交會殘差圖輔助分析DEM成果精度。ASP支持ISIS的cube文件格式，能夠讀取其相機模型信息。

但是與商業(yè)攝影測量軟件相比，現(xiàn)有開源行星攝影測量軟件在交互式編輯、線陣推掃式遙感影像處理效率等方面的功能稍弱。然而，商業(yè)攝影測量軟件通常無法直接用于處理月球、火星等地外天體遙感影像，主要原因是商業(yè)攝影測量軟件無法支持行星遙感影像的數(shù)據(jù)格式、相機模型、坐標基準等；另外，商業(yè)攝影測量軟件的核心算法（例如影像匹配）也通常結合對地觀測衛(wèi)星影像的特點進行了優(yōu)化，而這些優(yōu)化參數(shù)并不適合于處理行星遙感影像。

2 基于快速幾何糾正的行星攝影測量處理方法

本文基于線陣推掃式遙感影像快速反投影算法進行行星遙感影像的幾何糾正與控制網(wǎng)構建，并自主研發(fā)了相應的軟件模塊，應用于地外天體的攝影測量處理過程?；诳焖賻缀渭m正的線陣行星遙感影像攝影測量處理流程如圖1所示，具體處理步驟可以分為影像預處理、相機模型構建、連接點提取、光束法平差以及制圖產品生成等幾個步驟。經過幾何處理生成的正射影像圖與數(shù)字高程模型可以加載至ArcGIS、GlobalMapper、QGIS等地理信息系統(tǒng)（Geographic Information System，GIS）軟件中進行科學分析。由于USGS ISIS支持多個深空探測任務的科學數(shù)據(jù)，而且在行星測繪領域得到廣泛應用，尤其是ISIS內置了許多深空探測任務的相機參數(shù)，因此本文的行星攝影測量處理流程也使用USGS ISIS軟件進行格式轉換、輻射校正以及光束法平差等處理。

圖1 基于快速幾何糾正的行星遙感影像攝影測量處理流程圖Fig.1 Flowchart of the photogrammetric processing of planetary remote sensing images based on fast geometric rectification

在影像預處理階段需要導入PDS格式的原始影像并轉換為USGS ISIS系統(tǒng)支持的cube文件格式。ISIS針對不同深空探測任務的傳感器開發(fā)了專門的數(shù)據(jù)導入接口，例如lronac2isis用于導入月球偵察軌道器（Lunar Reconnaissance Orbiter，LRO）窄角相機（Narrow Angle Camera，NAC）影像，hrsc2isis用于導入火星快車（Mars Express，MEX）高分辨率立體相機（Hight Resolution Stereo Camera，HRSC）影像。在預處理階段還需要對影像進行輻射檢校，例如使用lronaccal與lronacecho對LRO NAC影像進行輻射檢校。在預處理階段的一個重要步驟是調用spiceinit獲取與影像相關的各個SPICE庫文件，從而計算出影像在攝影時刻的位置、姿態(tài)信息，這一步相當于構建了相機模型，以便在像點與地面點之間進行坐標轉換。

2.1 線陣推掃式影像地面點反投影

構建遙感影像的相機模型是開展攝影測量處理的基礎[17-18]。相機模型主要用于地面點坐標P(X,Y,Z)與像點坐標p(x,y)的相互轉換。以MEX HRSC影像為例，地面點P(X,Y,Z)與像點p(x,y)的坐標轉換公式為

圖2 線陣推掃式影像地面點反投影示意圖Fig.2 Illustration of back projection of linear pushbroom images

為了提升最佳掃描線的迭代計算效率，本文采用基于物方投影面約束的快速反投影算法進行線陣推掃式影像的地面點反投影計算[18-19]。先利用線陣推掃式影像的嚴密幾何模型建立每一條掃描線的物方投影面（例如第i條掃描線的物方投影面為SiAiBi），然后在地面點反投影迭代計算時，利用地面點至投影面的距離作為幾何約束條件輔助確定地面點對應的最佳掃描線。地面點P與某一掃描線所構成的物方投影面之間的距離 L的計算公式為

其中：A、B、C、D是物方投影面的法線方向，可由空間3點Si、Ai、Bi的三維坐標計算得出。由于基于物方投影面的迭代計算過程是簡單的幾何計算，并且迭代收斂速度快，因此地面點反投影過程的計算效率可以得到顯著提升。

2.2 快速幾何糾正

線陣推掃式遙感影像的幾何糾正通常采用反解法，即對于正射影像圖上的每一個像素，先計算出其對應的地面點坐標，再由地面點坐標反求出像點坐標（即地面點反投影過程），然后將該像點的灰度值進行內插處理并賦值給正射影像圖。USGS ISIS系統(tǒng)提供了cam2map工具軟件用于生成DOM。由于ISIS系統(tǒng)的線陣推掃式相機模型采用計算效率偏低的二分法確定最佳掃描線，導致地面點反投影計算的效率較低，因此cam2map在處理線陣推掃式行星遙感影像（例如LRO NAC、MEX HRSC）時非常耗時。

本文通過優(yōu)化線陣推掃式遙感影像的相機模型提升地面點反投影迭代計算效率，并且采用多線程程序設計方法進一步提高正射影像生成的處理效率，與USGS ISIS軟件相比，單線程幾何糾正處理的計算效率可提升3～5倍，多線程幾何糾正處理的計算效率可提升10～20倍（4線程）。

2.3 連接點控制網(wǎng)構建

通過spiceinit獲取的影像初始位置、姿態(tài)信息（即外方位元素）的精度偏低，需要利用光束法平差進行精化。光束法平差的關鍵是獲取均勻分布的可靠連接點，在USGS ISIS系統(tǒng)中稱為控制網(wǎng)（control network）。高精度的連接點控制網(wǎng)是獲得高質量行星攝影測量成果的前提。USGS ISIS提供了autoseed、cnetref、pointreg等一系列控制網(wǎng)構建工具，但是其計算效率較低，并且經常出現(xiàn)匹配處理失敗的情況，不適合大數(shù)據(jù)量的連接點提取，而且相關匹配算法難以用于畸變較為嚴重的立體影像（例如MEX HRSC的前、后視圖像）。

本文先利用線陣推掃式遙感影像的快速反投影算法生成近似正射影像，然后在近似正射影像上采用歸一化相關系數(shù)方法匹配連接點，再利用地面點反投影算法將匹配的連接點轉換至原始影像空間。由于在近似正射影像空間進行匹配，立體影像的分辨率變?yōu)橐恢?，同時也消除了地形起伏與攝影角度不同帶來的幾何畸變影響，并且能夠利用影像初始位置、姿態(tài)信息輔助確定匹配點的起始位置，因此能夠提升連接點提取步驟的計算效率與匹配精度。

2.4 光束法平差

USGS ISIS提供了jigsaw工具用于光束法平差，jigsaw軟件的計算效率較高，對深空探測任務各種傳感器的支持也比較好。因此，本文將匹配的連接點轉換至ISIS軟件支持的PVL （Parameter Value Language）控制網(wǎng)文件格式進行光束法平差處理。光束法平差的一個難點是設置不同類型觀測值的權值，主要涉及探測器位置與姿態(tài)觀測值、控制點觀測值等。權值設置可以依據(jù)探測器定軌、定姿精度進行設置，也需要一定的工程經驗。光束法平差通常需要經過幾次迭代，每一次平差時需要剔除連接點中的粗差，通常認為當像點觀測值殘差小于1個像素，平差指標Sigma0值小于0.5時可以認為獲得了理想的平差結果（實際精度情況受探測數(shù)據(jù)質量影響）。若要提升絕對幾何定位精度，則需要在光束法平差過程中引入控制數(shù)據(jù)。但是，現(xiàn)有行星表面控制數(shù)據(jù)的分辨率與精度較低，尤其是缺乏高精度的平面控制數(shù)據(jù)，實際工程應用中通常引入激光測高數(shù)據(jù)或者由激光測高數(shù)據(jù)生成的DEM作為高程控制信息。

2.5 數(shù)字高程模型與正射影像圖生成

在光束法平差后，利用立體像對進行影像匹配可以獲取密集的匹配點，再利用立體影像空間前方交會計算出點云的三維坐標，從而自動構建出DEM。通常使用ISIS軟件進行行星遙感影像的預處理與光束法平差，再利用ASP軟件生成DEM。ASP軟件自動生成DEM的分辨率與原始影像分辨率基本相同，但是存在較多噪聲，實際應用時通常需要將自動生成的DEM重采樣為3～5倍原始影像分辨率。

通過幾何糾正可以消除各種畸變的影響，生成含有精確地理信息的DOM。生成DOM時需要利用影像的外方位元素以及參考DEM，因此光束法平差的精度以及DEM的精度會直接影響DOM的成果精度。在生成正射影像圖時，可以利用USGS ISIS提供的cam2map軟件以及自主研發(fā)的快速幾何糾正軟件。

3 試驗與分析

試驗選取4幅LRO NAC影像與5幅MEX HRSC影像進行處理（影像基本信息見表1）。LRO NAC影像數(shù)據(jù)位于“嫦娥四號”著陸區(qū)，影像分辨率為0.88～1.18 m；MEX HRSC影像數(shù)據(jù)覆蓋Mars 2020任務預選著陸區(qū)的Jezero撞擊坑，影像分辨率為14～22 m。利用USGS ISIS 3.5.2軟件進行格式轉換、輻射校正、光束法平差等處理，使用自主研發(fā)軟件進行相機模型構建、幾何糾正、連接點提取等處理，使用ASP 2.6.0軟件生成數(shù)字高程模型。試驗軟件環(huán)境為64位Ubuntu 14.04操作系統(tǒng)（安裝在VMWare虛擬機上）；硬件環(huán)境為Intel Core i7-7 500 2.7 GHz CPU與8 GB內存。月球與火星制圖分別采用1 737.4 km與3 396.19 km的正球體，地圖投影采用Equirectangular投影方式。

表1 測試影像基本信息Table 1 Basic information of the test images

3.1 光束法平差精度分析

基于初始的影像位置、姿態(tài)信息以及USGS ISIS軟件內置的月球、火星全球DEM數(shù)據(jù)，使用自主研發(fā)的幾何糾正軟件模塊將試驗影像處理生成近似正射影像圖，在近似正射影像圖上提取連接點，再轉換為ISIS系統(tǒng)支持的PVL控制網(wǎng)文件格式用于光束法平差。MEX HRSC影像共提取495個連接點，LRO NAC測試影像共獲取187個連接點，兩組實驗數(shù)據(jù)的連接點分布情況見圖3～4。

圖3 MEX HRSC影像連接點分布示意圖Fig.3 Illustration of tie points distribution of MEX HRSC images

圖4 LRO NAC影像連接點分布示意圖Fig.4 Illustration of tie points distribution of LRO NAC images

使用jigsaw進行光束法平差，兩組實驗數(shù)據(jù)的Simga0值均小于0.3，光束法平差后影像像方殘差值見表2，殘差分布見圖5～6。分析像方殘差結果可知，LRO NAC與MEX HRSC影像光束法平差的像方殘差值均小于半個像素，表明影像匹配精度較高，不同條帶之間外方位元素的不一致性得到消除。連接點幾何定位精度見圖7～8。分析結果可知，LRO NAC影像平面定位精度為0.5～1.0 m，高程定位精度為1.0～1.5 m；MEX HRSC影像平面定位精度為4～10 m，高程方向上大部分連接點的定位精度為10～20 m，少量連接點的定位精度接近30 m。結合影像分辨率可知，LRO NAC與MEX HRSC影像均可達到平面約0.5～1個像素，高程約1～2個像素的定位精度，與對地觀測衛(wèi)星遙感影像的幾何定位精度基本一致。

圖5 LRO NAC影像光束法平差后像方殘差分布圖Fig.5 Image residuals of the bundle adjustment for LRO NAC images

圖6 MEX HRSC影像影像光束法平差后像方殘差分布圖Fig.6 Image residuals of the bundle adjustment for MEX HRSC images

圖7 LRO NAC影像光束法平差連接點幾何定位精度Fig.7 The posterior geometric accuracy of tie points for LRO NAC images

圖8 MEX HRSC影像光束法平差連接點幾何定位精度Fig.8 The posterior geometric accuracy of tie points for MEX HRSC images

表2 光束法平差影像殘差值Table 2 Image residuals of bundle adjustment results

3.2 數(shù)字高程模型精度分析

光束法平差之后，可以利用ASP的stereo、point2dem等工具軟件自動生成DEM。LRO NAC原始影像之間幾何變形較小，可以直接用于密集匹配；而MEX HRSC影像由于攝影角度差異大，影像之間幾何變形較大，直接在原始影像上匹配比較困難，因此可以先糾正為正射影像，然后在正射影像上匹配生成DEM。LRO NAC測試影像選用M1303619844LE與M1303640934LE組成立體像對，而MEX HRSC測試影像選用h7289 nd2、s12、s22（即下視、前視、后視）組成立體像對。自動生成的DEM與交會殘差見圖9～10，LRO NAC立體影像生成的DEM重采樣為5 m格網(wǎng)間距，MEX HRSC立體影像生成的DEM重采樣為50 m格網(wǎng)間距。分析LRO NAC與MEX HRSC立體影像的DEM交會殘差圖可知，大部分區(qū)域的交會殘差值為1～2個像素，表明影像匹配精度與光束法平差精度較高。受光照不足、陰影、紋理稀疏等因素影響，立體影像中也有少部分區(qū)域未匹配出同名點，例如圖9（a）左下角的撞擊坑（圖中顯示為白色無效值區(qū)域）。另外，圖9 （b）中DEM交會殘差圖顯示出有規(guī)律的條紋形狀，可能是由衛(wèi)星平臺的高頻顫振引起，本文未對顫振做進一步處理。

圖9 LRO NAC立體影像自動提取DEM結果Fig.9 The DEM results of LRO NAC stereo images

圖10 MEX HRSC立體影像自動提取DEM結果Fig.10 The DEM results of MEX HRSC stereo images

由于月球、火星表面通常缺乏高精度的絕對控制數(shù)據(jù)，行星攝影測量的絕對定位精度通常依賴于探測器定軌、定姿數(shù)據(jù)的精度。利用行星遙感影像生成DEM時的交會殘差及連接點幾何定位精度能夠表明立體影像之間的內部符合精度較好，實際上是對光束法平差后相對定位精度的評價。為了提升絕對定位精度，可以使用激光測高數(shù)據(jù)獲取的DEM添加高程控制，或者將激光測高足印數(shù)據(jù)與影像進行融合處理，即將高精度的激光測距信息作為觀測值與影像觀測值進行聯(lián)合平差，進一步提升絕對定位精度。

3.3 正射影像圖精度分析

利用光束法平差后的精化外方位元素以及自動生成的高分辨率DEM糾正生成DOM，結果如圖11～12所示。圖11（a）是M1303619844LE/RE構成的正射影像鑲嵌圖，影像分辨率為0.9 m，圖中藍色方框標注了圖11（b）、（c）局部區(qū)域影像拼接效果（M1303619844LE與M1303640934LE重疊區(qū)域），紅色五角星標注了“嫦娥四號”著陸點以及“玉兔二號”月球車在影像獲取時刻的位置，具體可參見圖11（d）。圖12（a）是利用h7289 nd2影像制作的正射影像圖，影像分辨率為14 m，影像中間是Jezero撞擊坑，圖中藍色方框標注了右側圖12（b）、（c）局部區(qū)域影像拼接效果（h7289 nd2與s12重疊區(qū)域），左側為nd2通道影像，右側為s12通道影像。由正射影像拼圖結果可知，使用光束法平差精化的外方位元素以及生成的高分辨率DEM可以顯著降低正射影像圖的幾何拼接誤差。

圖11 “嫦娥四號”著陸區(qū)LRO NAC正射影像圖Fig.11 The digital orthophoto map of LRO NAC images for the landing site of Chang’E-4

圖12 Mars 2020預選著陸區(qū)MEX HRSC正射影像圖Fig.12 The digital orthophoto map of MEX HRSC images for the candidate landing site of Mars 2020

4 結 論

本文針對線陣推掃式行星遙感影像攝影測量處理過程中的正射影像生成與控制網(wǎng)構建步驟進行算法優(yōu)化，自主研發(fā)了快速幾何糾正與連接點提取軟件模塊，顯著提升了線陣推掃式行星遙感影像的攝影測量處理效率。但是，自主研發(fā)的軟件模塊目前僅能支持LRO NAC與MEX HRSC的相機模型，下一步需要結合月球、火星全球遙感制圖的需求，支持更多的深空探測任務數(shù)據(jù)，并采用GPU計算、集群處理等技術進一步優(yōu)化算法性能，以適用于大數(shù)據(jù)量處理的需求。

中國的深空探測將持續(xù)向著“開放、合作、共建、共享”的方向發(fā)展，而科學數(shù)據(jù)的高效處理是深空探測工程任務與行星科學研究的基礎。當前，綜合利用國內外多探測任務數(shù)據(jù)制作更高精度、更高分辨率的行星制圖產品仍然面臨諸多技術難題。針對行星遙感影像特點，自主研發(fā)一系列的攝影測量軟件工具是深空探測遙感測繪中必不可少的基礎性工作。