本發(fā)明涉及激光遙感技術領域，尤其是涉及衛(wèi)星激光測高儀針對地球表面陸地目標觀測時大氣延遲測距誤差的修正方法。

背景技術：

星載激光測高儀是一種安置于低軌衛(wèi)星平臺，用于獲取地表高程信息的新型衛(wèi)星載荷。星載激光測高儀通過記錄激光脈沖從發(fā)射到被反射接收的渡越時間，計算衛(wèi)星平臺與地面激光腳點之間的距離；通過衛(wèi)星平臺的姿態(tài)、位置進而計算被測地物目標的平面和高程三維坐標。大氣延遲測距誤差是指當激光脈沖穿過地球大氣層時，由于大氣折射率不均勻產生的激光脈沖傳輸延遲，這對激光渡越時間的測量精度以及地球表面三維輪廓的測量造成顯著的影響。

根據(jù)靜力學方程推導可得，大氣延遲測距誤差由地表氣壓主導，地面每1hpa大氣壓將導致2.30mm的誤差，而地表氣壓平均值約為1000hpa，由此帶來的大氣延遲測距誤差約為2.3m。地表氣壓數(shù)據(jù)通常根據(jù)氣象站測量數(shù)據(jù)通過相應的算法計算獲得；激光測高儀隨衛(wèi)星飛行過程中，其星下激光腳點分布于地球表面的不同位置，且測量時間間隔通常小于1秒；任何國家和地區(qū)氣象站實測數(shù)據(jù)都不足以滿足衛(wèi)星激光測高儀計算大氣延遲的時間和空間網(wǎng)格密度要求，因此需要專門研究地表氣壓時空內插算法計算每個激光脈沖對應時刻和位置的地表氣壓數(shù)值，進而計算大氣延遲測距誤差的改正值。

世界上首臺具有全波形記錄能力的對地觀測星載激光測高儀是美國2003年發(fā)射的icesat衛(wèi)星上搭載的地球科學激光測高系統(tǒng)glas(geosciencelaseraltimetersystem)。目前內插計算地表氣壓方法僅有美國glas系統(tǒng)使用的基于ncep氣象數(shù)據(jù)集的時空內插方法。該方法對于glas系統(tǒng)觀測核心要素相對平坦的南北極冰蓋區(qū)域精度足夠(南極冰蓋90％地區(qū)坡度小于3:200，使用ncep基于標準氣壓層的全球氣象數(shù)據(jù)時空線性/雙線性內插誤差的標準差為5hpa(對應干項延遲12mm))，但對于地表起伏復雜的陸地目標精度較低(約5cm，相對于glas系統(tǒng)整體設計精度15cm而言已經(jīng)嚴重超標)。

而我國已發(fā)射的zy3-02號衛(wèi)星與未來發(fā)射的gf-7號衛(wèi)星為激光/光學立體測繪衛(wèi)星，以生產大比例尺地形圖為目的，其主要以國內陸地地表為首要觀測目標，獲取高精度的地表高程數(shù)據(jù)。因此，針對陸地目標的高精度星載激光大氣測距誤差修正方法具有現(xiàn)實的應用意義。

技術實現(xiàn)要素：

本發(fā)明主要是針對陸地目標的星載激光測高系統(tǒng)大氣延遲誤差，提供一種新的修正方法。

本發(fā)明的上述技術問題主要是通過下述技術方案得以解決的：

一種針對陸地目標的星載激光測高系統(tǒng)大氣延遲誤差修正方法，包括如下步驟：

步驟1，以激光腳點經(jīng)緯度坐標作為氣壓待插值點；以激光腳點為中心，確定距離激光腳點較近的k個氣象站，作為采樣點。

步驟2，根據(jù)步驟1中確定的氣象站，利用氣象站位置的海平面高度氣壓pi0、氣象站高度氣壓pi以及氣象站高度hi，計算氣象站所在位置氣壓衰減因子αi。

步驟3，通過步驟2中采樣點即每個氣象站的氣壓衰減因子αi，依據(jù)待插值點的海拔高度h，計算第i個氣象站在海拔高度h處的氣壓值pih。

步驟4，重復步驟2-3，計算所有k個采樣氣象站位置在海拔高度h處的氣壓值pih(i＝1,2,3…k)。

步驟5，利用步驟4計算的所有采樣氣象站的氣壓值pih，根據(jù)空間反距離加權算法，計算t1時刻待插值位置處的氣壓值pt1。

步驟6，重復步驟2-5，以1小時為間隔，依據(jù)t2時刻氣象站所采集的氣象數(shù)據(jù)計算t2時刻待插值位置處氣壓pt2。

步驟7，根據(jù)星載激光測高系統(tǒng)測量utc時間，在pt1和pt1之間進行時間線性內插，計算星載激光測高系統(tǒng)激光腳點位置的地表大氣壓psurf。

步驟8，計算天頂方向的大氣延遲測距誤差數(shù)值，并將該數(shù)值與光束指向角確定的映射函數(shù)相乘，計算任意指向角時激光測高大氣延遲改正值δl。

所述的步驟1中，作為采樣點的氣象站按照如下方法選?。?/p>

a.以待插值點作為中心，計算距離待插值點最近的l個氣象站。

b.根據(jù)待插值點的海拔高度h，確定l個氣象站中海拔高度位于0.5h至1.5h之間的m個氣象站，并選取m個氣象站中距離待插值點中心最近的k個氣象站。

所述步驟2，氣象站氣壓衰減因子αi求解方法如下：

根據(jù)步驟1中k個氣象站的海拔高度，氣象站平均海平面處氣壓pi0，氣象站高度處的氣壓pi，可以利用公式(17)計算每個氣象站所在位置處氣壓隨高度的衰減因子αi，其中i＝1,2,…k；

α＝1/z·ln(pz/p0)(17)

其中，α為衰減因子，p0是氣象站在海平面高度的氣壓，pz為氣象站在高度為z處的氣壓。

所述步驟3中，根據(jù)公式(18)，計算每個氣象站在海拔高度為h處的氣壓值pih，其中i＝1,2,…k；

其中，pi0為第i個氣象站在海平面高度上的氣壓。

本發(fā)明主要有以下優(yōu)點：

1)在氣象數(shù)據(jù)來源方面，ncep氣象數(shù)據(jù)空間網(wǎng)格間隔為1°×1°，時間采樣間隔為6小時，以模式函數(shù)計算的基于標準氣壓層格式，并不是當?shù)貙牡孛鎸崪y數(shù)據(jù)。國內氣象站無論是在時空覆蓋密度上，還是在氣象數(shù)據(jù)精度方面都高于ncep數(shù)據(jù)。

2)在數(shù)據(jù)處理算法方法方面，傳統(tǒng)基于ncep標準氣壓層的方法，在計算給定高度處的氣壓時需要使用4階數(shù)值積分方式，算法復雜度大于本發(fā)明壓高方程的方法。

3)在空間數(shù)據(jù)內插方法方面，傳統(tǒng)方法僅考慮距離激光腳點地理位置最近的4個網(wǎng)格頂點的ncep氣象數(shù)據(jù)，使用空間雙線性內插方法，而實際上，地表氣壓主要與海拔高度相關，本發(fā)明綜合考慮海拔高度和氣象站距離兩個因素，采用反距離加權空間內插方法。以上這些優(yōu)點，都使得采用本發(fā)明的方法修正精度比傳統(tǒng)基于ncep的修正方法有明顯提高。

附圖說明

圖1是本發(fā)明中涉及到的屬于基本站的氣象站分布，氣象站的位置以及海拔高度數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn)，從示意圖中可以看出，采用地面氣象站實測數(shù)據(jù)，具有更高的空間密度。

圖2是采用步驟1中方法選取采樣點的結果。

圖3與圖2中待插值點以及選取的6個氣象站相對應。其中，p點表示的是待插值點位置以及實測氣壓值；p1至p6表示的是圖2中的6個被選中的氣象站信息。數(shù)據(jù)為中國氣象數(shù)據(jù)網(wǎng)公布的河南省區(qū)域內氣象站2016年3月16日2時測量的氣壓。

圖4是通過交叉驗證的方法，根據(jù)反距離加權內插計算出來的位于河南省的114個氣象站所在位置的氣壓值與實際值比較。其中‘o’點表示氣象站實測氣壓值，方框表示采用交叉驗證方法，進行內插計算的結果。

圖5根據(jù)河南省2016年9月24日2時、8時、14時、20時的氣壓進行空間內插交叉驗證計算得到的氣壓內插誤差。

圖6是根據(jù)公式計算得到的采用本發(fā)明方法進行星載激光測高系統(tǒng)大氣延遲誤差修正，其中的平均誤差是每個氣象站2016年9月24日到2016年9月28日之間的20次誤差平均值。

圖7是將河南省范圍內氣象站實測的2016年9月24日2時和8時的氣壓作為已知值，分別通過時空內插方法計算這期間的3時、4時、5時的氣壓值并與實際氣壓比較得誤差如圖所示。

圖8是根據(jù)不同的時間密度進行插值的誤差。其中dt＝2表示根據(jù)2016年9月24日2時和4時的實測數(shù)據(jù)進行時空內插計算3時的氣壓誤差值；其中dt＝4表示根據(jù)2016年9月24日2時和6時的實測數(shù)據(jù)進行時空內插計算4時的氣壓誤差值；其中dt＝6表示根據(jù)2016年9月24日2時和8時的實測數(shù)據(jù)進行時空內插計算5時的氣壓誤差值；

圖9是分別計算每個區(qū)域內氣象站空間內插的誤差平均值以及對應的大氣延遲修正誤差。

圖10是在河北地區(qū)分別使用本方法以及使用ncep地表氣壓內插方法的誤差比較。圖示結果是每個氣象站2016年9月27日14時大氣延遲修正誤差。絕對誤差。采用ncep數(shù)據(jù)進行大氣延遲修正的誤差均值為：46mm；采用地面數(shù)據(jù)進行大氣延遲修正的誤差均值為：1.7mm。

具體實施方式

下面通過實施例，并結合附圖，對本發(fā)明的技術方案做進一步的具體的說明。

實施例：

首先介紹一下本發(fā)明所需要的理論基礎：

1.星載激光測距大氣延遲誤差與被測地表氣壓的理論關系式

當激光脈沖穿過地球大氣層時，由于大氣折射率不均勻產生的激光脈沖傳輸延遲，對激光測距精度造成顯著影響。激光測高系統(tǒng)單程大氣折射延遲表示為：

其中，δl是單程大氣延遲測距誤差值，n(s)是沿激光脈沖傳播路徑s對應的大氣折射率，satm是實際發(fā)射激光脈沖從衛(wèi)星傳播到地球表面所穿越的曲線路徑，svac是從衛(wèi)星沿著激光脈沖方位角方向到地球反射表面的直線路徑。第二項積分計算只需已知衛(wèi)星測量參考點和激光腳點的坐標位置；而第一項積分計算需要已知沿發(fā)射激光脈沖傳播曲線路徑沿途的大氣折射率，用光線追蹤方法都可以得到非常精確的結果。但對于大量的數(shù)據(jù)而言，這種方法需要的時間太長。因此，通常使用沿天頂方向的大氣延遲測距誤差值與映射函數(shù)乘積的方法來計算不同高度角時的大氣延遲：

δl＝m(ε)δlz(2)

其中，δlz為天頂方向大氣延遲測距誤差，積分項是由沿天頂方向從地面激光腳點loc到衛(wèi)星測量參考點z的直線積分得到；m(ε)是與高度角ε相關的映射函數(shù)。

大氣延遲修正解算模型中的積分變量(n-1)是光線在當前大氣中的折射率與在真空中的折射率兩者之差，它通常以百萬分之一為單位來給定，即(n-1)＝10^-6n。對于1.064μm的近紅外波段，n可以表示為：

其中k1(λ)和k2(λ)是波長λ的通用經(jīng)驗公式，對于1064nm波長的激光測高儀系統(tǒng)，k1＝0.7866070k/pa，k2＝0.6644364k/pa。pd和pw是干空氣和水汽的分壓強。zd和zw表征了干空氣和水汽的壓縮率，并遵從非理想氣體法則。ρd是干空氣氣體密度，ρw是空氣中水汽氣體密度，總氣體密度ρ＝ρd+ρw，得出：

式中，md和mw分別為干空氣和水汽的分子量，t為溫度，r為氣體常數(shù)，第一項為大氣的靜力學部分(與大氣中的干空氣有關，又稱干項延遲)，該項值對應的測距誤差值可以精確計算出。運用靜力學方程dp/dz＝-ρ(z)g(z)，其中dp/dz是以高度z為函數(shù)的大氣壓強p的微分，ρ(z)為大氣密度隨高度變化的函數(shù)，g(z)為重力加速度隨高度變化的函數(shù)。因此，大氣靜力延遲項δlh可以表示為：

其中，gm是大氣中的平均重力加速度。重力加速度隨著高度增加而緩慢減小，可以簡單估計為緯度和高度z的函數(shù)。因為積分項僅僅是高度z處的表面壓強，因此，天頂方向的大氣延遲可以表示為：

其中psurf是地表大氣壓強。天頂延遲的另一部分為靜力項延遲沒有包括的水汽部分，也稱“濕項”δlw，可以表示為：

其中k2ˊ＝k2-k1mw/md。積分項為總可降水氣量pw，是大氣模型中經(jīng)常用到的大氣變量。

參數(shù)k1和k2可以根據(jù)經(jīng)驗方程來解算，對于激光波長為1.064μm的激光測高儀系統(tǒng)而言，參數(shù)k1＝0.80277k/pa，k2＝0.66388k/pa。干空氣分子量md＝28.9644g/mol，水汽分子量mw＝18.0152g/mol，則k2ˊ＝0.16458k/pa。因此，最后干項、濕項分別以及天頂方向的總延遲為：

δlw＝(7.620×10^-5m/(kg/m²))pw(10)

δlz＝δlh+δlw(11)

假定地表平均氣壓為1000hpa，則總延遲的主要部分：天頂方向大氣干項延遲測距誤差約為2.35m?？山邓康闹祻臉O地的不到10mm到熱帶地區(qū)的50mm不等，天頂方向大氣濕項延遲測距誤差在1到4mm之間變化。相對于由干空氣引起的干項延遲，濕項延遲幾乎可以忽略；因此，天頂方向大氣延遲測距誤差可以直接使用(9)式替代(11)式計算得出，且只與地表氣壓psurf有關。

目前常見的復雜映射函數(shù)有marini連分式模型、nmf模型、cfa2.2模型，假設對流層大氣折射率為球對稱，marini給出一種常參數(shù)連分數(shù)形式的映射函數(shù)：

其中a，b，c…為待定常數(shù)，通常用氣象數(shù)據(jù)來估計獲得，目前根據(jù)高度角ε取值范圍已有詳細的查表法來確定這些參數(shù)。通常情況下，激光測高儀系統(tǒng)工作在近似天底入射的狀態(tài)，高度角ε近似等于90°，可以用此方程最簡單的形式近似表示，即a＝b＝c＝0。

在此種工作條件下，(12)式映射函數(shù)m(ε)約等于1，根據(jù)(2)式結論，激光測高儀的大氣延遲測距誤差可以使用(9)式近似計算得出。2.壓高方程理論及其簡化模型

為了精確地獲得氣壓與高度的對應關系，通常將靜力學方程從氣層底部到頂部進行積分，即得出壓高方程：

式中：z1為第1點高度值；z2為第2點高度值；p1為z1點高度的氣壓值；p2為z2點高度的氣壓值。該式表示任意兩個高度上的氣壓差等于這兩個高度間單位截面積空氣柱的重量。

將理想氣體狀態(tài)方程p＝ρrt代入(14)式可得：

當海拔高度z0＝0，則p0為海平面氣壓，(15)式給出了氣象站所在海拔高度z處的氣壓pz與平均海平面高度的氣壓p0的理論關系式。根據(jù)(15)式可知，衰減系數(shù)是由重力加速度g、溫度t以及氣體常數(shù)r沿海拔高度方向積分得到。假設重力加速度g、氣溫t以及氣體常數(shù)r在以海拔高度z附近的一段高度范圍內保持不變，則可認為衰減系數(shù)在該段范圍內保持不變，(15)式可以簡化為：

pz＝p0e^αz(16)

其中，α為衰減因子，如果已知氣象站在海平面高度氣壓p0和高度為z處的氣壓pz，則可以通過(17)式計算海拔高度z附近的衰減因子。

α＝1/z·ln(pz/p0)(17)

根據(jù)待插值點的海拔高度h計算第i氣象站在海拔高度h處的氣壓為：

其中，αi第i個氣象站的衰減因子；pi0為第i個氣象站在海平面高度上的氣壓。根據(jù)以上方法計算待插值位置的氣壓，比glas使用高階數(shù)值積分方式更加簡潔、高效。

3.反距離加權內插算法與時間線性內插

地球是一個近乎標準橢球體，它的赤道半徑為6378.140km，極半徑為6356.755km，平均半徑為6371.004km。假設地球是一個完美的球體，則它的半徑就是地球的平均半徑，即為rd。設第一點a的經(jīng)緯度為第二點b的經(jīng)緯度為按照0度經(jīng)線的基準，東經(jīng)取經(jīng)度的正值，西經(jīng)取經(jīng)度負值。北緯取(90-latitude)，南緯取(90+latitude)。經(jīng)過上述處理后的兩點記為和則兩點距離d表示為：

d＝rd·arccos(c)·π/180(20)

根據(jù)以上公式可以計算待插值點與氣象站之間的距離d，從而根據(jù)步驟1中說明，選擇最近的k個氣象站。

反距離加權內插是根據(jù)采樣氣象站與插值點之間的距離計算采樣點的權重值，并且距離插值點越近的采樣點，其權重越大。采樣點權重計算公式為：

式中，wi表示第i個采樣點的權重，di表示第i個采樣點與待插值點的距離。權指數(shù)u通常取2，這里每個采樣點即所選取的k個氣象站。若待插值位置海拔高度為h，距離所選取k個氣象站中第i個氣象站的距離為di，那么第i個氣象站對應海拔高度為h的氣壓pih可以通過(18)式計算。待插值位置處的氣壓p通過與反距離加權內插的方式計算得出。

以上過程可以計算待插值點在氣象數(shù)據(jù)測量時刻的氣壓值。采用反距離加權算法，與glas采用的空間線性內插相比，空間內插誤差更小。

由于氣象數(shù)據(jù)測量頻率遠大于星載激光測高系統(tǒng)的脈沖發(fā)射頻率，因此通常需要進一步進行氣壓的時間內插，才可以獲得激光腳點在t時刻的氣壓。若氣象數(shù)據(jù)相鄰兩次的測量時刻分別為t1和t2，并且t處于這兩個時刻之間，則根據(jù)以上公式可以計算得到待插值點在t1和t2時刻的氣壓值pt1和pt2，則通過線性內插可知，待插值位置處在t時刻的氣壓值psurf為：

4.實際大氣延遲修正流程

星載激光大氣延遲是由于折射率不均勻，通過靜力學方程推導建立大氣延遲測距誤差與地表氣壓之間的定量關系式，對地面氣象站測量的氣壓值進行時間和空間內插，可以獲得激光腳點處的大氣表面氣壓值，進而對星載激光大氣延遲測距誤差進行修正。具體流程如下：

a.根據(jù)激光腳點的經(jīng)緯度根據(jù)公式(19)和(20)計算每個氣象站與待插值點之間的距離di，并根據(jù)步驟1說明，選擇最近6個氣象站。

b.每個氣象站有所在地理位置和海拔高度的測量氣壓pz，以及所在地理位置海平面氣壓p0，根據(jù)公式(17)分別計算第i個氣象站所在位置的氣壓衰減因子αi，并根據(jù)公式(18)計算氣象站在待插值點海拔高度h上的氣壓pih。

c.根據(jù)公式(21)分別計算所選取各個氣象站的權重，并根據(jù)公式(22)計算分別計算待插值點在t1和t2時刻的氣壓值pt1和pt2，并根據(jù)公式(23)計算t時刻待插值點地表氣壓psurf。

d.在通過加權內插計算得到地表氣壓psurf后，根據(jù)公式(9)計算星載激光測高系統(tǒng)在天頂方向大氣延遲。

e.忽略大氣濕項延遲的影響，根據(jù)公式(2)和(13)計算得到大氣延遲測距誤差的修正量。

本文中所述的具體實施例僅是對本發(fā)明精神做舉例說明。本發(fā)明所屬技術領域的技術人員可以對所描述的具體實施例做各種各樣的修改或補充或采用類似的方法代替，但并不會偏離本發(fā)明的精神或者超越所附權利要求書所定義的范圍。