CN103926589A - 星载激光测高系统固体地表目标平面和高程精度检测方法 - Google Patents

Abstract

本发明涉及激光遥感技术领域，尤其是涉及星载激光测高系统固体地表目标平面和高程精度检测方法。本发明通过建立星载激光测高系统测量固体地表的距离精度和高程精度数学模型，设计了完整的测高系统精度快速检测方法，并推导了具体实施过程中所涉及的回波能量计算方法和设计了地表目标参数计算方法。本发明不仅能够检测已知测高系统硬件参数和地表目标参数时的卫星激光参考点与地表目标之间的测距精度，而且结合卫星平台姿态和位置数据可以检测激光测高系统的平面和高程精度。这种平面和高程精度快速检测方法对激光测高系统参数设计和精度评估有较高的参考意义。

Description

星载激光测高系统固体地表目标平面和高程精度检测方法

技术领域

本发明涉及一种测高系统固体地表目标平面和高程精度检测方法，尤其是涉及一种星载激光测高系统固体地表目标平面和高程精度检测方法。

背景技术

星载激光测高系统是一种主动式测量系统，通过接收卫星平台激光器发出的激光脉冲经地表反射的微弱回波，精确计算激光脉冲在卫星和地面之间的渡越时间，得到卫星和地表的距离；结合卫星精密轨道和姿态数据，生成激光脚点精确地理位置和高程结果。通过卫星连续运行，最终得到覆盖地球表面DEM模型。

激光测高系统发射的激光脉冲能穿透植被获取三维地形，具有传统摄影测量方式无法取代的优点；其激光发散角在亚毫弧量级，水平定位精度和分辨率远远高于微波雷达方式；位于大气透过窗口的1064nm波长激光可以直接被冰盖和海洋表面反射，几乎没有穿透效应，高程精度可达15cm。这些优点使其广泛用于南北极冰盖变化监测、极区附近海冰变化监测、植被年际变化监测，少量用于海洋环境监测等领域；如果作为遥感影像的地面高程控制点，可以生成1:10000的大比例尺地形图。

ICESat卫星搭载的GLAS系统于2003年发射升空，在轨间断运行7年，是目前为止唯一一颗用于全球地表观测的星载激光测高系统。对地球观测的GLAS系统的观测成果不仅包含器件本身所带来的误差，也包含如光束穿越大气产生的散射和折射等环境误差，以及由地表斜率或粗糙度等引起的目标误差。因此，建立测高系统测距精度模型，以及完整的激光脚点平面和高程精度快速检测方法对激光测高系统参数设计、误差分析和数据产品的有效性和准确性评估至关重要。

然而，目前已有测距精度模型较为复杂，还没有全面分析器件、环境和目标对测距影响的简化数学模型，GLAS的技术文档中也没有给出完整描述星载激光测高系统平面和高程精度的数学模型，而是直接列出影响数值大小。

发明内容

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

一种星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，包括以下步骤：

步骤1，根据测高系统硬件参数计算受器件影响的距离误差ΔR_n；并设定大气透过率T_a、目标表面反射率β_r、以及其它测高系统参数如发射能量Q，接收望远镜面积A_R；

步骤2，输入地表目标高度模型数据，即沿X_D和Y_D均匀采样地表目标相对于X_DO_DY_D平面的高度；筛选得到光斑直径D范围内的所有高程点(x_n,y_n,z_n)；使用LM算法拟合筛选后所有高程点(x_n,y_n,z_n)的最接近平面，并计算平面法向量(a₁,a₂,-1)，进而计算目标表面斜率S；并通过直接求解所有的z_n标准差得到筛选后的光斑内目标表面粗糙度Std(Δξ)；最后通过步骤1设定的参数计算回波光子数N；这里，需要解释一下筛选的目的，由于激光测高系统发射的激光脉冲照射到地球表面时呈近似圆形分布光斑，其光斑直径为D，该原则是保留光斑内的高程点，因为光斑之外的点不对观测造成影响，如果不剔除相当于混入了噪声。

步骤3：使用步骤2得到的表面斜率S和粗糙度Std(Δξ)，回波光子数N，并设定测高系统探测器噪声系数F，以及光束指向角β，即天底方向与激光光轴夹角，计算目标影响导致的测距误差ΔR_pr、ΔR_ps和ΔR_pp；

步骤4：基于一个受环境影响的修正残差值ΔR_e，实际上环境误差影响有几米量级，而经过大气延迟修正和固体潮汐修正等通用算法修正后，剩余的环境残差即ΔR_e，通常为2-3cm；步骤1中得到的距离误差ΔR_n，以及步骤3中得到的测距误差ΔR_pr、ΔR_ps和ΔR_pp，根据误差合成理论，利用RSS方法(Root-Sum-Squares，平方和的根)计算由ΔR_e、ΔR_n、ΔR_e、ΔR_pr、ΔR_ps和ΔR_pp组成的整体测距误差ΔR；其中，ΔR_e通常取2-3cm；

步骤5：根据系统进行正常测量工作时三个姿态角ω、φ和κ；测高系统偏航Δω、横滚Δφ和俯仰角误差Δκ；指向角误差Δβ；三个方向GPS卫星定位误差ΔX_O、ΔY_O和ΔZ_O；以及步骤5得出整体测距误差ΔR；得到激光测高系统平面XY和高程Z方向精度。

本发明旨在简化并完善针对固体地表的激光测距精度模型，建立完整的激光脚点平面和高程精度模型，用于快速检测对地观测星载激光测高系统的精度指标。并且，利用丹麦科技大学提供的高程精度和空间分辨率更高的机载Lidar数据评估美国ICESat卫星GLAS星载激光测高系统实测数据的高程偏差，用于验证本发明方法的正确性和可行性。

在上述的星载激光测高系统固体地表目标平面和高程精度检测方法，所述步骤1中，距离误差ΔR_n是基于以下公式得到：

$\mathrm{\Δ}{R}_n=\frac{c}{2}\sqrt{(\frac{N_Q{T}_r^3}{12N^2\mathrm{\Δt}}+\frac{{\mathrm{\Δt}}^2}{12})+[\frac{{\mathrm{FN}}_B{T}_r^3}{12N^2\mathrm{\Δt}}+\frac{({\mathrm{FN}}_D+N_T+N_A)T_r^3}{12N^2\mathrm{\Δt}}]},$

其中，c为真空中光速，T_r为回波信号的全脉宽，通常取4倍回波信号脉宽σ_S，N_B、N_D、N_T、N_A和N_Q分别表示辐射背景噪声、探测器暗电流噪声、热噪声、前置放大噪声和量化噪声，F为探测器噪声系数，N为回波光子数。

在上述的星载激光测高系统固体地表目标平面和高程精度检测方法，所述步骤2中，筛选光斑直径D范围内所有高程点(x_n,y_n,z_n)是基于(x-X_D)²+(x-X_D)²≤D²/4原则来进行，所述目标表面斜率S是基于以下公式得到：

其中，a₁和a₂分别是步骤2中使用LM算法拟合所得平面法向量(a₁,a₂,-1)的前两项；

所述回波光子数N是基于以下公式得到：

其中，A_R为接收望远镜面积，T_a为大气透过率，Q为激光发射能量，β_r为目标表面反射率，η为接收系统能量利用率，β为光束指向角，h为普朗克常数，v为光波频率，S为目标表面斜率。

在上述的星载激光测高系统固体地表目标平面和高程精度检测方法，所述步骤3中，测距误差ΔR_pr、ΔR_ps和ΔR_pp，是基于以下公式得到：

其中，N为回波光子数；F为探测器噪声系数，通常在3～7之间；K_S为接收孔径面积与散斑相关面积的比值，K_S＝πA_R(2tanθ_T/λ)²，A_R为接收望远镜面积，θ_T为光束发散角，λ为激光波长；Std(Δξ)为目标表面粗糙度；S为目标表面斜率；β为光束指向角；

${\mathrm{\ΔR}}_{\mathrm{ps}}=\sqrt{(\frac{F}{N}+\frac{1}{2K_s})}\frac{z\tan {\mathrm{\θ}}_T}{\cos \mathrm{\β}}\tan (\mathrm{\β}+S),$ 其中，N为回波光子数；F为探测器噪声系数；K_S为接收孔径面积与散斑相关面积的比值；z为卫星轨道高度；θ_T为光束发散角；S为目标表面斜率；β为光束指向角；

其中，z为卫星轨道高度；S为目标表面斜率；β为光束指向角；Δβ为光束指向角误差；

ΔR_pr表示由地表粗糙度引起的测距误差；

ΔR_ps表示由地表斜率引起的测距误差；

ΔR_pp表示当地表斜率不为零时由指向角引起的额外距离误差。

在上述的星载激光测高系统固体地表目标平面和高程精度检测方法，所述步骤4中，使用RSS方法计算整体测距误差ΔR，即计算ΔR_e、ΔR_n、ΔR_pr、ΔR_ps和ΔR_pp平方和的根，具体公式为：

$\mathrm{\ΔR}=\sqrt{\mathrm{\Δ}{R}_e^2+\mathrm{\Δ}{R}_n^2+\mathrm{\Δ}{R}_{\mathrm{pr}}^2+\mathrm{\Δ}{R}_{\mathrm{ps}}^2+\mathrm{\Δ}{R}_{\mathrm{pp}}^2}$

在上述的星载激光测高系统固体地表目标平面和高程精度检测方法，所述步骤5中，测高系统的平面精度和高程精度基于以下公式获得：

平面精度的X方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δω、Δφ和Δκ分别表示卫星平台的偏航、俯仰和横滚角误差；ΔX_O表示卫星平台X方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值；

平面精度的Y方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δω、Δφ和Δκ分别表示卫星平台的偏航、俯仰和横滚角误差；ΔY_O表示卫星平台Y方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值；

高程精度的Z方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δφ和Δκ分别表示卫星平台俯仰和横滚角误差；ΔZ_O表示卫星平台Z方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值。

因此，本发明具有如下优点：1、综合考虑对地观测星载激光测高系统受环境、器件和地表目标影响，建立完整的平面和高程精度数学模型。2、简化了测距精度模型，使之能够快速计算和检测已知系统参数和地表目标时的距离精度、平面和高程精度。3、使用精度更高的机载Lidar数据对ICESat卫星GLAS系统进行了该方法可行性和正确性校验。

附图说明

图1是星载激光测高系统的工作原理示意图。

图2是本发明部分公式中变量所代表的参数名称。

图3是假设地表斜率一维分布时测距误差示意图。

图4是激光测高系统瞬时测量坐标系示意图。

图5是使用机载Lidar数据生成ICESat卫星GLAS测高系统脚点DEM示意图。

图6a是图5中其中一个目标表面较为平坦且粗糙度较小的激光脚点DEM放大显示，其具体目标参数和高程结果显示在图8中第一行。

图6b是图5中其中一个目标表面较为平坦且粗糙度较小的激光脚点DEM放大显示，其具体目标参数和高程结果显示在图8中第二行。

图7a是图5中其中一个目标表面有斜度且粗糙度较大的激光脚点DEM放大显示，其具体目标参数和高程结果显示在图8中第三行。

图7b是图5中其中一个目标表面有斜度且粗糙度较大的激光脚点DEM放大显示，其具体目标参数和高程结果显示在图8中第四行。

图8是图6和图7中脚点信息的详细数据分析。

图9是对所有脚点结果进行统计分析后的偏差曲线图。

具体实施方式

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

实施例：

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

1.回波能量理论推导

根据Gardner的理论，发射激光脉冲经过一次菲涅尔衍射，入射到地球表面经过目标反射，再经过一次菲涅尔衍射，到达接收望远镜视场前的光波复振幅表示为(1)式，测量地表时的几何位置关系如图1所示。

$u_r(r,z,t)=\frac{T_a}{\mathrm{\λz}}\exp \left[i(\mathrm{\ωt}-2\mathrm{kz}+\frac{{\mathrm{kr}}^2}{2z})\right]{\∫\∫\mathrm{\β}}^{1/2}\left(\mathrm{\ρ}\right)a(\mathrm{\ρ},z)f[t-\frac{2z}{c}-\frac{{\mathrm{\ρ}}^2}{\mathrm{cz}}+\frac{2\mathrm{\ξ}\left(\mathrm{\ρ}\right)}{c}]\exp \{-\mathrm{ik}[\frac{{\mathrm{\ρ}}^2}{2z}-2\mathrm{\ξ}\left(\mathrm{\ρ}\right)+\frac{\mathrm{\ρ}.r}{z}\left]\right\}{d}^2\mathrm{\ρ}$

(1)

其中参数含义如图2中所示，望远镜接收到的回波脉冲信号功率表示为(2)式，尖括号表示均值计算符号。

$P\left(t\right)=A_R<u_r(r,z,t){u^{*}}_r(r,z,t)>=\frac{A_R{T_a}^2}{z^2}\&Integral;{\&Integral;}_{\mathrm{\&Sigma;}}{\mathrm{\&beta;}}_r\left(\mathrm{\&rho;}\right)\&CenterDot;{\left|a(\mathrm{\&rho;},z)\right|}^2\&CenterDot;{\left|f\right[t-\frac{2z}{c}-\frac{{\mathrm{\&rho;}}^2}{\mathrm{cz}}+\frac{2\mathrm{\&xi;}\left(\mathrm{\&rho;}\right)}{c}\left]\right|}^2{d}^2\mathrm{\&rho;}$

(2)

通常情况下，激光发射脉冲在空域横截面和时域波形都近似为高斯函数，由(3)式和(4)式表示：

${\left|a(\mathrm{\&rho;},z)\right|}^2=\frac{Q}{2\mathrm{\&pi;}{\left(z\tan {\mathrm{\&theta;}}_T\right)}^2}\exp (-{\mathrm{\&rho;}}^2/2z^2{\tan }^2{\mathrm{\&theta;}}_T)---\left(3\right)$

${\left|f\left(t\right)\right|}^2=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_f}\exp (-\frac{t^2}{2{\mathrm{\&sigma;}}_f^2})---\left(4\right)$

(2)式中目标表面高程轮廓ξ(ρ)在平坦平面为0，在平坦斜面其均值为0；对于陆地和冰盖表面，可以近似认为是朗伯反射体，其反射率与入射角度有关。假设激光测高系统由天底方向入射地表平坦平面，且地表反射率β_r分布相同，测高系统探测器输出光子数N满足：

$N=\&Integral;\frac{\mathrm{\&eta;}}{\mathrm{h\&upsi;}}P\left(t\right)*h\left(t\right)\mathrm{dt}=\frac{A_R{T_a}^{2}Q{\mathrm{\&beta;}}_r\mathrm{\&eta;}}{\mathrm{\&pi;h\&upsi;}{z}^2}---\left(5\right)$

式中，h(t)表示接收系统脉冲响应函数，η表示接收系统能量利用率，η＝η_dη_r。对于非天顶方向入射的激光测高系统，重新建立坐标系以入射光轴方向为z’轴，此时与天底方向z轴夹角为β，光轴与地表平面交点为x’y’轴原点o’，y’与y轴平行，x’轴垂直于y’o’z’平面与卫星飞行方向夹角为β。其坐标转换关系为z’＝z/cosβ，ξ’(ρ)＝ρxtanβ+ξ(ρ)/cosβ。而此时由于入射角度发生变换，陆地和冰盖表面的朗伯体反射率β_r’＝cos(β+S)。因此，非天顶入射时探测器输出光子数N满足：

$N=\frac{A_R{T_a}^{2}Q{\mathrm{\&beta;}}_r\mathrm{\&eta;}{\cos }^2\mathrm{\&beta;}\cos (\mathrm{\&beta;}+S)}{\mathrm{\&pi;h\&upsi;}{z}^2}---\left(6\right)$

2.影响距离精度的理论模型

由于接收回波时间重心与光速乘积的一半即激光测距值，则时间重心标准差与光速乘积的一半则表示距离精度。

2.1目标特性影响

根据激光测高系统回波信号理论模型，利用概率统计方法可以得到由目标特性参数引起的回波信号时间重心的方差，Gardner已经建立包含目标粗糙度、倾斜度和在目标是斜面时指向角对测距精度影响的解析式，如(7)-(9)式所示，其中Var(T_pr)、Var(T_ps)和Var(T_pp)分别表示粗糙度、倾斜度和指向角影响时间重心方差的表达式。

$\mathrm{Var}\left(T_{\mathrm{pr}}\right)=(\frac{F}{N}+\frac{1}{K_s})\frac{4\mathrm{Var}\left(\mathrm{\&Delta;\&xi;}\right){\cos }^2{S}_{//}}{c^2{\cos }^2(\mathrm{\&beta;}+S_{//})}---\left(7\right)$

$\mathrm{Var}\left(T_{\mathrm{ps}}\right)=(\frac{F}{N}+\frac{1}{2K_s})\frac{4z^2{\tan }^2{\mathrm{\&theta;}}_T}{c^2{\cos }^2\mathrm{\&beta;}}[{\tan }^2(\mathrm{\&beta;}+S_{//})+\frac{{\tan }^2{S}_{\&perp;}{\cos }^2{S}_{//}}{{\cos }^2(\mathrm{\&beta;}+S_{//})}]---\left(8\right)$

$\mathrm{Var}\left(T_{\mathrm{pp}}\right)=\frac{4z^2{(1+{\tan }^2{\mathrm{\&theta;}}_T)}^2}{c^2{\cos }^2\mathrm{\&beta;}}[{\tan }^2(\mathrm{\&beta;}+S_{//})\mathrm{Var}\left(\mathrm{\&Delta;}{\mathrm{\&beta;}}_{//}\right)+\frac{{\tan }^2{S}_{\&perp;}{\cos }^2{S}_{//}{\cos }^2\mathrm{\&beta;}}{{\cos }^2(\mathrm{\&beta;}+S_{//})}\mathrm{Var}\left(\mathrm{\&Delta;}{\mathrm{\&beta;}}_{\&perp;}\right)]$

(9)

式中N为回波光子数；F为APD探测器噪声系数，通常在3～7之间；K_S为接收孔径面积与散斑相关面积的比值，K_S＝πA_R(2tanθ_T/λ)²；Std(Δξ)为目标表面粗糙度；β为光束指向角；z为卫星轨道高度；c为光速；θ_T为光束发散角；S_//和S_⊥分别表示沿轨和垂轨方向地表斜率，Var(Δβ_//)和Var(Δβ_⊥)分别表示沿轨和垂轨方向光束指向角误差。通过(7)-(9)式可以看出，地表斜率影响所有三项目标特性误差，然而上述表达式过于复杂，如果假设地表斜率为一维分布，且与指向角、指向角误差方向都相同，如图3所示，即S＝S_//，S_⊥＝0；Δβ＝Δβ_//，Δβ_⊥＝0；在地表为斜面时指向角误差将导致额外的测距误差ΔR_pp，而如果指向角很小或地表为平面时，这项影响可以忽略。上述三项对距离影响为ΔR＝std(T)*c/2；事实上，对于发散角只有亚毫弧量级的激光测高系统，(9)式tan²θ_T影响很小，可以忽略，则由地表参数影响的测距误差可以简化为：

$\mathrm{\&Delta;}{R}_{\mathrm{pr}}=\sqrt{(\frac{F}{N}+\frac{1}{K_s})}\frac{\mathrm{Std}\left(\mathrm{\&Delta;\&xi;}\right)\cos S}{\cos (\mathrm{\&beta;}+S)}---\left(10\right)$

$\mathrm{\&Delta;}{R}_{\mathrm{ps}}=\sqrt{(\frac{F}{N}+\frac{1}{2K_s})}\frac{z\tan {\mathrm{\&theta;}}_T}{\cos \mathrm{\&beta;}}\tan (\mathrm{\&beta;}+S)---\left(11\right)$

$\mathrm{\&Delta;}{R}_{\mathrm{pp}}=\frac{z\tan (\mathrm{\&beta;}+S)\mathrm{\&Delta;\&beta;}}{\cos \mathrm{\&beta;}}---\left(12\right)$

2.2器件噪声影响

由测高系统硬件噪声影响导致测距误差可以表示为ΔR_n＝Std(T_n)*c/2，表达式为：

$\mathrm{\&Delta;}{R}_n=\frac{c}{2}\sqrt{(\frac{N_Q{T}_r^3}{12N^2\mathrm{\&Delta;t}}+\frac{\mathrm{\&Delta;}{t}^2}{12})+[\frac{{\mathrm{FN}}_B{T}_r^3}{12N^2\mathrm{\&Delta;t}}+\frac{({\mathrm{FN}}_D+N_T+N_A)T_r^3}{12N^2\mathrm{\&Delta;t}}]}---\left(13\right)$

式中T_r为回波信号的全脉宽，通常取4倍回波信号脉宽σ_S，以保证回波信号能够被完整采集；Δt表示数字化仪的采样间隔；N_B、N_D、N_T、N_A和N_Q分别表示辐射背景噪声、探测器暗电流噪声、热噪声、前置放大噪声和量化噪声，具体表达式满足：

$N_B=\frac{2{\mathrm{\&eta;}}_d\mathrm{B\&Delta;}{t}^2}{\mathrm{hv}}\&CenterDot;P_b=\frac{2{\mathrm{\&eta;}}_d\mathrm{B\&Delta;}{t}^2}{\mathrm{hv}}\&CenterDot;\frac{S\left(\mathrm{\&lambda;}\right)\mathrm{\&Delta;\&lambda;}{\mathrm{\&theta;}}_{\mathrm{FOV}}^2{A}_R{\mathrm{\&beta;}}_r{T}_a^2{\mathrm{\&eta;}}_r}{4}---\left(14\right)$

$N_D=\frac{2B{I}_{\mathrm{dk}}\mathrm{\&Delta;}{t}^2}{q}---\left(15\right)$

$N_T=\frac{4\mathrm{KTFB}}{R_L}\&CenterDot;\frac{\mathrm{\&Delta;}{t}^2}{G^2{q}^2}---\left(16\right)$

$N_A=\frac{I_a^2\mathrm{B\&Delta;}{t}^2}{G^2{q}^2}---\left(17\right)$

$N_Q=\frac{A^2\mathrm{\&Delta;}{t}^2}{12G^2{q}^2}---\left(18\right)$

其中，B为硬件带通滤波器带宽，η_d探测器量子效率，h为普朗克常数，v为光波频率，P_b为背景噪声平均功率，S(λ)为激光波长波段的太阳光辐射亮度，Δλ为接收滤光片带宽，θ_FOV为接收望远镜视场角，A_R为接收望远镜面积，β_r为光斑内目标反射率，T_a为大气透过率，η_r为接收光学系统透过率，I_dk为探测器暗电流，R_L为等效负载电阻，T为探测器等效噪声温度，K为波尔兹曼常数，A为A/D转换器量化步长，I_a为前置放大器等效噪声电流谱密度，q为电子电量，G为探测器增益。

器件噪声影响主要表现在接收回波波形信号方面，通过相应的信号处理算法处理后，由器件影响的距离误差较小，对于GLAS系统，起主导作用为采样噪声，约4cm左右；对于固体地表目标，由大气和潮汐等环境因素影响的距离误差经过相应算法修正后残差约为2-3cm；而地表斜率和粗糙度这类没有先验知识的环境误差则没有消除方法，需要根据上述误差模型分析在不同目标特性参数时影响程度，特别的，当地表斜率1°和指向角误差1.5″时，则三项影响分别为1cm，1.5cm和7.5cm，整体RSS距离精度 $\mathrm{\&Delta;R}=\sqrt{\mathrm{\&Delta;}{R}_e^2+\mathrm{\&Delta;}{R}_n^2+\mathrm{\&Delta;}{R}_{\mathrm{pr}}^2+\mathrm{\&Delta;}{R}_{\mathrm{ps}}^2+\mathrm{\&Delta;}{R}_{\mathrm{pp}}^2},$ 约9cm。

3.影响高程精度的理论模型

在星载激光测高系统测量参考坐标系中，如图4所示，激光发射参考点为原点O，Z轴指向地心方向，X轴指向卫星的飞行方向，Y轴垂直于卫星轨道平面，XYZ-O构成右手系，其激光脚点坐标[X,Y,Z]^T满足：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_O\\ Y_O\\ Z_O\end{array}\right]+M_{\mathrm{Altitude}}\left[\begin{array}{c}0\\ -R\sin \mathrm{\&beta;}\\ R\cos \mathrm{\&beta;}\end{array}\right]+\left[\begin{array}{c}\mathrm{\&Delta;}{X}_{\mathrm{SC}}\\ \mathrm{\&Delta;}{Y}_{\mathrm{SC}}\\ \mathrm{\&Delta;}{Z}_{\mathrm{SC}}\end{array}\right]-\left[\begin{array}{c}\mathrm{\&Delta;}{X}_O\\ \mathrm{\&Delta;}{Y}_O\\ \mathrm{\&Delta;}{Z}_O\end{array}\right]---\left(19\right)$

式中向量[X₀,Y₀,Z₀]^T为由GPS系统确定的卫星质心坐标，M_Altitude表示卫星飞行姿态矩阵，满足(20)式，[0,-Rsinβ,Rcosβ]^T表示激光指向矢量，第三项表示激光参考点与卫星质心的固定偏移量，第四项表示星上GPS相位中心与卫星质心固定偏移量，最后两项是常数向量。(20)式中ω、φ和κ分别代表卫星飞行中Z轴偏航角(Yaw)、Y轴俯仰角(Pitch)和X轴横滚角(Roll)；而实施例中ICESat卫星飞行姿态角在绝大多数情况趋近于零，只在需要变换参考轨道(8天或91天重复轨道)或者过境标定场时进行姿态机动进行调整。

根据(19)式将X、Y和Z按泰勒公式展开并保留至一次项，且假定在各项误差影响相互独立的情况下，根据误差传播定律和误差合成理论，推导得出X、Y和Z三个方向的精度表达式：

式中，Δω、Δφ和Δκ分别表示卫星平台的偏航、俯仰和横滚角误差；ΔX_O、ΔY_O和ΔZ_O分别表示卫星平台GPS X、Y和Z方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；ΔX、ΔY和ΔZ则表示三个方向精度。模拟数据表明当系统进行正常测量工作时，即三个姿态角趋近于零和1°指向角时，按照GLAS系统1″偏航、横滚和俯仰角误差、1.5″指向角误差、5cmGPS卫星定位误差，以及前文所述的典型地表约9cm距离误差等影响脚点定位精度因素，使用该模型计算激光脚点定位精度，其平面精度约5m，高程精度约15cm，在本步骤中，以激光发射参考中心建立了直角坐标系，其中Z轴指向天顶(垂直向上)，X、Y轴分别指向平行卫星前进方向(沿轨方向)和垂直于前进方向(垂轨方向)，XOY平面平行于地表，平面精度指光斑照射中心地表位置的偏移大小，高程精度指垂直Z方向的误差大小。

4.地表目标参数求解

以激光光轴与地表目标交点为原点O_D，X_D轴与激光测高系统测量参考坐标系X轴同向，Y_D和Z_D与Y和Z轴反向，建立被照射地表目标坐标系。沿X_D和Y_D均匀采样地表目标相对于X_DO_DY_D平面高度，即前文中ξ(ρ)，其中ρ²＝x²+y²。激光光斑照射区域直径D＝zθ_T，对于GLAS系统其光斑直径D＝66m；利用(x-X_D)²+(x-X_D)²≤D²/4原则筛选光斑直径范围内所有高程点(x_n,y_n,z_n)。

利用最小二乘法(Leverberg-Marquardt)对光斑内所有高程点(x_n,y_n,z_n)拟合最接近的平面方程z＝a₁x+a₂y+a₃，求出平面参数a₁、a₂和a₃。LM算法兼具梯度下降法的全局特性和高斯-牛顿法的局部收敛特性，是较常用的非线性拟合方法，其目标是使拟合数据与实际数据之间的余量尽量接近白噪声。所拟合平面的法向量为(a₁,a₂,-1)，X_DO_DY_D平面法向量为(0,0,1)，利用余弦角公式可得拟合平面斜率S表达式，单位为°。

$S=180-\arccos \frac{-1}{\sqrt{{a_1}^2+{a_2}^2+1}}---\left(24\right)$

光斑内目标表面粗糙度Std(Δξ)计算方法为直接求解所有的z_n标准差。

5.实际精度计算流程

由于在计算高程精度之前，必须先计算测距精度，而测距精度又包含由地表目标、环境和硬件器件影响因素，因此具体步骤为：

A.代入公式(13)中所需测高系统硬件参数，计算辐射背景噪声N_B、探测器暗电流噪声N_D、热噪声N_T、前置放大噪声N_A和量化噪声N_Q；对于如GLAS系统有回波波形记录功能和波形处理算法的测高仪这五项影响几乎可以忽略，主要受采样噪声影响，可以直接使用ΔR_n＝3^1/2cΔt/12简化计算；

B.输入地表目标高度模型数据，即沿X_D和Y_D均匀采样地表目标相对于X_DO_DY_D平面高度。利用(x-X_D)²+(x-X_D)²≤D²/4原则筛选光斑直径范围内所有高程点(x_n,y_n,z_n)，使用LM算法拟合平面参数并计算目标表面斜率S；计算z_n标准差得到表面粗糙度Std(Δξ)；

C.设定大气透过率T_a，对1064nm激光波长其通常值为0.7；设定目标表面反射率β_r，其固体地表通常值为0.3；以及其他测高系统硬件参数，利用(6)式计算回波光子数N；

D.使用B和C过程所计算表面斜率S和粗糙度Std(Δξ)，以及回波光子数N，代入测高系统硬件参数F，光束指向角β，使用(10)-(12)计算目标影响导致的测距误差ΔR_pr、ΔR_ps和ΔR_pp。

E.对地观测激光测高系统受环境影响包含如光束穿越大气产生的散射和折射，由章动和岁差等产生的固体潮汐等环境误差，目前经过相应的大气延迟修正算法和固体潮汐修正算法修正，环境影响误差修正后的残差ΔR_e约2-3cm；代入ΔR_n、ΔR_e、ΔR_pr、ΔR_ps和ΔR_pp利用RSS方法计算整体测距误差ΔR；

F.代入当系统进行正常测量工作时的三个方向姿态角ω、φ和κ，通常三个数值趋近于0；代入测高系统偏航、横滚和俯仰角误差Δω、Δφ和Δκ；指向角误差Δβ；三个方向GPS卫星定位误差ΔX_O、ΔY_O和ΔZ_O；以及计算得出的距离误差ΔR；根据(21)-(23)式计算激光测高系统平面X和Y方向，以及高程Z方向精度。

在图5中，本发明使用了丹麦科技大学(DTU,Danmarks TekniskeUniversitet)提供的于2004年5月25日ALS(机载激光扫描测量系统,airborne laser scanning)对格陵兰岛沿岸附近冰盖和其以北的北冰洋海域海冰其中一段20km长轨迹的测量数据。其测量轨迹地理位置范围为纬度[83.21°,83.34°]，经度[311.84°,312.67°]，与ICESat卫星L2C工作周期Track0031星下点轨迹重合，数据平面精度约1m，高程精度约5cm，沿轨和垂轨方向点密度几乎均匀(约1.2m)；因此，可以使用ALS数据生成GLAS系统激光光斑内的DEM图，并对GLAS高程数据进行精度分析。由于GLAS高程数据基于TOPEX椭球，需转换至WGS-84椭球坐标；每个GLAS脚点直径70m光斑内包含约2000个ALS激光雷达数据。

在图8中，测量时指向角都为0.319°；ALS高程为GLAS脚点内2000个ALS测量点的平均高程，偏差为GLAS与ALS平均高程之差；上下两栏数据分别对应图6和图7，其粗糙度和斜率计算结果很好的反应了图中DEM的真实情况。利用所得脚点内斜率和粗糙度信息，以及GLAS系统参数，代入前文中(10-13)式计算各项距离误差影响大小，并使用RSS方法计算整体测距精度；最后将距离误差和姿态、位置误差代入(23)式高程精度模型，计算的精度显示在最后一列，模型精度与偏差结果有较好的一致性。

在图9中，统计全部106个GLAS脚点内斜率和粗糙度信息，所有脚点斜率在0-2°之间，粗糙度在0-1m之间，并做精度计算和偏差统计。统计结果显示ALS高程数据与GLAS高程数据绝对值平均偏差13.5cm，均方根偏差15.8cm，与典型地表15cm的高程精度有很好的一致性，高程偏差均在模型精度3σ以内，很好的符合了本发明所建立的高程精度模型。

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

Claims (6)

1.一种星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，包括以下步骤：

步骤1，根据测高系统硬件参数计算受器件影响的距离误差ΔR_n；并设定大气透过率T_a、目标表面反射率β_r、以及其它测高系统参数如发射能量Q，接收望远镜面积A_R；

步骤2，输入地表目标高度模型数据，即沿X_D和Y_D均匀采样地表目标相对于X_DO_DY_D平面的高度；筛选得到光斑直径D范围内的所有高程点(x_n,y_n,z_n)；使用LM算法拟合筛选后所有高程点(x_n,y_n,z_n)的最接近平面，并计算平面法向量(a₁,a₂,-1)，进而计算目标表面斜率S；并通过直接求解所有的z_n标准差得到筛选后的光斑内目标表面粗糙度Std(Δξ)；最后通过步骤1设定的参数计算回波光子数N；

步骤3：使用步骤2得到的表面斜率S和粗糙度Std(Δξ)，回波光子数N，并设定测高系统探测器噪声系数F，以及光束指向角β，即天底方向与激光光轴夹角，计算目标影响导致的测距误差ΔR_pr、ΔR_ps和ΔR_pp；

步骤4：基于一个受环境影响的修正残差值ΔR_e(是否需要稍微展开论述一下ΔR_e)，步骤1中得到的距离误差ΔR_n，以及步骤3中得到的测距误差ΔR_pr、ΔR_ps和ΔR_pp，根据误差合成理论，利用RSS方法(Root-Sum-Squares，平方和的根)计算由ΔR_e、ΔR_n、ΔR_pr、ΔR_ps和ΔR_pp组成的整体测距误差ΔR；其中，ΔR_e通常取2-3cm；

步骤5：根据系统进行正常测量工作时三个姿态角ω、φ和κ；测高系统偏航Δω、横滚Δφ和俯仰角误差Δκ；指向角误差Δβ；三个方向GPS卫星定位误差ΔX_o、ΔY_o和ΔZ_o；以及步骤5得出的整体测距误差ΔR；得到激光测高系统平面XY和高程Z方向精度。

2.根据权利要求1所述的星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，所述步骤1中，距离误差ΔR_n是基于以下公式得到：

$\mathrm{\&Delta;}{R}_n=\frac{c}{2}\sqrt{(\frac{N_Q{T}_r^3}{12N^2\mathrm{\&Delta;t}}+\frac{{\mathrm{\&Delta;t}}^2}{12})+[\frac{{\mathrm{FN}}_B{T}_r^3}{12N^2\mathrm{\&Delta;t}}+\frac{({\mathrm{FN}}_D+N_T+N_A)T_r^3}{12N^2\mathrm{\&Delta;t}}]},$

其中，c为真空中光速，T_r为回波信号的全脉宽，通常取4倍回波信号脉宽σ_S，N_B、N_D、N_T、N_A和N_Q分别表示辐射背景噪声、探测器暗电流噪声、热噪声、前置放大噪声和量化噪声，F为探测器噪声系数，N为回波光子数。

3.根据权利要求1所述的星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，所述步骤2中，筛选光斑直径D范围内所有高程点(x_n,y_n,z_n)是基于(x-X_D)²+(x-X_D)²≤D²/4原则来进行，所述目标表面斜率S是基于以下公式得到：

其中，a₁和a₂分别是步骤2中使用LM算法拟合所得平面法向量(a₁,a₂,-1)的前两项；

所述回波光子数N是基于以下公式得到：

其中，A_R为接收望远镜面积，T_a为大气透过率，Q为激光发射能量，β_r为目标表面反射率，η为接收系统能量利用率，β为光束指向角，h为普朗克常数，v为光波频率，S为目标表面斜率。

4.根据权利要求1所述的星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，所述步骤3中，测距误差ΔR_pr、ΔR_ps和ΔR_pp，是基于以下公式得到：

其中，N为回波光子数；F为探测器噪声系数，通常在3～7之间；K_S为接收孔径面积与散斑相关面积的比值，K_S＝πA_R(2tanθ_T/λ)²，A_R为接收望远镜面积，θ_T为光束发散角，λ为激光波长；Std(Δξ)为目标表面粗糙度；S为目标表面斜率；β为光束指向角；

${\mathrm{\&Delta;R}}_{\mathrm{ps}}=\sqrt{(\frac{F}{N}+\frac{1}{2K_s})}\frac{z\tan {\mathrm{\&theta;}}_T}{\cos \mathrm{\&beta;}}\tan (\mathrm{\&beta;}+S),$ 其中，N为回波光子数；F为探测器噪声系数；K_S为接收孔径面积与散斑相关面积的比值；z为卫星轨道高度；θ_T为光束发散角；S为目标表面斜率；β为光束指向角；

其中，z为卫星轨道高度；S为目标表面斜率；β为光束指向角；Δβ为光束指向角误差；

ΔR_pr表示由地表粗糙度引起的测距误差；

ΔR_ps表示由地表斜率引起的测距误差；

ΔR_pp表示当地表斜率不为零时由指向角引起的额外距离误差。

5.根据权利要求1所述的星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，所述步骤4中，使用RSS方法计算整体测距误差ΔR，即计算ΔR_e、ΔR_n、ΔR_pr、ΔR_ps和ΔR_pp平方和的根，具体公式为：

$\mathrm{\&Delta;R}=\sqrt{\mathrm{\&Delta;}{R}_e^2+\mathrm{\&Delta;}{R}_n^2+\mathrm{\&Delta;}{R}_{\mathrm{pr}}^2+\mathrm{\&Delta;}{R}_{\mathrm{ps}}^2+\mathrm{\&Delta;}{R}_{\mathrm{pp}}^2}.$

6.根据权利要求1所述的星载激光测高系统固体地表目标平面和高程精度检测方法，其特征在于，所述步骤5中，测高系统的平面精度和高程精度基于以下公式获得：

平面精度的X方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δω、Δφ和Δκ分别表示卫星平台的偏航、俯仰和横滚角误差；ΔX_O表示卫星平台X方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值；

平面精度的Y方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δω、Δφ和Δκ分别表示卫星平台的偏航、俯仰和横滚角误差；ΔY_O表示卫星平台Y方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值；

高程精度的Z方向基于以下公式：

其中，ω、φ和κ分别代表卫星飞行中Z轴偏航角、Y轴俯仰角和X轴横滚角，Δφ和Δκ分别表示卫星平台俯仰和横滚角误差；ΔZ_O表示卫星平台Z方向定位误差；Δβ表示光束指向角误差；ΔR表示整体测距误差；R表示激光测距值。