CN1959431A - 一种大脚印激光脉冲回波波形的仿真算法 - Google Patents

Abstract

一种大脚印激光脉冲回波波形仿真算法，包括步骤：1.结合星载激光高度计测高相关参数，推导出激光测高理论公式；2.考虑测高激光的时间和空间高斯特性，把发射激光脉冲分解为无数多个入射到大脚印上各点的高斯脉冲，经大脚印对应点发射后被激光高度计接收，得到对应激光脉冲回波功率的积分公式；3.结合大脚印地形特点，进行地形分割，得到激光脉冲回波功率的简单公式，编程实现得到大脚印激光脉冲回波的波形。本发明不仅可用于对太空卫星测高的理论进行仿真分析，为实际测高提供先验性的参考，而且仿真得到的输入地形脉冲回波对于纠正测高精度和反演三维地形有一定的参考价值。

Description

一种大脚印激光脉冲回波波形的仿真算法

技术领域：

本发明涉及航空探测中星载激光高度计测高仿真，特别是针对大脚印激光脉冲回波波形的一种仿真算法。

背景技术：

在航空航天科技中，通过激光测高得到探测目标的高度数据是必不可少的技术。比如中国的“嫦娥”探月工程中第一阶段的航天器“嫦娥一号”卫星的重要有效载荷之一激光高度计将利用主动激光测高技术，测量卫星到月球表面天底点的距离，得到的数据经过地面应用系统处理后可获得月球表面的高度数据，再结合拍摄的月球纹理照片获得月球三维立体地形等。这不仅进一步加深对月球信息的了解，而且为后期探月阶段中卫星的软着陆中着陆点的选取，提供了强有力的依据。

但是，限于实验条件的限制和目前国际在探月方面共享的信息很少，也不能轻易做一次卫星飞行试验来验证一些理论研究，对卫星飞行中测高回波波形的仿真就有助于预先了解一些特性，从理论上掌握一些先验知识。但是目前国内外对激光测高回波波形的仿真主要集中在机载激光测高方面，这种情况下激光高度计与被测目标间距离小，只有几公里，导致激光光斑(即激光脚印)小，一般面积不到1m²，而且容易试验，可以结合现实情况对其仿真模型进行验证和修改，所以小脚印激光脉冲回波的仿真比较成熟，而在大脚印方面，因为激光高度计距离探测天体一般有几百公里，对应的激光光斑面积大到上万平方米，而且试验成本高，在机载方面的理论研究不再适应大脚印远距离的情况，据调查，国内尚无这方面的专门研究，国外在这方面的研究也透露甚少。所以，本仿真算法对大脚印激光脉冲回波的仿真，不仅有利于对大脚印情况下的回波波形做一些理论研究，得到先验知识，而且仿真得到的相应地形对应的脉冲回波对于纠正测高精度和反演三维地形有一定的参考价值。

发明内容：

本发明的目的是提供一种大脚印激光脉冲回波波形的仿真算法，使输入一定的大脚印地形，通过该算法能得出理论上经过激光高度计测高其对应的脉冲回波，可以用于根据不同的典型地形(如裂缝、山峰、斜坡、丘陵等)对应的脉冲回波，建立地形-回波的相应数学模型，用于修正实际测高中测高精度的修正，和根据实际回波波形对其对应地形定性分析等。

本发明仿真的过程为，从激光高度计发射激光到目标地形，再经目标地形反射，激光回波被激光高度计接收，对于激光高度计的内部发射结构和接收结构不做考虑。

本发明的技术方案如下：

对于星载激光高度计，在探测过程中由于一次测高时间比较短，故可以近似为缓慢移动的收发合置系统，同时在太空中没有大气衰减的情况，探测目标为表面漫反射的大目标，故可以得到综合以上情况下的激光测高公式如下：

$P_r=\frac{\mathrm{\πf}{T_A}^2{D}^2{\mathrm{\η}}_t{\mathrm{\η}}_r}{{4r}^2}{P}_s---\left(1\right)$

其中，Pr为接收的回波功率，Ps为发射的激光功率，f为地形后向散射率，TA为大气传输系数，D为高度计的接收孔径，ηt为发射机的光学效率，ηr为接收机的光学效率，r为目标到激光高度计的距离。

由于星载高度计采用的激光脉冲特性是在时间和空间上都是高斯分布的，即在脉冲持续时间上脉冲功率随时间的的变化是一个高斯波形，在空间上发射激光脉冲随光束与中心轴线的夹角成高斯分布的。那么发射激光脉冲功率与时间的关系是：

$P_s\left(t\right)=\frac{A}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp (-\frac{{(t-\frac{\mathrm{\τ}}{2})}^2}{{2\mathrm{\σ}}^2})---\left(2\right)$

其中σ是标准方差，A为一个调整峰值大小的乘积因子，而τ为发射的激光脉冲半峰值功率脉宽。工程上选取峰值功率ps_max和半峰值功率点脉宽τ作为激光脉冲的相关参数，经过推导可以得到A，σ与ps_max，τ对应关系为：

$A=\frac{(\sqrt{\mathrm{\π}}\mathrm{\τps}\_\max )}{2\sqrt{\ln 2}},\mathrm{\σ}=\frac{\mathrm{\τ}}{2\sqrt{2\ln 2}}---\left(3\right)$

因为激光光束本身存在发散，当激光垂直入射到地面时，会在其上面形成一个圆形的激光光斑，其半径大小近似等于激光高度计到天底点的距离H和1/2束宽α的乘积。那么，激光入射到地形上，激光光斑上的每一点都会反射回激光回波，所以可以理解为高斯脉冲入射时被分成无数多个波形相同的高斯小脉冲，每一个脉冲对应入射到激光光斑上的一点，这些入射脉冲的能量应该满足高斯分布，以此来体现发射高斯脉冲在空间上的高斯分布特性。所以，返回的脉冲回波实际上就是这些无数个高斯脉冲经过反射后叠加形成的。故我们可以假设天底点入射的高斯激光脉冲为p_s0(t)，那么光斑内任意一点到激光高度计的路径与激光光轴夹角为φ的入射高斯脉冲为

          p_sxy＝p_s0(t)*exp(-2tan²/α²)            (4)

同时，根据能量守恒，在这个激光光斑内的这些入射高斯激光小脉冲积分起来就是入射的总激光脉冲p_s(t)，所以有下式：

上式即表示在激光光斑内对各点入射的这些高斯小脉冲积分，因为涉及到某一点对应的发散角，为了便于计算，采用极坐标系的积分方式，用(ρ，θ)来表示某一点。

那么显然有关系tan＝ρ/H，而这个激光光斑的半径为H*tanα≈H*α，所以可以把(5)变换为：

$p_s\left(t\right)={\∫}_0^{2\mathrm{\π}}{\∫}_0^{H\tan \mathrm{\α}}{p}_{s0}\left(t\right)*\exp (-{2\mathrm{\ρ}}^2/{(H*\mathrm{\α})}^2)\mathrm{\ρd\ρd\θ}---\left(6\right)$

经过简单的积分计算，结合(2)式，就可以得到天底点的入射高斯小脉冲为：

$p_{s0}\left(t\right)=\frac{\sqrt{2}A}{(1-\exp (-2)){\mathrm{\π}}^{3/2}\mathrm{\σ}{H}^2{\mathrm{\α}}^2}\exp (-\frac{{(t-\frac{\mathrm{\τ}}{2})}^2}{{2\mathrm{\σ}}^2})---\left(7\right)$

根据公式(1)有，瞬时的回波功率应该满足该方程，只是回波脉冲相对于发射脉冲来说，有一个2r/c的时间延迟，r为激光高度计到激光光斑内对应点的距离，所以结合(7)式，激光光斑内任一点(ρ，θ)的瞬时脉冲回波应该为：

同时，假设点(ρ，θ)相对于天底点的高度差为h(ρ，θ)，令

$k=\frac{\mathrm{\πf}{T_A}^2{D}^2{\mathrm{\η}}_t{\mathrm{\η}}_r}{4}*\frac{\sqrt{2}A}{(1-\exp (-2)){\mathrm{\π}}^{3/2}{\mathrm{\σH}}^2{\mathrm{\α}}^2}---\left(9\right)$

而该点到测高计的距离r为：

$r=\sqrt{{\mathrm{\ρ}}^2+{(H-h(\mathrm{\ρ},\mathrm{\θ})}^2}---\left(10\right)$

所以(8)式又可以表示为：

$p_r=(t,\mathrm{\ρ},\mathrm{\θ})=k*\exp (-\frac{2{\mathrm{\ρ}}^2}{H^2{\mathrm{\α}}^2})*\frac{\exp (-\frac{{(t-\frac{\mathrm{\τ}}{2}-\frac{2\sqrt{{\mathrm{\ρ}}^2+{(H-h(m,n))}^2}}{c})}^2}{{2\mathrm{\σ}}^2})}{{\mathrm{\ρ}}^2+{(H-h(\mathrm{\ρ},\mathrm{\θ}))}^2}---\left(11\right)$

那么整个激光光斑内产生的脉冲回波为：

因为实际获得的地形数据为离散的，而且上式虽然形式简单，但是积分很复杂，需要再对地形进行建模分析，并从中找到一种近似的算法，简单地计算出激光回波波形。

在本发明中，地形建模采取方法是：采用极坐标来描述大脚印地形的情况，并且用以天底点为圆心的一系列半径从

开始，相邻圆半径差距为

的

个同心圆和从天底点出发的一系列相邻夹角均为π/60的射线来分割地形，那么这

个同心圆和120条相邻夹角相等的射线将地面分成120 个小的网格。

光斑地形内任一点(ρ，θ)属于网格(m，n)的条件是：

$(n-1)\frac{\mathrm{\&pi;}}{60}\&le;\mathrm{\&theta;}<n\frac{\mathrm{\&pi;}}{60},(n=\mathrm{1,2},...,120)---\left(13\right)$

在这120 个小网格中，每个网格的面积不一样，靠近激光光斑边缘的最大的网格面积是靠近天底点的最小面积的几百倍。但是，这是合理的，因为当H取200000米，α取0.3mrad时，最大的网格面积才3m2多一点，在这么小的面积内其高度变化不大，可以认为在该网格内的高度是一定的，即近似成为平面。在π(Ha)²面积内划分120

个网格，把每个网格近似成一个平面，用120

个点来代替全部点的高度值，根据地形变化的特点来说，是完全可行的。当然，本仿真算法主要考虑到卫星对月球测高一般在200000米的情况，没有考虑更远的测高情况，选择了7200个网格来分割地形，当距离更远时，激光光斑的面积增大，选择7200个网格很可能不能足够准确的表现光斑内的地形高度特征，那么只需要把其中的地形分割更密集一些，得到更多的网格就可以了，算法本身不需做什么改动。

因此，由于高空探测时激光照射的面积比较大，在照射区域内的地形地貌、植被比较复杂，为了简化地面目标模型，假设如下：

1、把激光光斑内的地形分割为120 个网格，假设每个网格内的高度值一样，即该网格为平面；

2、激光光斑内各部分地形的后向散射率是一样的，即认为其是一个常数，实际上不同地形的后向反射率不是一个常数，其是激光参数(频率、极化方式、入射角)和地面参数(如湿度，粗糙度等)的复杂函数。但是相对于地球地形上分布着植被，冰川，海洋，山地等的复杂，在月球上地质情况相对比较单一，

所以这么假设是合理的。如果为了更加精确，可以加入对后向散射率的考虑。

根据(12)，结合对激光光斑地形的划分，可以得到：

其中，

$p_r=(t,\mathrm{\&rho;},\mathrm{\&theta;})=k*\exp (-\frac{2{\mathrm{\&rho;}}^2}{H^2{\mathrm{\&alpha;}}^2})*\frac{\exp (-\frac{{(t-\frac{\mathrm{\&tau;}}{2}-\frac{2\sqrt{{\mathrm{\&rho;}}^2+{(H-h(m,n))}^2}}{c})}^2}{{2\mathrm{\&sigma;}}^2})}{{\mathrm{\&rho;}}^2+{(H-h(\mathrm{\&rho;},\mathrm{\&theta;}))}^2}---\left(14\right)$

只要求出所有网格(m，n)内的回波，就可以求出整个地形的激光脉冲回波波形。根据前面的假设，在这个网格内选一点的高度值来代替整个网格的高度，相当于在该网格中h(m，n)是一个常数，可以设定一个60×120的地形采样高度矩阵h，矩阵中元素h(m，n)的值就是网格(m，n)的高度。

此外，在网格(m，n)中，因为网格内高度差不大，采取的是其中满足下述条件的一点(ρ_mn，θ_mn)，以它的高度值来代替整个网格的高度值h(m，n)。

${\mathrm{\&rho;}}_{\mathrm{mn}}=\frac{2m-1}{2}*\frac{\mathrm{Ha}}{60}$

${\mathrm{\&theta;}}_{\mathrm{mn}}=\frac{2n-1}{2}*\frac{\mathrm{\&pi;}}{60}---\left(15\right)$

故有网格(m，n)产生的脉冲激光回波为：

在网格(m，n)中，因为ρ²＜＜(H-h(m，n))²，故为了简化积分，可把积分号里面的

$\frac{\exp (-\frac{{(t-\frac{\mathrm{\&tau;}}{2}-\frac{2\sqrt{{\mathrm{\&rho;}}^2+{(H-h(m,n))}^2}}{c})}^2}{{2\mathrm{\&sigma;}}^2})}{{\mathrm{\&rho;}}^2+{(H-h(m,n))}^2}$

项中的(ρ，θ)直接用点(ρ_mn，θ_mn)来代替，这样的误差很小，所以结合(14)和(16)有

经过积分就有：

$\frac{\exp (-\frac{{(t-\frac{\mathrm{\&tau;}}{2}-\frac{2\sqrt{{(\frac{2m-1}{2}*\frac{\mathrm{Ha}}{60})}^2+{(H-h(m,n))}^2}}{c})}^2}{{2\mathrm{\&sigma;}}^2})}{{(\frac{2m-1}{2}*\frac{\mathrm{Ha}}{60})}^2+{(H-h(m,n))}^2}---\left(18\right)$

这样就可以把网格(m，n)的回波脉冲表示出来了，在利用(14)式就能够求出整个激光光斑内的脉冲回波。

所以，本发明所介绍的大脚印激光脉冲回波波形仿真算法主要思路如下：

结合星载激光高度计测高的具体条件，得到缓慢移动的收发合置测高系统对太空中没有大气衰减影响的漫反射大目标测高情况下的激光测高理论公式；

结合激光高度计采用的时间和空间上呈高斯分布的激光的特点，假设激光脉冲有无数多个小脉冲组成而得到激光脉冲回波波形积分公式；

结合大脚印地形的特点，对地形采用相应方法建模，得到易于编程实现的简单测高公式。

由以上介绍的仿真算法可知，本发明思路清晰，易于计算和编程实现，提出了目前在星载激光测高中大脚印激光脉冲回波波形的一种可实现仿真算法，对于实际大脚印激光测高的精度分析和反演地形有一定的借鉴意义。

附图说明：

图1为大脚印激光脉冲回波波形仿真流程图；

图2为大脚印地形分割示意图；

图3为大脚印为平地的回波波形；

图4为大脚印为斜坡(斜率为0.20)的回波波形；

图5为大脚印为四周平地中间有凸起的回波波形；

具体实施方式：

下面结合图1-5介绍本仿真算法的一个实施例。

参考中国科学院上海技术物理研究所制作的嫦娥一号探月卫星激光高度计的基本参数，选取高度计距离天底点的高度H＝200km，激光1/2束宽为0.3mrad，峰值功率为50mw，半峰值功率脉宽7ns，因为在月球表面不存在大气，故大气传输参数TA＝1。在不影响实验效果的情况下，其它参数一律取1。采用matlab软件编程实现该算法，本实施例的仿真流程如图1所示。

1、地形数据采样

本实施例对应大脚印为一个半径为60m的圆，其分割算法如图1所示。对于分割的地形网格(m，n)，选取点(ρ_mn，θ_mn)的高度代替整个网格，它并不一定恰好就存在于高程数据矩阵DEM所表示的地形离散点集中，所以还需要一个近似的过程，取该采样点(ρ_mn，θ_mn)距离高程数据矩阵DEM所表示地形中最近的那点高度DEM(i，J)作为网格(m，n)的高度h(m，n)，因为这两点的距离很近，200000米最大的距离也才几厘米，所以高度变化不可能太大，这种近似也是合理的。

本仿真实施例中限定高程数据矩阵大小DEM[1000][1000]，而地形分割后所得到的地形网格(m，n)的高度存储在采样矩阵h[60][120]中，由DEM[1000][1000]得到h[60][120]的方法简单叙述如下：

对于地形数据矩阵DEM[1000][1000]，取点(500，500)，(500，501)，(501，500)和(501，501)四点组成的小正方形中心点为天底点，这就是该高程图所描述的地形中心点。同时取该高程图像对应实际地形的大小为2Hα×2Hα，故图像的分辨率为2Hα/999，即相邻点最近距离Δh为2Hα/999。计算的脉冲回波对应的地形也不是高程数据矩阵对应的完整地形。而是这个正方形地形的内切圆，半径显然为Hα。

以天底点为圆心建立的直角坐标系，对于网格(m，n)中采样点(ρ_mn，θ_mn)，有

x_mn＝ρ_mn cosθ_mn，y_mn＝ρ_mn sinθ_mn

因为(x_mn，y_mn)一般不能直接对应直角坐标系中的某点整数，所以不能直接取到高程图中某一点高度作为采样矩阵h(m，n)的值。故先求出(x_mn，y_mn)与

DEM[1000][1000]所表示地形中最接近的离散点，就近似取那点的高程值。运用简单的坐标变换公式，可以找到这两者之间的对应关系如下：

h(m，n)＝DEM(500-i，501+j)；

其中

这样就完成了对地形数据的采样，得到了地形采样矩阵h[60][120]。

2、回波数据产生

根据公式(18)，可以得出任一网格(m，n)的脉冲回波功率，但是，在编程实现中，无法表示连续时间，所以只能用等间隔的离散点表示，而且，为了减少数据量和提高仿真算法运算效率，在激光脉冲发送到大脚印，大脚印上最高点最早反射脉冲回波回激光高度计之前这段时间内没有脉冲回波返回，所以返回功率为0，故可以不用计算。所以，回波数据产生算法简单叙述如下：

a.从h[60][120]中找出高度最大和最小的网格，分别记为h(i_max，J_max)、h(i_min，J_min)；

b.计算h(i_max，j_max)和h(i_min，j_min)对应激光小脉冲从发射返回到激光高度计的时间延迟，记为t_min，t_max；

c.确定计算脉冲回波的时间采样区间为[t_min-τ，t_max+τ]，τ为激光脉冲的半峰值功率脉宽，确定采样间隔为0.01ns，这时测高分辨率为1.5mm；

d.运用(18)式计算网格(m，n)在不同时间采样点下对应的回波数据，然后再根据(14)式把整个大脚印中的网格回波数据叠加起来，即得到整个脉冲回波的数据。

3、回波数据的处理和显示

根据上面得到的回波数据，为了后期激光测高精度分析等应用的方便，计算出理论上天底点在波形中对应的点和工程应用中采用半功率点测高方法得出的高度值，显示对应的回波波形并用绿色小三角形标出天底点在回波波形对应的位置和红色圆点标出回波的半功率点。具体结果如图3显示了一平地的回波波形，图4显示了一斜率为0.20的斜坡对应的回波波形，图5显示了一中间有凸起的地形的回波。

这样，就完成了大脚印激光脉冲回波波形的仿真，从结果我们可以看到，不同地形的回波波形有显著的不同，这对于实际测高中分析大脚印地形的类型和改善测高的精度有很好的借鉴作用。

Claims (1)

1.一种大脚印激光脉冲回波波形的仿真算法，其主要步骤为：

1)结合星载收发合置的激光高度计在基本无大气衰减的外太空中对表面漫反射的大目标一次测高过程中由于时间很短可近似为静止的特点，总结推导出其理论测高公式为：

$P_r=\frac{\mathrm{\&pi;f}{T_A}^2{D}^2{\mathrm{\&eta;}}_l{\mathrm{\&eta;}}_r}{4r^2}P,$

其中，Pr为接收的回波功率，Ps为发射的激光功率，f为地形后向散射率，TA为大气传输系数，D为高度计的接收孔径，ηt为发射机的光学效率，ηr为接收机的光学效率，r为目标到激光高度计的距离；

2)考虑到测高激光在时间和空间上的高斯特性，把发射的激光脉冲分解为无数多个小的入射到大脚印上各点的高斯脉冲，经大脚印对应点发射后被激光高度计接收，得到对应大脚印上激光脉冲回波功率的积分公式为：

其中，

$p_{s0}\left(t\right)=\frac{\sqrt{2}A}{{(1-\exp (-2))\mathrm{\&pi;}}^{3/2}\mathrm{\&sigma;}{H}^2{\mathrm{\&alpha;}}^2}\exp (-\frac{{(t-\frac{\mathrm{\&tau;}}{2})}^2}{2{\mathrm{\&sigma;}}^2}),$

$A=\frac{(\sqrt{\mathrm{\&pi;}}\mathrm{\&tau;ps}\_\max )}{2\sqrt{\ln 2}},$ $\mathrm{\&sigma;}=\frac{\mathrm{\&tau;}}{2\sqrt{2\ln 2}}$

与1)中公式相同的参数如无重复说明其物理意义与1)中一样，ps_max为激光峰值功率，τ为激光半峰值功率脉宽，α为激光1/2束宽，H为激光高度计到天底点的距离，r为大脚印中一点到激光高度计的距离；

3)结合大脚印地形的特点，进行一定的地形分割，得到激光脉冲回波功率的简单公式：

$\frac{\exp (-\frac{{(t-\frac{\mathrm{\&tau;}}{2}-\frac{2\sqrt{{(\frac{2m-1}{2}*\frac{\mathrm{Ha}}{60})}^2+{(H-h(m,n))}^2}}{c})}^2}{2{\mathrm{\&sigma;}}^2})}{{(\frac{2m-1}{2}*\frac{\mathrm{Ha}}{60})}^2+{(H-h(m,n))}^2}$

其中， $k=\frac{\mathrm{\&pi;f}{T_A}^2{D}^2{\mathrm{\&eta;}}_l{\mathrm{\&eta;}}_r}{4}*\frac{\sqrt{2}A}{(1-\exp (-2)){\mathrm{\&pi;}}^{3/2}\mathrm{\&sigma;}{H}^2{\mathrm{\&alpha;}}^2}$

(m，n)为地形分割中网格的序数，h(m，n)为网格(m，n)的高度值，经过计算即得到大脚印激光脉冲回波的波形。

Images