CN101915912A - 一种全面的激光测高回波仿真方法 - Google Patents

Abstract

本发明涉及遥感探测技术领域，尤其涉及一种全面的激光测高回波仿真方法。本发明不仅能够精确仿真倾斜或者起伏的裸地面所产生的回波信号，而且能够精确仿真有植被覆盖的地面所产生的回波信号。在此基础之上，本发明还能够通过所提出的植被覆盖有效面积比概念指导仿真实验设计，从而推导地面起伏或者倾斜、植被厚度和植被覆盖有效面积等因素变化，对回波所产生的影响，深入探索如何从真实回波信号反演地形和地物参数的方法。最后，本发明还能够考虑激光指向角不是指向天底方向，而与天底方向有一定的夹角时所产生的回波。

Description

一种全面的激光测高回波仿真方法

技术领域

本发明涉及遥感探测技术领域，尤其涉及一种全面的激光测高回波仿真方法。

背景技术

天基遥感探测技术是地球科学领域的重要技术，其中天基激光测高是利用卫星搭载激光测高仪、从太空对星球表面轮廓进行探测。天基激光测高的工作原理为，天基激光器持续地向地面发射激光，激光光束穿越大气到达地面后产生微弱的后向散射回波，该散射回波再次穿越大气被测高系统中的望远镜所接收；通过光电探测器将激光回波转换成电脉冲回波，通过分析该电脉冲回波的渡越时间TOF(time of flight)得到激光脉冲的渡越时间，从而获得地表足迹光斑与测高仪间的距离R_m；然后通过下式计算激光光斑处的地表高度：

$h_s={[R_s^2+R_m^2-{2R}_s{R}_m\cos \mathrm{\φ}]}^{1/2}-R_{\mathrm{ref}}$

式中，R_s是卫星轨道半径，φ是激光光轴与卫星至地心连线的夹角(即天底偏离角)，R_ref是地球参考表面的半径，通常为大地水准面；激光测高仪的联结方程(link equation)如下：

$E_r=E_t\frac{A_r}{R_m^2}\frac{r_s}{\mathrm{\π}}{\mathrm{\τ}}_a^2{\mathrm{\τ}}_s$

式中，E_r是回波脉冲能量，E_t是发射脉冲能量，A_r是接收望远镜的面积(平方米)，R_m是地表足迹光斑与测高仪间的距离，r_s是目标表面反射率(朗伯)，τ_a是单程的大气透过率，τ_a是系统透过率。显然，地表足迹光斑内的高度分布和地表反射率都会对激光回波信号产生影响。因此，在A_r、τ_a和τ_s已知的情况下，如果能够获得发射脉冲信号与回波脉冲信号，就能够通过数据处理与分析获取地表足迹光斑内的高度分布和反射率信息，从而为天基激光测高仪增加雷达探测功能。

渡越时间的测量精度除了受发射脉冲宽度和模拟/数字转换器ADC(Analog-to-digitalConverter)的采样效应(量化精度和采样间隔)影响之外，还要受到地理特征的影响。地理形状决定了脉冲的拓宽以及回波的整个形状，地表起伏也会引入进一步的波形拓宽，地表反射率决定了回波能量。因此，激光测高仪的回波波形包含了地表结构、起伏和反射率等信息，如何从回波波形中提取这些信息是相关领域研究者所普遍关注的问题。

在没有噪声，被测表面又具有相同的反射率的情况下，回波波形只取决于发射波形和表面几何形状。在简单情况下，如平面或均匀斜坡面，可以推导出描述回波形状的闭合表达式。但在复杂的表面地形下，解析推导是无法进行的。数值模拟有助于分析这些复杂的因素，研究测高仪设计参数、测高性能与地理特征之间的关系，理解复杂表面特征对回波形状的影响，提高测距精度分析。对激光测距的精度分析以及提取回波所蕴含信息的需要推动了激光测高仪仿真软件的研制。仿真软件的用途在于：研究激光测高仪各个组成子系统的设计参数、被测表面复杂的地形特征与回波信号之间的关系；分析测距精度以及各种因素对测距精度的影响。

目前，世界主要发达国家都在大力发展天基激光测高技术，美国明显处于世界领先地位。受美国海军和导弹防御组织支持的位于著名的麻省理工大学的林肯实验室，为美国军方研制了一系列的军用激光测高系统。而由美国国家航天航空局NASA(National Aeronautics andSpace Administration)组织研制的一系列已经出色完成太空任务的民用天基激光测高系统，同样能够反映这一研究领域的世界最新水平。其中最具代表性的是美国地球观测系统EOS(EarthObserving System)中的地球科学激光测高系统GLAS(Geoscience Laser Altimeter System)。这是第一个用于连续全球观测的激光测高系统，能够测量冰被地形和相应的温度变化，同时也监测云层和大气的特性。此外，由NASA地球科学探路者计划资助研制出的植被覆盖雷达VCL(Vegetation Canopy Lidar)能够通过测量植被的垂直和水平结构以及土壤表面地形来描述地球的三维结构。航天激光测高系统SLA(Shuttle Laser Altimeter)计划受到由NASA总部和Goddard空间飞行中心GSFC(Goddard Space Flight Center)发起的行星地球任务计划的支持。由于采用了高速的模数采集电路，该装置能够分析光斑范围内的地表高度变化。于1996年11月7号升空的火星勘探号MGS(Mars Global Surveyor)宇宙飞船所携带的激光测高系统MOLA(Mars Orbiter LaserAltimeter)-2的距离分辨率为37cm，能够以约300m的间距探测火星表面轮廓。

发明内容

本发明的目的是提供一种全面的激光测高回波仿真方法，在已知激光测高仪的硬件参数、激光测高仪的飞行参数、大气与地表的光学参数、地面高度模型的条件下，计算出激光测高仪的足迹光斑的回波波形。

为达到上述目的，本发明采用如下的技术方案：

输入激光测高仪的硬件参数；

输入激光测高仪的飞行参数；

输入大气与地表的光学参数；

输入地表高度模型数据：地表高度模型用矩形点集(x，y，t，cod)来表示，其中x和y为网格点的横坐标与纵坐标，表征网格点在与天底方向垂直平面上的位置，t为渡越时间，cod为网格点的高度模型类型；

对地面高度数据进行初步的处理，剔除光斑以外的数据，保留光斑以内的数据；

在xoy平面组成三角形网格，确定每个三角单元的高度模型类型；

根据发射脉冲功率面密度的空间分布、三角形单元在光斑内的位置和高度模型类型，计算每个三角单元的功率面密度；

选择计算回波合适的时间单元(time bin)，根据高度模型类型，确定每个三角单元渡越时间的分布范围，根据每个时间单元所占据的面积及其功率面密度，计算每个时间单元的返回光子数目；

将激光光斑内每个三角单元的回波信号叠加起来，得到整个光斑内的回波信号H(n)；

根据激光器发射脉冲功率的时间分布I(n)，利用卷积和求系统对激励I(n)的响应Y(n)：Y(n)＝I(n)*H(n)；

利用LM方法对有植被覆盖地面的仿真回波进行高斯拟合，获得最后一个高斯分量之外的每个高斯分量的能量与总回波能量的比，将它定义为该层植被在地表足迹光斑内所占据的有效面积比。

所述激光测高仪的硬件参数包括激光波长、激光发散角、激光发射能量、激光脉冲宽度、望远镜接收面积和系统透过率。

所述激光测高仪的飞行参数包括轨道高度R_h和天底偏离角φ。

所述大气与地表的光学参数包括大气透过率和地表反射率。

步骤(5)所采用的判据为(x-x₀)²+(y-y₀)²≤a²，其中a为光斑半径，这里取地表足迹光斑的3倍空间均方根宽度。

本发明具有以下优点和积极效果：

1)不仅能够仿真高度分布连续的裸地面的回波信号，也能够仿真高度分布断续的植被的回波信号；

2)采用规范化的矩形网格给出地面高度模型，能够非常方便地进行三角单元的划分；

3)提出的有效面积比概念能够指导仿真实验的设计，深入研究植被覆盖有效面积比变化对回波所产生的影响。

4)能够考虑激光器指向偏离天底方向对回波的影响。

附图说明

图1是本发明提出的一种全面的激光测高回波仿真方法的流程图。

图2是三角单元高度模型类型的确定示意图。

图3是飞行器与地表足迹光斑网格点的距离计算示意图。

图4是连续类型三角单元的时间单元划分示意图。

具体实施方式

下面以具体实施例结合附图对本发明作进一步说明：

本发明提供的一种全面的激光测高回波仿真方法，如图1所示，采用如下步骤：

(1)输入激光测高仪的硬件参数；

激光测高仪的硬件参数包括激光波长、激光发散角、激光发射能量、激光脉冲宽度、望远镜接收面积和系统透过率。例如，对于航天测高仪Shuttle Laser Altimeter II(SLA-02)，相应的硬件参数如下：激光波长为1064nm；激光2倍均方根发散角为0.25mrad；激光发射单脉冲能量为35mJ；激光脉冲半高全宽为15nS；望远镜接收面积为πr²＝0.113m²；系统透过率为50％。

(2)输入激光测高仪的飞行参数。

激光测高仪的飞行参数包括轨道高度R_h和天底偏离角φ。例如，对于航天测高仪SLA-02，相应的飞行参数如下：轨道高度为281.945-291.999km；天底偏离角是变化的，每个激光脉冲对应于一个特定的天底偏离角，对于SLA-02的第一次观测，天底偏离角在0.020-1.708度之间变化。

(3)输入大气与地表的光学参数。

大气与地表的光学参数包括大气透过率和地表反射率。例如，对于航天测高仪SLA-02，大气透过率和地表反射率都是变化的，每个激光脉冲对应于一个特定的大气透过率和地表反射率。

(4)输入地表高度模型数据：地表高度模型用矩形点集(x，y，t，cod)来表示，其中x和y为网格点的横坐标与纵坐标，表征网格点在与天底方向垂直平面上的位置，t为渡越时间，cod为网格点的高度模型类型。

地表高度模型数据采用一个四维点集来表示，高度模型类型只有两种——连续和间断：对于裸地面，高度模型类型为连续，用数字1表示；对于植被，高度模型类型为断续，用数字0表示。

(5)对地面高度数据进行初步的处理——剔除光斑以外的数据，保留光斑以内的数据。

所采用的判据为(x-x₀)²+(y-y₀)²≤a²，其中a为光斑半径，这里取地表足迹光斑的3倍空间均方根宽度。

(6)在xoy平面组成三角形网格，确定每个三角单元的高度模型类型。

每个三角单元用三个点集A、B和C来表示。将离内切圆圆心距离最近的三角形顶点处的高度模型类型作为该三角单元的高度模型类型。如图2所示，选择顶点B的高度模型类型作为该三角单元的高度模型类型。

每个点采用四维坐标(x，y，t，cod)来表示，渡越时间t是根据飞行器与地表足迹光斑网格点的距离R_m计算出来的：

其中光速c＝299,792,458ms^-1。而飞行器与地表足迹光斑网格点的距离R_m则是根据激光测高仪轨道高度R_h、地表足迹光斑网格点的高度h_s、偏离天底角φ和地球参考表面的半径R_ref(通常为大地水准面)计算出来的，如图3中所示；

$R_m^2-{2R}_s{R}_m\cos \mathrm{\φ}+R_s^2={(h_s+R_{\mathrm{ref}})}^2$

式中，R_s＝R_h+R_ref。这里考虑了激光指向角与天底方向的偏离。

(7)根据发射脉冲功率面密度的空间分布(xoy平面)、三角单元在光斑内的位置和高度模型类型，计算每个三角单元的功率面密度。

假定发射脉冲功率面密度的空间分布(xoy平面)也呈高斯型：

$E_A(x,y)=\frac{E_T}{{2\mathrm{\π\σ}}^2}{e}^{-\frac{{(x-x_0)}^2+{(y-y_0)}^2}{{2\mathrm{\σ}}^2}}$

E_T为发射能量，(x₀，y₀)为光斑中心点坐标，σ为发射脉冲的空间均方根宽度。根据三角单元在光斑内的位置和高度模型类型，计算每个三角单元的功率面密度：当高度模型类型为连续时，将每个三角单元的内切圆圆心处(如图2中的D点)的功率面密度作为整个三角形网格的功率面密度；当高度模型类型为间断时，将离内切圆圆心距离最近的三角形顶点处的功率面密度作为整个三角单元的功率面密度。如图2所示，选择顶点B处的功率面密度作为该三角单元的功率面密度。

(8)选择计算回波合适的时间单元长度；根据高度模型类型，确定每个三角单元渡越时间的分布范围；根据每个时间单元所占据的面积及其功率面密度，计算每个时间单元的返回光子数目。

首先选择计算回波合适的时间单元长度，也即time bin宽度。一般说来，时间单元的宽度越小，回波的计算就越精确，但同时计算量也就越大。通常选择发射脉冲的3倍均方根宽度的千分之一到百分之一之间。对于天基激光测高仪，选择0.1nS通常是合适的。随后根据三角单元高度模型类型，确定每个三角单元的渡越时间分布范围：当高度模型类型为连续时，比较三角单元三个顶点处的渡越时间t_A、t_B和t_C，找出最小值和最大值。这里假定最小值为t_A，最大值为t_C。以选定的时间单元为单位划分三角单元(图4)。每个三角单元内采用线性插值；计算三角单元内每个所对应的面积A_N，从而求得每个时间单元所返回的光子数目。当高度模型类型为断续时，将离内切圆圆心D距离最近的三角形顶点B处的渡越时间作为整个三角单元的渡越时间(如图2所示)。这样一来，整个三角单元的返回光子处于同一个时间单元，根据三角单元的面积计算该三角单元的返回光子数目。

(9)将激光光斑内每个三角单元的回波信号叠加起来，得到整个光斑内的回波信号H(n)。

(10)根据激光器发射脉冲功率的离散时间分布I(n)，利用卷积和求系统对激励I(n)的响应Y(n)：Y(n)＝I(n)*H(n)。

假定激光器发射脉冲功率的时间分布呈高斯型：

$I\left(t\right)=\frac{A}{\sqrt{2\mathrm{\π}}\mathrm{\δ}}{e}^{-\frac{{(t-t_0)}^2}{{2\mathrm{\δ}}^2}}$

式中，A为单脉冲能量，δ为发射脉冲的时间均方根宽度。然后利用卷积和求系统对I(n)的响应y(n)：y(n)＝I(n)*H(n)。I(n)的时间长度我们取3倍的均方根宽度。

(11)利用LM(Leverberg-Marquardt)方法对有植被覆盖地面的仿真回波进行高斯拟合，获得最后一个高斯分量之外的每个高斯分量的能量与总回波能量的比，将它定义为该层植被在地表足迹光斑内所占据的有效面积比。

LM方法是梯度下降法与高斯一牛顿法的结合，既具有高斯-牛顿法的局部收敛特性又具有梯度法的全局特性，而且LM方法利用了误差指标函数近似的二阶导数信息，所以运算速度比梯度法快得多。设w为拟合的参数向量，r(w)为拟合数据的余量向量，J(w)为r关于w的Jacobian矩阵，那么LM方法的形式为：

Δw＝-[J^T(w)J(w)+uI]^-1J(w)^Tr(w)

式中，u是比例系数，I是单位矩阵。从上式可看出，如果比例系数u＝0，则为高斯-牛顿法；如果u取值很大，则LM算法接近梯度下降法。每迭代成功一步，则u减小一些，这样在接近误差目标的时候，逐渐与高斯-牛顿法相似。高斯-牛顿法在接近误差的最小值的时候，计算速度更快，精度也更高。LM方法是地球物理学数据分析中常用的非线性拟合方法。

植被在地表足迹光斑内所占据的有效面积比不仅与植被在足迹光斑内所占据的面积有关，还与植被在光斑中的位置，也即与光斑中心的距离远近有关。

本发明提供的一种全面的激光测高回波仿真方法所具备的功能如下：

能够计算各种不同参数的激光测高仪所产生的回波波形；不仅能够计算高度分布连续的裸地面的回波，也能够计算高度分布断续的植被的回波；提出的植被覆盖有效面积比概念能够指导仿真实验的设计，从而研究其变化对回波所产生的影响；能够考虑指向角偏离天底方向对回波的影响。

本发明所涉及的激光测高回波仿真方法的最大特点是功能全面：不仅能够精确仿真倾斜或者起伏的裸地面所产生的回波信号，而且能够精确仿真有植被覆盖的地面所产生的回波信号。对于倾斜或者起伏的裸地面，其表面高度的分布通常是连续的；而对于有植被覆盖的地面，植被的高度变化通常是不连续的。因此在处理地表足迹光斑的数字地形数据时，必须充分考虑这一点，才能够得到好的仿真结果。其次，在本仿真方法中，还提出了植被和裸地面的有效面积比概念。我们定义植被和裸地面的有效面积比分别为由植被和地面单独产生的高斯信号的能量与回波总能量之比。有效面积比虽然不能给出植被在光斑中所占据的真实面积与光斑总面积之比，但却是植被在光斑中所占据的面积以及该面积在光斑中的位置的综合体现。本发明既能够精确仿真有植被覆盖的地面所产生的回波信号，又能够通过所提出的植被覆盖有效面积比概念指导仿真实验设计，从而推导地面起伏或者倾斜、植被厚度和植被覆盖有效面积等因素变化，对回波所产生的影响，深入探索如何从真实回波信号反演地形和地物参数的方法。最后，本发明还能够考虑激光指向角不是指向天底方向，而与天底方向有一定的夹角时所产生的回波。

Claims (5)

1.一种全面的激光测高回波仿真方法，其特征在于，包括以下步骤：

(1)输入激光测高仪的硬件参数；

(2)输入激光测高仪的飞行参数；

(3)输入大气与地表的光学参数；

(4)输入地表高度模型数据：地表高度模型用矩形点集(x，y，t，cod)来表示，其中x和y为网格点的横坐标与纵坐标，表征网格点在与天底方向垂直平面上的位置，t为渡越时间，cod为网格点的高度模型类型；

(5)对地面高度数据进行初步的处理，剔除光斑以外的数据，保留光斑以内的数据；

(6)在xoy平面组成三角形网格，确定每个三角单元的高度模型类型；

(7)根据发射脉冲功率面密度的空间分布、三角形单元在光斑内的位置和高度模型类型，计算每个三角单元的功率面密度；

(8)选择计算回波合适的时间单元长度，根据高度模型类型，确定每个三角单元渡越时间的分布范围，根据每个时间单元所占据的面积及其功率面密度，计算每个时间单元的返回光子数目；

(9)将激光光斑内每个三角单元的回波信号叠加起来，得到整个光斑内的回波信号H(n)；

(10)根据激光器发射脉冲功率的时间分布I(n)，利用卷积和求系统对激励I(n)的响应Y(n)：Y(n)＝I(n)*H(n)；

(11)利用LM方法对有植被覆盖地面的仿真回波进行高斯拟合，获得最后一个高斯分量之外的每个高斯分量的能量与总回波能量的比，将它定义为该层植被在地表足迹光斑内所占据的有效面积比。

2.根据权利要求1所述的全面的激光测高回波仿真方法，其特征在于：

所述激光测高仪的硬件参数包括激光波长、激光发散角、激光发射能量、激光脉冲宽度、望远镜接收面积和系统透过率。

3.根据权利要求1所述的全面的激光测高回波仿真方法，其特征在于：

所述激光测高仪的飞行参数包括轨道高度R_h和天底偏离角φ。

4.根据权利要求1所述的全面的激光测高回波仿真方法，其特征在于：

所述大气与地表的光学参数包括大气透过率和地表反射率。

5.根据权利要求1所述的全面的激光测高回波仿真方法，其特征在于：

步骤(5)所采用的判据为(x-x₀)²+(y-y₀)²≤a²，其中a为光斑半径，这里取地表足迹光斑的3倍空间均方根宽度。

Images