天宮二號衛星寬波段成像光譜儀地理定位方法

楊濤，何玉青，胡秀清

(1.北京理工大學 光電學院，北京 100081；2.國家衛星氣象中心，北京 100081)

0 引言

對地觀測衛星遙感圖像的地理定位是建立星載遙感儀器觀測圖像的每個像元與地基坐標系中位置關系的過程。傳統的遙感圖像地理定位方法分為參數法和非參數法[1]。非參數法主要是通過地面點與遙感數據來建立對應關系，因此選取合適的地面控制點是該方法的關鍵。但是對于運行中的遙感衛星來說，非參數法需要大量的地面點來進行輔助，并且在云覆蓋較多的情況下不易分辨地面點，因此不適合用于業務運行中的衛星。該方法主要用于早期星上測量手段缺乏或者星上測量數據丟失的情況，如QuickBrid、印度資源衛星、天繪一號等[2-4]。此外早期的星載高分辨率成像儀器常采用非參數定位方法，發展出了直接線性變換(direct lineaer transform，DLT)[5]、有理多項式系數(rational polynomial coefficients，RPC)[6]等一系列的模型。參數法是根據遙感儀器觀測幾何及其空間位置和指向，建立觀測像元與地面觀測位置之間的模型，其主要是利用精確測定的衛星位置、速度、姿態以及遙感儀器的掃描幾何和時序等參數來計算出每個像元的經度、緯度和觀測角度。目前絕大部分衛星都是使用參數法進行地理定位，如NOAA/AVHRR3[7]、Terra/MODIS[8]、Landsat-7/TM、NPP/VIIRS[9]，其不同之處主要體現在各個遙感儀器內部光學幾何結構以及各個衛星軌道計算方法。

天宮二號(tiangong No.2，TG-2)是我國首個真正意義上的空間實驗室，于2016年9月15日發射升空，其中包含了空間冷原子鐘、三維成像微波高度計和多角度寬波段成像光譜儀等14個精密載荷。寬波段成像光譜儀是其中的有效載荷之一，其主要用于觀測獲取海洋、大氣和陸地的各類參數，用于觀測海洋水色海溫以及優化天氣預報[10]。其對天氣預報進行優化主要是利用觀測數據可以改善氣溶膠反演數據，而這些科學和行業應用需要高質量的地理數據進行輔助，如經緯度和觀測角。因此對衛星定位數據的質量和精度提出了較高要求。

寬波段成像光譜儀采用的是多電荷耦合元件(charge-coupled device，CCD)拼接直視推掃的掃描方法。本文采用參數法的原理，針對其特殊的掃描特點，分析了寬波段成像光譜儀的工作原理與幾何模型，對整個儀器的地理定位方法進行分析與實驗，并對定位結果進行了初步評估。

1 寬波段成像光譜儀的工作原理與幾何模型

寬波段成像光譜儀根據工作波段分為三大模塊：可見光近紅外模塊、短波紅外模塊和熱紅外模塊。寬波段成像光譜儀總視場為42°，在考慮到可選用的探測器不能滿足條件的情況下，各模塊采用了單相機小視場合成總視場的方案：可見光近紅外采用3個子相機；短波和熱紅外模塊各采用2個子相機。同一波段的各個相機在同一平面內扇形展開，以可見光近紅外波段的相機為例，單個相機的視場略大于14°，3個相機兩兩成14°的夾角，合成一個42°的總視場。可見光近紅外波段相機的焦平面是由一個1 024×512(空間維×光譜維)面陣CCD構成，3個相機拼接成一幅完整圖像。因為相機的視場比夾角略大，因此3個相機成像時會出現重疊區域。TG-2軌道高度約為400 km，寬波段成像光譜儀可見光近紅外模塊的地面分辨率為100 m，瞬時視場為250 μrad。相機之間存在著32個暗電平。在去處暗電平之后，2個相機之間的重疊區域觀測角約為0.4°。根據瞬時視場計算可得重疊區域約為28個像元。

寬波段成像光譜儀安裝在衛星的對地面，開口朝向地球，3個模塊分層安置，同時進行推掃掃描，掃描方向就是衛星飛行方向，借助于衛星繞地球運行，獲取地球的二維景象。寬波段成像光譜儀可見光近紅外相機觀測地球景象的原理圖如圖1所示。

圖1 寬波段成像光譜儀可見光近紅外相機觀測地球景象原理圖

寬波段成像光譜儀單次采樣可以獲取4個數據：可見光近紅外相機獲取的15個通道的圖像數據、短波紅外相機獲取的2個通道的圖像數據、熱紅外相機獲取的2個通道的圖像數據以及相同時間段內平臺探測設備獲取的輔助參數數據。圖像信息數據包含的是各通道的圖像數據以及該數據起始時間、截止時間以及一個總計數時間。輔助參數數據包含定位所需的參數，即衛星軌道GPS位置與速度、衛星姿態角度信息、與GPS姿態數據相對應的協調世界時(UTC)以及其他平臺參數，如溫度、船下點、軌道高度等。

3個模塊的成像原理大致相同，不同之處只是相機數量、地面空間分辨率和各個相機之間安裝角度。短波和熱紅外模塊因為采用2個子相機，每個相機的光學視場為22°，兩兩夾角為21°。短波紅外相機的CCD為800×2，地面分辨率為200 m；熱紅外相機的CCD尺寸為400×2，地面分辨率為400 m。以下所述的定位方法均是以3個可見光近紅外波段相機為例，一個可見光近紅外數據即為可見光近紅外相機獲取的圖像數據以及與其匹配的平臺輔助參數數據。本次開機共獲取了600多個可見光近紅外數據，其中一半數據為海洋觀測區域，另一半數據為陸地區域，以中國為主。

2 寬波段成像光譜儀地理定位方法及流程

本節對寬波段成像光譜儀地理定位方法的整個流程都進行了詳細的設計，主要包括讀取原始數據并對數據做預處理、通過處理后的數據建立焦平面坐標系視向量與復合轉換矩陣、利用視向量計算地理定位信息、通過調整安裝矩陣的方法對定位結果進行了修正。定位方法具體流程圖如圖2所示。

圖2 地理定位方法流程圖

2.1 原始數據讀取與預處理

傳統的參數法定位大致需要2組參數，一組是與衛星相關的，主要是衛星軌道位置、速度以及運行時的姿態等；另一組與儀器相關，主要是儀器內部安裝幾何、探測器焦距像元尺寸等。其中儀器相關的數據在儀器研制過程中就已經進行了設計和測量，可以直接獲取。而衛星軌道相關數據，傳統上是通過衛星兩行報中的描述衛星運行的6個參數，利用衛星軌道模型來計算得到衛星實時位置信息。TG-2自身攜帶了GPS接收器以及星敏感器，可以實時測量得到衛星的位置以及姿態。因此對于該儀器的地理定位方法中不需要計算軌道參數，直接通過讀取L0級數據中的位置和姿態參數進行地理定位。

寬波段成像光譜儀的輔助信息和圖像信息分別存在2種不同格式的L0級數據中。因為2種數據的L0級文件格式不同，時間的采樣頻率與表示形式也均不相同，所以2種數據的時間不是一一對應的關系，而進行衛星定位時需要每幀圖像精確的觀測時間。所以需要進行預處理，即對輸出的數據時間進行統一，保證每一幀圖像數據都有輔助數據對應。

儀器主要的時間問題是由于圖像信息采樣頻率、輔助信息采樣頻率、平臺輔助參數的更新頻率以及時間表示不同造成的。其中平臺輔助參數的更新頻率最小，其時間是UTC時間。輔助信息采樣頻率是平臺輔助參數的更新頻率的一倍，可見光圖像信息采樣頻率又是輔助信息采樣頻率的40倍，二者的時間均是以開機時刻作為起始時間零的計數時間。

由于輔助數據包與圖像數據包采用的記錄時間不同，所以先要進行時間的統一，將輔助數據包采用的計數時間和圖像數據包采用的UTC時間都轉換為國際原子時。在統一時間單位后，通過時間來對輔助參數數據進行插值，實現輔助數據與圖像數據一一對應，插值的流程如下：

①利用圖像信息起始時間與截止時間獲得當前數據包的總觀測時間；

②根據采樣頻率比以及總時間得到每幀圖像的具體精確時間；

③利用每幀圖像的具體時間和平臺輔助信息的時間進行插值，獲得每幀圖像對應的具體輔助信息。

2.2 焦平面坐標系視向量及復合轉換矩陣構建

地理定位方法的基本原理就是通過儀器焦平面引出對地觀測的向量，通過該向量與地球模型的交點來得到XYZ坐標，再經過坐標系轉換最終得到經緯度等地理信息。首先要建立焦平面模型，來獲取成像光譜儀焦平面坐標系下的視向量。

雖然寬波段成像光譜儀的可見光波段由3個相機組成，但是在地理定位的過程中并不相互影響，因此本文把3個相機作為3個獨立的個體，分別建立焦平面模型。單個寬波段成像光譜儀的可見光相機是由一個1 024×512(空間維×光譜維)面陣CCD掃描成像，本文把掃描方向也就是飛行方向作為x軸正方向，穿軌方向為y軸方向，y軸正方向與圖像的數據存儲方向一致，如圖3所示。

圖3 成像光譜儀焦平面坐標系

由此可以得到焦平面上每個像元的x、y坐標以及焦距f，可以建立焦平面模型，也就是焦平面視向量ufoc。

(1)

焦平面視向量是根據焦平面坐標系建立，需要將其轉換到測地坐標系下，再結合地球模型才可以得到視向量與地球交點，利用交點的坐標信息通過坐標系轉換可以獲取所需的經緯度信息。其中涉及到望遠系統、儀器、衛星本體、軌道、地心慣性(earth centered inertial，ECI)、地心旋轉(earth centered rotate，ECR)和大地測量8個坐標系以及一系列的坐標轉換矩陣，如圖4所示，把這一系列坐標轉換矩陣組合可以得到一個復合轉換矩陣。

圖4 多個坐標系轉換圖

由焦平面坐標系轉換到望遠系統坐標系的轉換矩陣Ttel/foc主要與儀器的光學成像系統有關，一般由儀器研制方在發射前提前測量得到，通常使用單位矩陣代替。但是對于寬波段成像光譜儀來說，由于其是多CCD拼接成像，對應的焦平面有3個，這3個焦平面模型與望遠系統坐標系不盡相同。其中垂直掃描的相機可以使用單位矩陣，而兩側的相機與中間的相機成14°角，根據焦平面坐標系，其掃描方向為x方向，因此兩側的相機等于繞x軸旋轉正負14°，繞x軸旋轉的公式為：

(2)

式中：θ為旋轉的角度。由望遠系統坐標系轉換到儀器坐標系的轉換矩陣Tinst/tel主要與儀器掃描機理有關，該儀器成像方式為推掃成像，沒有掃描鏡調整方向，因此可以使用單位矩陣來表示這個旋轉矩陣。

(3)

儀器坐標系與衛星本體的坐標系之間的關系是由儀器的安裝矩陣確定，其是由儀器研制方在場地測試中計算得到的，用于調整儀器安裝時帶來的誤差。可以先將轉換矩陣Tsat/inst表示為單位矩陣，在衛星在軌運行之后需要繼續進行調整，以修正定位誤差[11]。

衛星的姿態決定了衛星坐標系與軌道坐標系之間的關系，該姿態由衛星上的星敏感器提供，表示形式為歐拉角，即滾動(roll)角、俯仰(pitch)角與偏航(yaw)角3個姿態角。其在直角坐標系下對應的就是繞X軸、Y軸與Z軸旋轉，而在衛星地理中一般按照ZXY的順序計算轉換矩陣[12]，因此該轉換矩陣Torb/sat可以表示為：

(4)

式中：y為偏航角；p為俯仰角；r為滾動角。

軌道坐標系與地心慣性系之間的關系是由衛星的瞬時位置與速度計算得到的，可以通過在ECI坐標系中形成軌道坐標系軸來構建從軌道轉換為ECI的旋轉矩陣Teci/orb，軌道坐標系的原點為衛星質心，Z軸方向對應的是衛星星下點方向，Y軸方向是衛星角動量方向的負方向，X軸方向由右手準則確定，即Z軸方向與Y軸方向的叉積。衛星星下點方向就是衛星在ECI坐標系下的坐標歸一化后的方向，而衛星角動量方向的負方向是垂直于衛星星下點方向與衛星速度方向的所構成的平面，公式如下：

(5)

式中：p表示衛星在ECI坐標系下的瞬時位置；v表示衛星在ECI坐標系下的瞬時速度。

從ECI坐標系到ECR坐標系的變換主要是由于時間不同、地球旋轉帶來的。此外還要考慮歲差、章動與極移這些變化的影響。這些矩陣可以通過國際天文聯合會提供的歲差章動極移模型得到。

將這些轉換矩陣構成一個復合轉換矩陣T：

T=Tecr/eciTeci/orbTorb/satTsat/instTinst/telTtel/foc

(6)

通過復合轉換矩陣與焦平面視向量，就可以得到在ECR坐標系下的視向量模型uecr

uecr=T×ufoc

(7)

2.3 獲取地理坐標

X=p+d·u

(8)

式中：p為衛星在ECR坐標系下的位置矢量；d為觀測點到衛星的距離，其是由通過地球橢圓模型重新縮放的視向量與衛星位置矢量構成的：

(9)

式中：u′與p′就是根據地球橢圓模型長短軸重新縮放的視向量與衛星位置矢量；a為長軸，b為短軸：

(10)

衛星運行觀察地球的幾何關系如圖5所示。

圖5 觀測交點幾何模型

最后將此位置矢量X(x1,x2,x3)換算為大地測量坐標(lat，lon，h)，此外地形和大氣折反射也會對遙感圖像地理定位精度造成影響，利用全球數字高程模型(digital elevation model，DEM)進行修正，本文地理定位采用來自美國國家影像與制圖局和世界數字化圖的30″分辨率DEM數據[13]。利用DEM找到對應區域高度信息，通過迭代的方法利用高度信息對ECR坐標系下視向量進行不斷修正直到高度信息相符，把修正后的視向量轉換到大地坐標系下得到修正后的經緯度和高度信息。

(11)

2.4 定位修正

因為3個相機的定位數據是單獨計算的，所以定位數據也會出現重疊區域。把3個相機獲取的定位數據拼接在一起時，需要考慮到重疊區域，拼接方法和圖像的拼接方法一致。2個相機之間有28個像元的重疊。考慮到邊緣像元分辨率下降的問題，于是在剔除重疊區域的時候，剔除各個相機更邊緣的區域，保留了靠近中心的區域。即兩邊的相機各自剔除14個像元，而中心相機兩邊各剔除14個像元。

利用拼接后的圖像與定位數據，主要選取具有特征的地面點，如海陸邊界線和湖泊，進行投影輸出。投影結果與高精度海陸分界線模板進行對比，可以發現存在定位結果誤差，效果如圖6所示。通過調整安裝矩陣可以減少誤差提高定位精度。安裝矩陣對圖像的修正效果體現在圖像上是衛星圖像繞坐標軸旋轉后的結果，坐標軸方向與焦平面坐標系的正方向一致，繞X軸旋轉會使圖像左右偏移，繞Y軸旋轉會使圖像前后偏移，繞Z軸旋轉會使圖像產生一定的旋轉偏移。相應的我們可以通過計算圖像中偏移像元數量來對這幾個偏移角度進行迭代調整，直到圖像定位誤差達到最小。

圖6 TG-2寬波段成像光譜儀未修正定位結果局部投影圖

可見光近紅外波段探測器的瞬時視場角ωIPOV為0.014 3°。本文挑選了50個帶有明顯地標特征的圖像數據進行統計，主要統計數據中心區域海岸線模板與圖像中海岸線的像素偏差。統計結果為在星下點附近左右方向平均偏移約13.41個像元，上下方向平均偏移約6.53個像元，整幅圖像兩側邊緣平均有個5.26像元的轉動，把像元偏移量帶入角度旋轉公式：

(12)

式中：Δny為左右平移量；Δnx為上下平移量；Δnz為旋轉平移量；L為每行的采樣點3 016個。最終計算得到繞X軸旋轉角φ為0.192°、繞Y軸旋轉角μ為0.093°、繞Z軸旋轉角ψ為0.099°。

把相應的角度帶入公式(4)可以計算出修正后的安裝矩陣，修正前后的定位改進效果可以參看圖7。

圖7 TG-2寬波段成像光譜儀局部放大圖像定位結果圖

在定位精度提高方面，除了對安裝矩陣進行調整的修正方法外，還可以采用基于視向量的修正方法[14]，即通過統計像元在掃面和穿軌方向上的誤差來建立像元與誤差的模型，利用模型修正視向量，從而減少定位誤差提高定位精度。本文也利用這一方法對第一種方法采用的50個數據進行了修正，其中，修正模型為視向量誤差隨影像列號成線性變化的一次線性模型[15]。通過分別統計2種方法的修正結果，我們發現修正安裝矩陣后平均誤差像元數為2.07，采用視向量修正模型后平均誤差像元數為1.99，進而每幅定位圖像與海陸模板的誤差仍處于2個像元左右，由此這2種方法的修正效果差異不明顯。與此同時，由于修正安裝矩陣的方法更容易實現，因而本文在后續的分析中仍然采用了對安裝矩陣進行調整的修正方法。

3 定位實驗與結果

3.1 實驗數據

寬波段光譜成像儀在TG-2入軌不久就開機獲取數據，此次開機運行持續了6個多月的時間。主要包含有2016年9月、10月、12月以及2017年1月、2月、3月的數據。12月份的數據最多每天都有數據，其他月份僅有十幾天有數據，其中每天的數據也不是整天所有軌道的數據都存在，以幾個或幾十個可見光近紅外數據為主。

在這600多個可見光近紅外數據的基礎上，本文做了進一步的篩選。因為數據中存在了整個圖像都是海洋的這種數據，這類數據對定位結果的判斷沒有幫助，因此剔除了這類數據。剩下的共301個可見光近紅外數據以中國陸地區域為主或者是包含海岸線區域，把這些數據作為有效數據。因為一個數據包含的圖像數據的垂直像元數量很多，大約有為50 000多個。為了顯示與存儲的方便，本文以4 800列為單位，將一個數據分解成了若干份。即一份實驗用的圖像數據為3 016×4 800，如圖8所示。

圖8 TG-2可見光近紅外相機原始圖像

3.2 實驗結果與誤差分析

本文從6個月的300多個有效可見光近紅外數據中每個月都隨機挑選了幾個可見光近紅外數據進行地理定位，利用修正后的定位結果與高精度海陸分界線模板比較，結果如圖7所示，圖中白色線條為高精度海陸分界線模板的投影。從局部放大的圖像中可以看到海陸分界線與圖像基本契合，二者之間存在幾個像元的誤差。

為了檢驗定位算法和修正矩陣的適用性，本文對所有的300多個有效可見光近紅外數據都進行地理定位，并且進行投影輸出。與平臺輔助信息提供的船下點經緯度比較，誤差在0.01°以內。本文把圖7中對應數據星下點的經緯度與平臺輻射信息提供的經緯度進行了對比，如表1所示。

此外本文還把海陸模板疊加在投影圖像上，統計了二者相差的像元數量，考慮到海陸模板的精度問題，每幅定位圖像與海陸模板的誤差在2個像元左右，如表2所示。此外本文對寬波段光譜成像儀的其他2個波段的圖像也進行了地理定位，由于短波紅外波段相機與熱紅外波段的相機空間分辨率都低于可見光波段的相機，定位精度有所提升，熱紅外波段相機誤差在1個像元左右。

表1 實驗定位信息與平臺輔助信息比較表

表2 定位結果誤差統計表

定位結果中主要誤差是由衛星觀測時間的偏差、衛星姿態測量精度、衛星位置測量精度、發射后儀器安裝偏差等造成的。其中安裝偏差是誤差的主要來源，其是由于安裝誤差和衛星發射過程中的劇烈震動，以及發射后的環境的影響，導致儀器與衛星之間的位置關系發生改變，即使采用地面測試得到安裝矩陣依舊無法完全消除定位誤差，需要通過對圖像進行統計分析，調整安裝矩陣來進行修正。另外，衛星姿態測量和衛星位置測量精度還存在隨機誤差，它會隨著時間的變化而變化，造成圖像出現一定的平移與旋轉，這種變化不易進行模擬。

4 結束語

天宮二號寬波段成像光譜儀為多CCD拼接直視推掃模式。本文詳細設計了基于該儀器掃描模式的參數型定位方法，其中包括數據讀取與預處理、建立焦平面視向量與坐標系轉換矩陣、地理參數計算以及定位修正這幾個過程。本文通過衛星的空間位置姿態參數，建立了儀器焦平面與地理空間位置的對應模型。通過大量實際觀測的定位結果建立旋轉角度與變形數量的關系，修正安裝矩陣，提高模型定位精度。利用該模型實現了對TG-2衛星寬波段成像光譜儀圖像的地理定位。利用圖像的地理位置進行投影后，與高精度海陸模板進行比較，定位精度約為2個像元。該研究成果為后期的TG-2寬波段成像光譜儀的圖像應用，以及海洋、大氣和陸地地球物理參數的定量反演奠定了堅實的基礎。