基于卡爾曼濾波的Landsat-8衛星影像幾何精校正

基于卡爾曼濾波的Landsat-8衛星影像幾何精校正

郭少鋒1，2，李安1，李山山1，馮鐘葵1

(1.中國科學院遙感與數字地球研究所，北京 100094；2.中國科學院大學，北京 100049)

摘要：卡爾曼濾波是一種遞推線性最小方差估計算法，被廣泛應用于GPS動態數據處理、組合導航和通信與信號處理等領域。本文將卡爾曼濾波引入衛星影像幾何精校正處理中，以Landsat-8衛星影像為例，利用卡爾曼濾波方法和地面控制點構建運動方程和測量方程，實現了對衛星星歷和姿態數據的修正。實驗證明，相比于傳統最小二乘平差法，該方法使用較少的控制點即可達到穩定的精度，降低了幾何精校正對控制點數量的依賴，具有較高適用性。

關鍵詞：Landsat-8；卡爾曼濾波；幾何精校正

doi:10.3969/j.issn.1000-3177.2015.01.003

中圖分類號：TP751文獻標識碼：A

收稿日期：2013-11-28修訂日期：2014-03-06

基金項目：江西省數字國土重點實驗室開放基金(DLLJ201303)；航空遙感技術國家測繪地理信息局重點實驗室經費資助項目(2014B02)。

作者簡介：于英(1985～)，男，博士，主要研究基于移動平臺的視頻遙感圖像快速處理方法。

Geometric Precision Correction of Landsat-8 Imagery Based on Kalman Filter

GUO Shao-feng1，2，LI An1，LI Shan-shan1，FENG Zhong-kui1

(1.InstituteofRemoteSensingandDigitalEarth，ChineseAcademyofSciences，Beijing100094；

2.UniversityofChineseAcademyofSciences，Beijing100049)

Abstract：Kalman filter is a recursive linear minimum variance estimation algorithm，and it has been widely used in GPS dynamic data processing，integrated navigation，communication and signal processing fields.This paper applies kalman filter into geometric precision correction of Landsat-8 imagery，and completes the correction of satellite ephemeris and attitude data.In the experiments，the comparison of accuracy between kalman filter algorithm and least square method shows that the application of kalman filter in precision correction of landsat-8 imagery can reduce the dependence of GCPs.

Key words：Landsat-8；kalman filter；geometric precision correction

1引言

隨著航天技術、傳感器技術和計算機技術的蓬勃發展，遙感技術在近30年取得了很大的進步，并被越來越廣泛地應用于國民經濟建設中。為了適應多層次、多角度、全方位和全天候的全球立體對地觀測的要求，近兩年世界遙感強國相繼發射了眾多遙感衛星。其中，美國為了延續Landsat系列衛星的對地觀測任務，于2013年2月11日成功發射Landsat-8衛星。該衛星分辨率為30m，衛星重訪周期為16天，帶有OLI和TRIS兩個對地觀測傳感器。地物變化監測成為其最重要的研究目標之一，這對衛星的定位精度提出了較高的要求，因此針對Landsat-8圖像進行幾何精校正處理十分重要。

幾何精校正主要通過地面控制點來建立遙感圖像與地面之間的精確關系，利用數學平差方法，校正由于衛星星歷數據、姿態數據等測量誤差造成的衛星圖像的幾何誤差[1]。Yoichi提出了一種雙向掃描型傳感器幾何校正的方法，分別用高低頻分離方法建立了正掃和反掃視點模型以此校正幾何誤差[2]；Kratky等對共線方程加入軌道約束條件，提出了基于軌道約束的嚴格成像幾何模型[3]；Gruen等從航測角度出發，對SPOT影像應用共線方程，并開發出光束法平差模式求解衛星外方位元素，其校正精度可以達到亞像元級[4]；張過等人將有理函數模型應用于推掃式光學衛星影像的幾何精校正中，平原影像糾正精度達到2個像元以內，山區達到3個像元以內[5]。上述方法建立校正模型后，都是基于最小二乘的平差方法解算衛星圖像的幾何誤差。

對于Landsat-8衛星的幾何精校正，待估參數包括衛星三維位置和速度、衛星三維姿態及其變化率，所以建立最小二乘平差模型，至少需要6個以上的控制點，若獲得穩定的幾何精校正精度，則需要更多分布均勻、數量充足的地面控制點作為多余觀測參與平差過程。而對于Landsat系列遙感衛星，一方面，地面采集控制點獲取的成本較高，尤其對于偏遠地區，人力、物力成本也相應增加[6]；另一方面，中分辨率(30m)影像在進行目視解譯或自動匹配同名點時，達到較高匹配精度更加不易。為了解決上述問題，本文提出了采用卡爾曼濾波的方法來進行Landsat-8衛星的幾何精校正。實驗證明，基于卡爾曼濾波的Landsat-8幾何精校正減少了對控制點數量的依賴，取得了較好的精校正效果。

2基本原理及其實現

2.1卡爾曼濾波的基本原理

卡爾曼濾波是一種遞推線性最小方差估計算法。它的基本思想是以最小均方誤差為最佳估計準則，采用信號與噪聲的狀態空間模型，在不同狀態變化過程中，利用前一狀態的估計值和當前狀態的觀測值來更新狀態變量的估計，求出當前狀態的估計值。算法根據系統方程和測量方程，對需要處理的狀態矢量做出滿足最小均方誤差的估計[7]。

卡爾曼濾波不但計算簡便，運算速率快，而且在誤差模型具體形式不確定的情況下，仍然具有適用性。因此，該算法在平穩的動態運動系統中已經被廣泛應用，尤其在GPS動態數據處理、組合導航和通信與信號處理等領域得到了深入的研究和應用[8-10]。

2.2基于卡爾曼濾波的Landsat-8幾何精校正原理

在遙感衛星成像過程中，衛星狀態空間的變化可以近似與時間線性相關。因此卡爾曼濾波在修正遙感衛星狀態方面也具有適用性，梁澤環論證了將卡爾曼濾波應用于SPOT衛星幾何精校正的可行性[11]，Shin等采用擴展卡爾曼濾波進行ATSR圖像幾何精校正能夠達到一個像元內的幾何精度[12]。本文提出用卡爾曼濾波算法進行Landsat-8影像幾何精校正。

將卡爾曼濾波應用于Landsat-8衛星幾何精校正的總流程如圖1所示。在控制點數據和系統級校正的結果圖像匹配之后，利用卡爾曼濾波建立衛星星歷數據和姿態數據的運動方程和測量方程，通過引入控制點，將系統級幾何校正后的星歷數據和姿態數據中存在的系統誤差進行修正，并最終處理得到幾何精校正的結果。

圖1　基于卡爾曼濾波的Landsat-8衛星影像 幾何精校正總流程

圖2　基于卡爾曼濾波的Landsat-8幾何精校正具體流程

基于卡爾曼濾波的Landsat-8幾何精校正是一個迭代求解過程，其計算分為3個過程：預估過程、校正過程與迭代過程。

預估過程：建立運動方程，利用K狀態下Landsat-8衛星狀態矢量的估計值及對應的誤差協方差矩陣進行預估，計算K+1狀態下Landsat-8衛星狀態矢量的預測值及對應的誤差協方差矩陣。

校正過程：該過程分為兩部分：①建立測量方程。在K+1狀態下，利用Landsat-8衛星狀態矢量預測值及對應的誤差協方差矩陣和引入的控制點坐標，求出該狀態下的卡爾曼增益；②建立更新方程。利用Landsat-8衛星狀態的預測值、引入的控制點坐標及卡爾曼增益，求解K+1狀態下Landsat-8衛星狀態矢量的估計值及對應的誤差協方差矩陣。

迭代過程：若K+1狀態下的Landsat-8衛星狀態矢量的變化量低于預設閾值，則迭代終止，并利用該狀態下的衛星狀態矢量獲得精確幾何校正模型文件，處理得到精校正結果；否則，加入新的控制點，令K=K+1，將其作為下一狀態的輸入值，重復進行預估過程和校正過程。在實際應用中，迭代過程往往采用控制點個數控制迭代次數，作為迭代結束條件，在本文中采用了此做法。

2.3基于卡爾曼濾波的Landsat-8幾何精校正關鍵技術

基于卡爾曼濾波的Landsat-8幾何精校正算法包括3個基本方程：運動方程、測量方程和更新方程。其中，運動方程對應預估過程，測量方程和更新方程對應校正過程。

2.3.1運動方程

運動方程用于Landsat-8幾何精校正中不同狀態下的衛星狀態矢量的轉換，是預估過程的主要方程。運動方程利用K狀態下Landsat-8衛星狀態矢量的估計值及其誤差協方差，推算K+1狀態Landsat-8衛星狀態矢量的預測值及其誤差協方差，其過程如圖3所示。

衛星狀態矢量的預測值計算公式為：

X(K+1/K)=φ(K+1/K)·X(K/K)+

G(K+1)·W(K+1)

(1)

式中，X(K+1/K)是由K狀態下衛星狀態矢量的估計值預測得到的K+1狀態下衛星狀態矢量的預測值，φ(K+1/K)是由K狀態到K+1狀態的轉移矩陣，X(K/K)是K狀態下衛星狀態矢量的估計值，W(K+1)是白噪聲向量，G(K+1)是噪聲矢量系數。

預測值的誤差協方差矩陣計算公式為：

P(K+1/K)=φ(K+1/K)·P(K/K)·

φ(K+1/K)T+Q(K+1)

(2)

式中，P(K+1/K)是X(K+1/K)對應的誤差協方差矩陣，P(K/K)是K狀態下X(K/K)對應的誤差協方差矩陣，X(K/K)是噪聲的協方差矩陣。

圖3　運動方程

運動方程形式的推導不涉及到幾何處理過程，僅與衛星動力學方程和姿態動力學方程相關[11]。轉移矩陣包括兩部分：與衛星星歷有關的部分和與衛星姿態有關的部分，其形式推導的流程如圖4所示。

衛星動力學方程是基于萬有引力定律和開普勒衛星運動三定律，在地心慣性坐標系下建立的。由于建立的方程是矢量方程，需要將矢量方程分解為坐標系3個方向上的分量組成方程組，并求解方程組的通解，最后經過簡化后，得到由K狀態預測K+1狀態的轉移矩陣中與星歷有關的矩陣部分φ(K+1/K)ephemeris，其形式如下：

(3)

式中，tK是引入的第K個控制點對應的Landsat-8影像同名點的成像時刻與景中心成像時刻的時間差，控制點同名像點的成像時刻從控制點匹配文件中讀取，景中心成像時刻從系統級校正結果文件中獲取；ω0是Landsat-8衛星在慣性坐標系下的角速度，其數值可在Landsat-8衛星參數文件中獲取。

姿態動力學方程基于剛體的動量矩定理建立，即剛體對慣性空間某固定點的角動量的變化率，等于作用于剛體的所有外力對此點力矩的總和。同樣，將矢量方程分解為3個方向上的分量方程，建立方程組，求解方程，可得由K狀態預測K+1狀態的轉移矩陣中與星歷有關的矩陣部分φ(K+1/K)attitude，其形式如下：

(4)

式中，tK和ω0的含義和來源與式(3)相同。

由K狀態預測K+1狀態的轉移矩陣φ(K+1/K)形式如式。

結合式(3)和式(4)，即：

(5)

2.3.2測量方程

測量方程的作用是聯系衛星狀態矢量和控制點坐標，利用測量矩陣，將運動方程中獲得的當前狀態下Landsat-8衛星狀態矢量的預測值，轉化為大地坐標系下控制點坐標的計算值。其過程如圖5所示。

圖5　測量方程

(6)

測量矩陣H(K)中的元素，可根據其定義直接導出[13]。

H(K)=

(7)

式中，(x,y,z)是引入的控制點在地心坐標系的坐標值，(φ,λ)是引入的控制點大地經緯度，其數值都可以從控制點匹配結果文件中獲取；B是Landsat-8衛星軌道傾角，h是Landsat-8衛星的地面高度，兩者的數值都可在Landsat-8衛星參數文件中獲取。

2.3.3更新方程

更新方程的作用，是根據運動方程中獲得的當前狀態下Landsat-8衛星狀態矢量的誤差協方差矩陣和測量矩陣計算當前狀態的卡爾曼增益，并利用當前狀態下，運動方程中得到的Landsat-8衛星狀態矢量的預測值和控制點大地坐標測量值，計算該狀態下Landsat-8衛星狀態矢量的估計值及誤差協方差矩陣，處理過程如圖6所示。

Kg(K+1)=P(K+1/K)·H(K+1)T/

{H(K+1)·P(K+1/K)·H(K+1)T+R(K+1)}-1

(8)

X(K+1/K+1)=X(K+1/K)+Kg(K+1)·

(9)

P(K+1/K+1)=P(K+1/K)-Kg(K+1)·

H(K+1)·P(K+1/K)

(10)

式(8)中Kg(K+1)為K+1狀態下的卡爾曼增益，R(K+1)是白噪聲矩陣；式(9)中X(K+1/K+1)是該狀態下衛星狀態矢量的估計值，Z(K+1)是該狀態下引入的控制點大地坐標測量值；式(10)中P(K+1/K+1)是X(K+1/K+1)對應狀態下衛星狀態矢量的誤差協方差矩陣。

圖6　更新方程

2.4實現步驟

選取Landsat-8衛星影像景中心時刻的衛星星歷和姿態數據，作為Landsat-8衛星狀態矢量預測值的初始值，此時，K初始值為0。

②利用協方差矩陣P(K+1/K)和測量矩陣H(K+1)計算出卡爾曼濾波的增益值Kg(K+1)。

③將卡爾曼濾波增益Kg(K+1)帶入到更新方程中，利用第K+1狀態下Landsat-8衛星狀態矢量預測值X(K+1/K)和控制點大地坐標的測量值Z(K+1)，求出第K+1狀態下Landsat-8衛星狀態矢量的估計值X(K+1/K+1)及其對應的誤差協方差矩陣P(K+1/K+1)。

④判斷算法終止條件，若迭代結束，則可將誤差校正后的星歷及姿態數據生成幾何精校正模型文件，再通過后續處理獲得幾何精校正圖像；若不滿足迭代結束條件，則令K=K+1。

⑤計算轉移矩陣，將其帶入運動方程中，利用第狀態下Landsat-8衛星狀態矢量的估計值預估第狀態下Landsat-8衛星狀態矢量的預測值，并求出相應的誤差協方差矩陣，并轉入①繼續迭代過程。

3實驗與討論

實驗處理了多景Landsat-8衛星圖像，限于篇幅本文選取了其中兩景Landsat-8衛星的圖像數據的實驗做介紹。第一景圖像對應地理范圍內的主要地形為平原，第二景圖像對應地理范圍內的主要地形為山地。兩景Landsat-8衛星的圖像數據分別在不同控制點數量的條件下，利用卡爾曼濾波算法進行幾何精校正，并將該算法與傳統最小二乘平差法進行幾何精校正的算法做比較。兩景Landsat-8衛星圖像信息及選取的控制點信息如表1所示。精度評價所使用的50個檢查點來自GLS2005控制點庫[14]。

兩組數據中，為說明卡爾曼濾波在迭代過程中誤差修正的效果，將每次迭代處理生成結果圖像，均利用檢查點計算RMSE(均方根誤差)，評價精校正精度。

表1　實驗圖像和控制點信息表

在兩個實驗數據的圖像中選取的控制點的分布情況如圖7和圖10所示。

圖7　數據一選取控制點分布圖

在不同的控制點引入順序下，基于卡爾曼濾波的Landsat-8圖像幾何精校正結果圖像的精度隨控制點引入個數變化情況如表2、表3和表4所示。基于卡爾曼濾波算法的和最小二乘法的Landsat-8圖像幾何精校正結果圖像的精度隨控制點引入個數變化情況的對比如圖8、圖9和圖11所示。其中，控制點為0時表示初始輸入系統級校正圖像的幾何精度。

表2　數據一卡爾曼濾波算法結果圖像精度隨控制點

圖8　數據一中兩種算法結果圖像精度與控制點 個數關系對比圖一

控制點數引入控制點的編號和順序總體RMSE0-14.05977100913.0743292009、00311.0854123009、003、0079.7813644009、003、007、0018.4539875009、003、007、001、0057.6759416009、003、007、001、005、0047.1713687009、003、007、001、005、004、0106.8204658009、003、007、001、005、004、010、0026.6081269009、003、007、001、005、004、010、002、0086.60761910009、003、007、001、005、004、010、002、008、0066.607619

圖9　數據一中兩種算法結果圖像精度與控制點 個數關系對比圖二

對比表2和表3可知，基于卡爾曼濾波的Landsat-8幾何精校正結果精度與控制點的引入順序無關。如圖9，圖10和圖11所示，最小二乘算法由于受到參數個數估計限制，至少需要7個控制點才能進行計算。相比而言，卡爾曼濾波算法則不受控制點個數的限制，隨著控制點個數增加，誤差減少，校正精度逐漸提高。當平原區控制點數達到8個(山地區控制點數達到10個)時，圖像精度已基本趨于穩定，其精度較系統級幾何校正精度提高40%以上，相較于最小二乘平差法，精度平均提高20%以上。當控制點數目繼續增加時，卡爾曼濾波算法結果的精度已經無明顯變化，而最小二乘算法結果誤差仍未達到最低。最終，當最小二乘方法引入22個控制點時，精校正精度趨于穩定，并與卡爾曼濾波算法結果精度相當，表明卡爾曼濾波算法收斂速度快于最小二乘法。

圖10　數據二選取控制點分布圖

控制點數引入控制點的編號和順序總體RMSE0-22.66217100819.410992008、00217.187773008、002、00115.442784008、002、001、00913.590875008、002、001、009、01012.031256008、002、001、009、010、00610.285557008、002、001、009、010、006、0038.7645378008、002、001、009、010、006、003、0047.7732069008、002、001、009、010、006、003、004、0077.10131010008、002、001、009、010、006、003、004、007、0056.620061

圖11　數據二中兩種算法結果圖像精度與 控制點個數關系對比圖

4結束語

Landsat-8衛星數據幾何精校正中，本文將卡爾曼濾波算法用于衛星姿態及星歷中的誤差校正。實驗證明，當控制點匹配正確且分布均勻時，當其數目低于7個時，基于卡爾曼濾波的幾何精校正方法能夠獲得較系統級校正結果精度更高的幾何校正結果。隨著控制點引入數目的增多，在迭代結束，定位精度穩定情況下，基于卡爾曼濾波的幾何精校正較傳統最小二乘算法達到誤差最小時的速度更快，且結果的校正精度相當。因此，基于卡爾曼濾波的Landsat-8衛星影像的幾何精校正較傳統最小二乘平差方法有效降低了對控制點數量的依賴，并且適用于不同地形條件。理論上，對于其他遙感衛星，如果能夠提供算法所需要的參數，此幾何精校正技術也同樣適用于該衛星的幾何精校正。

參考文獻：

[1]孫家抦.遙感原理與應用(第二版)[M].武漢：武漢大學出版社，2009.

[2]YOICHI S，HOMMA K，KOMURA F.Geometric correction algorithms for satellite imagery using a bi-directional scanning sensor[J].Geoscience and Remote Sensing，1991，29(2)：292-299.

[3]KRATKY V.Rigorous photogrammetric processing of SPOT images at CCM Canada[J].ISPRS Journal of Photogrammetry and Remote Sensing，1989，(44)：53-71.

[4]GRUEN A.Potentian and limitation of high resolution satellite imagery[C].The 21stAsian Conference on Remote Sensing，Taipei，2000.

[5]張過，李揚，祝小勇，等.有理函數模型在光學衛星影像幾何糾正中的應用[J].航天返回與遙感，2010，31(4)：51-57.

[6]肖倩，馮鐘葵，李山山.GLS2005控制點庫介紹及應用[J].遙感信息，2013，28(5)：59-63.

[7]彭丁聰.卡爾曼濾波的基本原理及應用[J].軟件導刊，2009，8(11)：32-34.

[8]秦永元，張洪鉞，汪淑華.卡爾曼濾波與組合導航原理[M].西安：西北工業大學出版社，2012.

[9]王力循.無跡Kalman濾波在IMU和GPS組合導航系統中的應用研究[D].南昌：南昌大學，2010.

[10]邱鳳云.Kalman濾波理論及其在通信與信號處理中的應用[D].濟南：山東大學，2008.

[11]梁澤環.卡爾曼濾波在衛星遙感影像大地校正中的應用[J].環境遙感，1990，5(4)：297-307.

[12]SHIN D，POLLARD J K，MULLER J P.Accurate geometric correction of ATSR images[J].IEEE Transactions on Geoscience and Remote Sensing，1997，35(4)：997-1006.

[13]ARNOLD P.Application of a recursive distortion estimator to the geodetic correction of THEMATIC MAPPER imagery[C].Machine Processing of Remotely Sensed Data Symposium，332-385.

[14]肖倩，馮鐘葵，李山山.GLS2005控制點庫介紹及應用[J].遙感信息，2013，28(5)：59-63.

E-mail：yuying5559104@163.com