针对陆地目标的激光测高仪大气延迟测距误差修正方法

摘要

本发明属于激光遥感技术领域，特别涉及一种针对陆地目标的星载激光测高系统大气延迟误差修正方法。本发明首先推导了大气延迟误差与地表大气压强之间的理论关系式，进而简化了压高方程，根据地面气象站测得的气象站所在高度气压、海平面高度气压，计算每个气象站气压衰减因子，建立气压随高度变化的衰减模型。依据激光脚点的经纬度选取距离最近的k个气象站，根据气压衰减模型，计算每个气象站所在位置在待插值点海拔高度处的气压；根据反距离加权算法，计算待插值点地理位置处的气压。交叉验证结果表明，该方法的气压时空内插的精度高于目前基于NECP再分析资料的计算方法，在我国激光测高卫星的测量精度提高方面具有重要应用价值。

说明书

技术领域

本发明涉及激光遥感技术领域，尤其是涉及卫星激光测高仪针对地球表面陆地目标观测时大气延迟测距误差的修正方法。

背景技术

星载激光测高仪是一种安置于低轨卫星平台，用于获取地表高程信息的新型卫星载荷。星载激光测高仪通过记录激光脉冲从发射到被反射接收的渡越时间，计算卫星平台与地面激光脚点之间的距离；通过卫星平台的姿态、位置进而计算被测地物目标的平面和高程三维坐标。大气延迟测距误差是指当激光脉冲穿过地球大气层时，由于大气折射率不均匀产生的激光脉冲传输延迟，这对激光渡越时间的测量精度以及地球表面三维轮廓的测量造成显著的影响。

根据静力学方程推导可得，大气延迟测距误差由地表气压主导，地面每1hPa大气压将导致2.30mm的误差，而地表气压平均值约为1000hPa，由此带来的大气延迟测距误差约为2.3m。地表气压数据通常根据气象站测量数据通过相应的算法计算获得；激光测高仪随卫星飞行过程中，其星下激光脚点分布于地球表面的不同位置，且测量时间间隔通常小于1秒；任何国家和地区气象站实测数据都不足以满足卫星激光测高仪计算大气延迟的时间和空间网格密度要求，因此需要专门研究地表气压时空内插算法计算每个激光脉冲对应时刻和位置的地表气压数值，进而计算大气延迟测距误差的改正值。

世界上首台具有全波形记录能力的对地观测星载激光测高仪是美国2003年发射的ICESat卫星上搭载的地球科学激光测高系统GLAS(Geoscience Laser AltimeterSystem)。目前内插计算地表气压方法仅有美国GLAS系统使用的基于NCEP气象数据集的时空内插方法。该方法对于GLAS系统观测核心要素相对平坦的南北极冰盖区域精度足够(南极冰盖90％地区坡度小于3:200，使用NCEP基于标准气压层的全球气象数据时空线性/双线性内插误差的标准差为5hPa(对应干项延迟12mm))，但对于地表起伏复杂的陆地目标精度较低(约5cm，相对于GLAS系统整体设计精度15cm而言已经严重超标)。

而我国已发射的ZY3-02号卫星与未来发射的GF-7号卫星为激光/光学立体测绘卫星，以生产大比例尺地形图为目的，其主要以国内陆地地表为首要观测目标，获取高精度的地表高程数据。因此，针对陆地目标的高精度星载激光大气测距误差修正方法具有现实的应用意义。

发明内容

本发明主要是针对陆地目标的星载激光测高系统大气延迟误差，提供一种新的修正方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种针对陆地目标的星载激光测高系统大气延迟误差修正方法，包括如下步骤：

步骤1，以激光脚点经纬度坐标作为气压待插值点；以激光脚点为中心，确定距离激光脚点较近的k个气象站，作为采样点。

步骤2，根据步骤1中确定的气象站，利用气象站位置的海平面高度气压P_i0、气象站高度气压P_i以及气象站高度h_i，计算气象站所在位置气压衰减因子α_i。

步骤3，通过步骤2中采样点即每个气象站的气压衰减因子α_i，依据待插值点的海拔高度h，计算第i个气象站在海拔高度h处的气压值P_ih。

步骤4，重复步骤2-3，计算所有k个采样气象站位置在海拔高度h处的气压值P_ih(i＝1,2,3…k)。

步骤5，利用步骤4计算的所有采样气象站的气压值P_ih，根据空间反距离加权算法，计算t1时刻待插值位置处的气压值P_t1。

步骤6，重复步骤2-5，以1小时为间隔，依据t2时刻气象站所采集的气象数据计算t2时刻待插值位置处气压P_t2。

步骤7，根据星载激光测高系统测量UTC时间，在P_t1和P_t1之间进行时间线性内插，计算星载激光测高系统激光脚点位置的地表大气压P_SURF。

步骤8，计算天顶方向的大气延迟测距误差数值，并将该数值与光束指向角确定的映射函数相乘，计算任意指向角时激光测高大气延迟改正值ΔL。

所述的步骤1中，作为采样点的气象站按照如下方法选取：

a.以待插值点作为中心，计算距离待插值点最近的l个气象站。

b.根据待插值点的海拔高度h，确定l个气象站中海拔高度位于0.5h至1.5h之间的m个气象站，并选取m个气象站中距离待插值点中心最近的k个气象站。

所述步骤2，气象站气压衰减因子α_i求解方法如下：

根据步骤1中k个气象站的海拔高度，气象站平均海平面处气压P_i0，气象站高度处的气压P_i，可以利用公式(17)计算每个气象站所在位置处气压随高度的衰减因子α_i，其中i＝1,2,…k；

α＝1/z·ln(P_z/P₀)(17)

其中，α为衰减因子，P₀是气象站在海平面高度的气压，P_z为气象站在高度为z处的气压。

所述步骤3中，根据公式(18)，计算每个气象站在海拔高度为h处的气压值P_ih，其中i＝1,2,…k；

其中，P_i0为第i个气象站在海平面高度上的气压。

本发明主要有以下优点：

1)在气象数据来源方面，NCEP气象数据空间网格间隔为1°×1°，时间采样间隔为6小时，以模式函数计算的基于标准气压层格式，并不是当地对应的地面实测数据。国内气象站无论是在时空覆盖密度上，还是在气象数据精度方面都高于NCEP数据。

2)在数据处理算法方法方面，传统基于NCEP标准气压层的方法，在计算给定高度处的气压时需要使用4阶数值积分方式，算法复杂度大于本发明压高方程的方法。

3)在空间数据内插方法方面，传统方法仅考虑距离激光脚点地理位置最近的4个网格顶点的NCEP气象数据，使用空间双线性内插方法，而实际上，地表气压主要与海拔高度相关，本发明综合考虑海拔高度和气象站距离两个因素，采用反距离加权空间内插方法。以上这些优点，都使得采用本发明的方法修正精度比传统基于NCEP的修正方法有明显提高。

附图说明

图1是本发明中涉及到的属于基本站的气象站分布，气象站的位置以及海拔高度数据来源于中国气象数据网(http://data.cma.cn)，从示意图中可以看出，采用地面气象站实测数据，具有更高的空间密度。

图2是采用步骤1中方法选取采样点的结果。

图3与图2中待插值点以及选取的6个气象站相对应。其中，P点表示的是待插值点位置以及实测气压值；P1至P6表示的是图2中的6个被选中的气象站信息。数据为中国气象数据网公布的河南省区域内气象站2016年3月16日2时测量的气压。

图4是通过交叉验证的方法，根据反距离加权内插计算出来的位于河南省的114个气象站所在位置的气压值与实际值比较。其中‘o’点表示气象站实测气压值，方框表示采用交叉验证方法，进行内插计算的结果。

图5根据河南省2016年9月24日2时、8时、14时、20时的气压进行空间内插交叉验证计算得到的气压内插误差。

图6是根据公式计算得到的采用本发明方法进行星载激光测高系统大气延迟误差修正，其中的平均误差是每个气象站2016年9月24日到2016年9月28日之间的20次误差平均值。

图7是将河南省范围内气象站实测的2016年9月24日2时和8时的气压作为已知值，分别通过时空内插方法计算这期间的3时、4时、5时的气压值并与实际气压比较得误差如图所示。

图8是根据不同的时间密度进行插值的误差。其中dt＝2表示根据2016年9月24日2时和4时的实测数据进行时空内插计算3时的气压误差值；其中dt＝4表示根据2016年9月24日2时和6时的实测数据进行时空内插计算4时的气压误差值；其中dt＝6表示根据2016年9月24日2时和8时的实测数据进行时空内插计算5时的气压误差值；

图9是分别计算每个区域内气象站空间内插的误差平均值以及对应的大气延迟修正误差。

图10是在河北地区分别使用本方法以及使用NCEP地表气压内插方法的误差比较。图示结果是每个气象站2016年9月27日14时大气延迟修正误差。绝对误差。采用NCEP数据进行大气延迟修正的误差均值为：46mm；采用地面数据进行大气延迟修正的误差均值为：1.7mm。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案做进一步的具体的说明。

实施例：

首先介绍一下本发明所需要的理论基础：

1.星载激光测距大气延迟误差与被测地表气压的理论关系式

当激光脉冲穿过地球大气层时，由于大气折射率不均匀产生的激光脉冲传输延迟，对激光测距精度造成显著影响。激光测高系统单程大气折射延迟表示为：

其中，ΔL是单程大气延迟测距误差值，n(s)是沿激光脉冲传播路径s对应的大气折射率，S_ATM是实际发射激光脉冲从卫星传播到地球表面所穿越的曲线路径，S_VAC是从卫星沿着激光脉冲方位角方向到地球反射表面的直线路径。第二项积分计算只需已知卫星测量参考点和激光脚点的坐标位置；而第一项积分计算需要已知沿发射激光脉冲传播曲线路径沿途的大气折射率，用光线追踪方法都可以得到非常精确的结果。但对于大量的数据而言，这种方法需要的时间太长。因此，通常使用沿天顶方向的大气延迟测距误差值与映射函数乘积的方法来计算不同高度角时的大气延迟：

ΔL＝m(ε)ΔL_Z(2)

其中，ΔL_z为天顶方向大气延迟测距误差，积分项是由沿天顶方向从地面激光脚点loc到卫星测量参考点z的直线积分得到；m(ε)是与高度角ε相关的映射函数。

大气延迟修正解算模型中的积分变量(n-1)是光线在当前大气中的折射率与在真空中的折射率两者之差，它通常以百万分之一为单位来给定，即(n-1)＝10^-6N。对于1.064μm的近红外波段，N可以表示为：

其中k₁(λ)和k₂(λ)是波长λ的通用经验公式，对于1064nm波长的激光测高仪系统，k₁＝0.7866070K/Pa，k₂＝0.6644364K/Pa。P_d和P_w是干空气和水汽的分压强。Z_d和Z_w表征了干空气和水汽的压缩率，并遵从非理想气体法则。ρ_d是干空气气体密度，ρ_w是空气中水汽气体密度，总气体密度ρ＝ρ_d+ρ_w，得出：

式中，M_d和M_w分别为干空气和水汽的分子量，T为温度，R为气体常数，第一项为大气的静力学部分(与大气中的干空气有关，又称干项延迟)，该项值对应的测距误差值可以精确计算出。运用静力学方程dP/dz＝-ρ(z)g(z)，其中dP/dz是以高度z为函数的大气压强P的微分，ρ(z)为大气密度随高度变化的函数，g(z)为重力加速度随高度变化的函数。因此，大气静力延迟项ΔL_H可以表示为：

其中，g_m是大气中的平均重力加速度。重力加速度随着高度增加而缓慢减小，可以简单估计为纬度和高度z的函数。因为积分项仅仅是高度z处的表面压强，因此，天顶方向的大气延迟可以表示为：

其中P_SURF是地表大气压强。天顶延迟的另一部分为静力项延迟没有包括的水汽部分，也称“湿项”ΔL_W，可以表示为：

其中k₂ˊ＝k₂-k₁M_w/M_d。积分项为总可降水气量PW，是大气模型中经常用到的大气变量。

参数k₁和k₂可以根据经验方程来解算，对于激光波长为1.064μm的激光测高仪系统而言，参数k₁＝0.80277K/Pa，k₂＝0.66388K/Pa。干空气分子量M_d＝28.9644g/mol，水汽分子量M_w＝18.0152g/mol，则k₂ˊ＝0.16458K/Pa。因此，最后干项、湿项分别以及天顶方向的总延迟为：

ΔL_W＝(7.620×10^-5m/(kg/m²))PW(10)

ΔL_Z＝ΔL_H+ΔL_W(11)

假定地表平均气压为1000hPa，则总延迟的主要部分：天顶方向大气干项延迟测距误差约为2.35m。可降水量的值从极地的不到10mm到热带地区的50mm不等，天顶方向大气湿项延迟测距误差在1到4mm之间变化。相对于由干空气引起的干项延迟，湿项延迟几乎可以忽略；因此，天顶方向大气延迟测距误差可以直接使用(9)式替代(11)式计算得出，且只与地表气压P_SURF有关。

目前常见的复杂映射函数有Marini连分式模型、NMF模型、CfA2.2模型，假设对流层大气折射率为球对称，Marini给出一种常参数连分数形式的映射函数：

其中a，b，c…为待定常数，通常用气象数据来估计获得，目前根据高度角ε取值范围已有详细的查表法来确定这些参数。通常情况下，激光测高仪系统工作在近似天底入射的状态，高度角ε近似等于90°，可以用此方程最简单的形式近似表示，即a＝b＝c＝0。

在此种工作条件下，(12)式映射函数m(ε)约等于1，根据(2)式结论，激光测高仪的大气延迟测距误差可以使用(9)式近似计算得出。2.压高方程理论及其简化模型

为了精确地获得气压与高度的对应关系，通常将静力学方程从气层底部到顶部进行积分，即得出压高方程：

式中：Z₁为第1点高度值；Z₂为第2点高度值；P₁为Z₁点高度的气压值；P₂为Z₂点高度的气压值。该式表示任意两个高度上的气压差等于这两个高度间单位截面积空气柱的重量。

将理想气体状态方程P＝ρRT代入(14)式可得：

当海拔高度z₀＝0，则P₀为海平面气压，(15)式给出了气象站所在海拔高度z处的气压P_z与平均海平面高度的气压P₀的理论关系式。根据(15)式可知，衰减系数是由重力加速度g、温度T以及气体常数R沿海拔高度方向积分得到。假设重力加速度g、气温T以及气体常数R在以海拔高度z附近的一段高度范围内保持不变，则可认为衰减系数在该段范围内保持不变，(15)式可以简化为：

P_z＝P₀e^αz(16)

其中，α为衰减因子，如果已知气象站在海平面高度气压P₀和高度为z处的气压P_z，则可以通过(17)式计算海拔高度z附近的衰减因子。

α＝1/z·ln(P_z/P₀)(17)

根据待插值点的海拔高度h计算第i气象站在海拔高度h处的气压为：

其中，α_i第i个气象站的衰减因子；P_i0为第i个气象站在海平面高度上的气压。根据以上方法计算待插值位置的气压，比GLAS使用高阶数值积分方式更加简洁、高效。

3.反距离加权内插算法与时间线性内插

地球是一个近乎标准椭球体，它的赤道半径为6378.140km，极半径为6356.755km，平均半径为6371.004km。假设地球是一个完美的球体，则它的半径就是地球的平均半径，即为R_d。设第一点A的经纬度为第二点B的经纬度为按照0度经线的基准，东经取经度的正值，西经取经度负值。北纬取(90-Latitude)，南纬取(90+Latitude)。经过上述处理后的两点记为和则两点距离d表示为：

d＝R_d·arccos(C)·π/180(20)

根据以上公式可以计算待插值点与气象站之间的距离d，从而根据步骤1中说明，选择最近的k个气象站。

反距离加权内插是根据采样气象站与插值点之间的距离计算采样点的权重值，并且距离插值点越近的采样点，其权重越大。采样点权重计算公式为：

式中，w_i表示第i个采样点的权重，d_i表示第i个采样点与待插值点的距离。权指数u通常取2，这里每个采样点即所选取的k个气象站。若待插值位置海拔高度为h，距离所选取k个气象站中第i个气象站的距离为d_i，那么第i个气象站对应海拔高度为h的气压P_ih可以通过(18)式计算。待插值位置处的气压P通过与反距离加权内插的方式计算得出。

以上过程可以计算待插值点在气象数据测量时刻的气压值。采用反距离加权算法，与GLAS采用的空间线性内插相比，空间内插误差更小。

由于气象数据测量频率远大于星载激光测高系统的脉冲发射频率，因此通常需要进一步进行气压的时间内插，才可以获得激光脚点在t时刻的气压。若气象数据相邻两次的测量时刻分别为t₁和t₂，并且t处于这两个时刻之间，则根据以上公式可以计算得到待插值点在t₁和t₂时刻的气压值P_t1和P_t2，则通过线性内插可知，待插值位置处在t时刻的气压值P_SURF为：

4.实际大气延迟修正流程

星载激光大气延迟是由于折射率不均匀，通过静力学方程推导建立大气延迟测距误差与地表气压之间的定量关系式，对地面气象站测量的气压值进行时间和空间内插，可以获得激光脚点处的大气表面气压值，进而对星载激光大气延迟测距误差进行修正。具体流程如下：

a.根据激光脚点的经纬度根据公式(19)和(20)计算每个气象站与待插值点之间的距离d_i，并根据步骤1说明，选择最近6个气象站。

b.每个气象站有所在地理位置和海拔高度的测量气压P_z，以及所在地理位置海平面气压P₀，根据公式(17)分别计算第i个气象站所在位置的气压衰减因子α_i，并根据公式(18)计算气象站在待插值点海拔高度h上的气压P_ih。

c.根据公式(21)分别计算所选取各个气象站的权重，并根据公式(22)计算分别计算待插值点在t₁和t₂时刻的气压值P_t1和P_t2，并根据公式(23)计算t时刻待插值点地表气压P_SURF。

d.在通过加权内插计算得到地表气压P_SURF后，根据公式(9)计算星载激光测高系统在天顶方向大气延迟。

e.忽略大气湿项延迟的影响，根据公式(2)和(13)计算得到大气延迟测距误差的修正量。

本文中所述的具体实施例仅是对本发明精神做举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方法代替，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。