本发明涉及激光测高技术领域，具体涉及一种卫星激光测高仪回波波形模拟仿真方法及设备。
背景技术：
卫星激光测高是一种新型的主动遥感技术，具有快速获取大范围地表三维信息的优势，冰、云和陆地高程卫星icesat(ice,cloudandlandelevationsatellite)上搭载的地球科学激光测高系统glas(geo-sciencelaseraltimetersystem)作为世界上第一个进行全球观测的全波形星载激光测高系统，在全球高程控制点获取、极地冰盖变化监测以及林业生物量估算等方面有广泛的应用，其中回波波形是主要的数据源。
卫星激光测高仪通过激光发射器向地面发射激光脉冲，经由大气传输和地表反射，再穿过大气返回到星上的激光接收系统，接收系统通过光电转换、数字化、波形采样等过程得到激光回波的波形数据。回波波形中既包含了地物高程信息，也包含了光斑内的地形和地物特征，是大光斑激光测高仪获取的核心数据，也是开展激光测高数据应用的基础。开展激光测高卫星在轨回波波形模拟仿真对于深入理解激光测高卫星的工作特性，指导开展卫星参数论证、分析测量精度、反演地形地物特征以及激光测高仪在轨几何检校等均具有重要的科学和实际应用价值。现有的接收系统采用全波形数字化记录(fullwaveform)，即对激光的发射信号和回波信号以很小的时间采样间隔(一般为1ns)进行采样并记录，但是由于激光在传递过程中，受地表反射率、发射波形、大气参数等的影响，因此现有技术中的卫星激光测高系统得到的回波波形通常会与真实的回波波形之间存在较大误差，影响后续的激光测高数据应用。
技术实现要素：
有鉴于此，本发明实施例提供了一种卫星激光测高仪回波波形模拟仿真方法及设备，以解决现有技术中的卫星激光测高系统得到的回波波形不准确的问题。
根据第一方面，本发明实施例提供了一种卫星激光测高仪回波波形模拟仿真方法，包括获取实际发射波形参数和激光光斑影像强度数据；利用所述激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息；根据所述激光地面足印中心的地理坐标及形状信息、机载激光雷达数据，提取激光光斑范围内的激光雷达点云数据坐标及强度信息；利用所述激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型；利用所述地表回波响应模型与所述实际发射波形参数得到卫星激光测高仪的回波波形数据。通过上述步骤，综合考虑了实际发射波形参数、激光光斑影像强度数据和地表反射率，基于上述数据来模拟仿真获取卫星激光测高仪的回波波形数据，与现有技术相比，该过程能够得到更加准确的回波波形数据，使其更加接近真实的回波波形。
结合第一方面，在第一方面第一实施方式中，所述利用所述激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息的步骤，包括：利用所述激光光斑影像强度数据，计算激光光斑的质心；利用激光脉冲能量分布特点和激光光斑的区域，计算描述激光光斑形状的椭圆形态参数；根据所述激光光斑的质心、所述椭圆形态参数及卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息。
结合第一方面第一实施方式，在第一方面第二实施方式中，在所述利用所述激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息的步骤之后，还包括：计算由于大气折射造成的激光地面足印中心的偏移量；根据所述偏移量，校正所述激光地面足印中心的地理坐标。由于大气对激光的折射影响会造成激光地面足印中心平面位置偏移，因此需要考虑该偏移量，以对激光地面足印中心的地理坐标进行校正，使得激光地面足印中心的地理坐标更加准确。
结合第一方面第一实施方式，在第一方面第三实施方式中，所述利用激光脉冲能量分布特点和激光光斑的区域，计算描述激光光斑形状的椭圆形态参数的步骤，包括：利用激光脉冲能量分布特点，根据预定阈值将所述激光光斑二值化；遍历二值化后的所述激光光斑的区域边界；根据所述区域边界的坐标，采用最小二乘法拟合椭圆方程，计算得到描述激光光斑形状的椭圆形态参数。
结合第一方面，在第一方面第四实施方式中，所述利用所述激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型的步骤，包括：利用所述激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率；利用激光雷达方程，计算所述不规则三角网模型的每个三角面元的回波光子数；按预设高程间隔将地形均匀分层，根据每一层的面元数量和所述每个三角面元的回波光子数，得到每层所述地形的能量响应值；将每层所述地形的能量响应值按时间序列排列，得到地表回波响应模型。该过程考虑到了地表反射率，从而使求得的地表回波响应模型更加精确。
结合第一方面第四实施方式，在第一方面第五实施方式中，采用以下公式计算所述不规则三角网模型的每个三角面元的回波光子数：
结合第一方面第四实施方式，在第一方面第六实施方式中，通过以下公式得到每层所述地形的能量响应值：
结合第一方面，在第一方面第七实施方式中，通过以下公式得到卫星激光测高仪的回波波形数据：e(t)＝t(t)*n(t)，其中，e(t)表示时间t时卫星激光测高仪的回波波形数据；t(t)表示时间t时所述实际发射波形参数；*表示卷积。
根据第二方面，一种计算机设备，包括至少一个处理器；以及与所述至少一个处理器通信连接的存储器；其中，所述存储器存储有可被所述至少一个处理器执行的计算机程序，所述计算机程序被所述至少一个处理器执行，以使所述至少一个处理器执行第一方面或第一方面任意一种实施方式所述的卫星激光测高仪回波波形模拟仿真方法。
根据第三方面，本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如第一方面或第一方面任意一种实施方式所述的卫星激光测高仪回波波形模拟仿真方法。
本发明实施例技术方案，具有如下优点：
本发明实施例提供了一种卫星激光测高仪回波波形模拟仿真方法，包括获取实际发射波形参数和激光光斑影像强度数据；利用激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息；根据该激光地面足印中心的地理坐标及形状信息、机载激光雷达数据，提取激光光斑范围内的激光雷达点云数据坐标及强度信息；利用激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型；利用地表回波响应模型与实际发射波形参数得到卫星激光测高仪的回波波形数据。通过上述步骤，综合考虑了发射波形参数、激光光斑影像强度数据以及地表反射率，基于上述数据来模拟仿真获取卫星激光测高仪的回波波形数据，与现有技术相比，该过程能够得到更加准确的回波波形数据，使其更加接近真实的回波波形。
附图说明
通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：
图1是根据本发明实施例的卫星激光测高仪回波波形模拟仿真方法的流程图；
图2是根据本发明实时例的icesat卫星工作期间glas激光脉冲能量空间分布图；
图3是根据本发明实施例的激光光斑和地面足印对应示意图；
图4是根据本发明实施例的激光光斑影像和地面足印形状示意图；
图5是根据本发明实施例的大气折射引起的激光平面位置偏移示意图；
图6是根据本发明实施例的卫星激光测高仪回波波形模拟仿真方法的具体流程图；
图7是本发明实施例提供的计算机设备的硬件结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例提供了一种卫星激光测高仪回波波形模拟仿真方法，图1是根据本发明实施例的卫星激光测高仪回波波形模拟仿真方法的流程图，如图1所示，该卫星激光测高仪回波波形模拟仿真方法包括：
步骤s101：获取实际发射波形参数和激光光斑影像强度数据；具体地，从卫星激光测高仪测量的数据中获取实际发射波形参数和激光光斑影像强度数据。
步骤s102：利用激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息；具体地，首先根据激光光斑影像数据描述激光光斑的形状，然后结合卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息。
步骤s103：根据激光地面足印中心的地理坐标及形状信息、机载激光雷达数据，提取激光光斑范围内的激光雷达点云数据坐标及强度信息；具体地，基于激光地面足印中心的地理坐标及形状信息，从机载激光雷达数据中提取激光足印中心所覆盖的点云数据，点云数据信息包括三维空间坐标以及点云数据的强度信息。
步骤s104：利用激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型；
步骤s105：利用地表回波响应模型与实际发射波形参数得到卫星激光测高仪的回波波形数据；将地表回波响应模型和实际发射波形参数综合考虑来得到卫星激光测高仪的回波波形数据，使结果更加准确。
通过上述步骤，获取实际发射波形参数和激光光斑影像强度数据，利用激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息；根据该激光地面足印中心的地理坐标及形状信息、机载激光雷达数据，提取激光光斑范围内的激光雷达点云数据坐标及强度信息；利用激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型；利用地表回波响应模型与实际发射波形参数得到卫星激光测高仪的回波波形数据。通过上述步骤，综合考虑了实际发射波形参数、激光光斑影像强度数据及地表反射率，基于上述数据来模拟仿真获取卫星激光测高仪的回波波形数据，与现有技术相比，该过程能够得到更加准确的回波波形数据，使其更加接近真实的回波波形。
上述步骤s101涉及到获取实际发射波形参数和激光光斑影像强度数据；具体地，从卫星激光测高仪测量的数据中获取实际发射波形参数和激光光斑影像强度数据。
从卫星激光测高仪的硬件测量相关数据中获取实际发射波形参数和激光光斑影像强度数据，如icesat/glas的glas01波形数据中获得发射波形t、从glas04工程数据中获得激光轮廓阵列lpa(laserprofilearray)对应的激光光斑影像强度数据，如图2所示，图2是根据本发明实时例的icesat卫星工作期间glas激光脉冲能量空间分布图。
上述步骤s102涉及到利用激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息，在一个具体实施方式中，该步骤包括：
利用激光光斑影像强度数据，计算激光光斑的质心。具体采用高斯拟合法计算lpa激光光斑影像强度数据的光斑质心，计算过程如下，由于激光能量分布近似符合高斯椭圆方程，用数学公式表达为：
公式(1)中，a为激光能量的幅值，(x0,y0)为光斑的质心，
对公式(1)取对数得：
将公式(2)展开为公式(3)的形式：
z＝ax2+by2+cx+dy+f(3)
公式(3)中：
解方程组(5)，求得系数a、b、c、d、f，
由方程组(5)的解，求得激光光斑的质心坐标为x0＝-2c/a，y0＝-2d/b，x、y方向上的标准差
得到激光光斑的质心坐标后，利用激光脉冲能量分布特点和激光光斑的区域，计算描述激光光斑形状的椭圆形态参数。具体地，利用激光脉冲能量分布特点，根据预定阈值将所述激光光斑二值化，并将激光光斑的有效区域和背景区域分离，遍历二值化后的激光光斑的有效区域的边界，根据区域边界的坐标，采用最小二乘法拟合椭圆方程，如公式(6)所示，计算得到描述激光光斑形状的椭圆形态参数：
ax2+bxy+cy2+dx+ey+f＝0(6)
公式(6)中，a、b、c、d、e、f为系数，(x,y)为激光光斑影像对应的像元坐标。
取f＝1，误差方程为：
v＝g-(ax2+bxy+cy2+dx+ey+1)(7)
将公式(7)写成矩阵形式：
c＝[x2xyy2xy]
x＝[abcde]t(8)
l＝-(ax2+bxy+cy2+dx+ey+1)
将公式(8)进行变形：
ctcx＝ctl(9)
x＝(ctc)-1ctl(10)
解得椭圆参数：
上述椭圆参数包括椭圆轴向θ、长轴a和短轴b，这些参数能描述激光光斑的形状。
根据激光光斑的质心、椭圆形态参数及卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息。具体地，卫星轨道数据包括地球参考坐标下卫星的扫描时刻、位置、速度数据，通过天球参考坐标系中的轨道六参数来描述卫星轨道：轨道高度、轨道倾角、轨道偏心率、近地点幅角、真近点角、升交点赤经，卫星姿态数据包括卫星平台的俯仰角、滚转角以及偏航角。卫星轨道姿态信息的模拟过程为现有方法即可实现，在此不再赘述。基于光线直线传播理论，求得激光地面足印中心地理坐标及形状信息，激光光斑和地面足印对应的示意图如图3所示，激光光斑影像和地面足印形状的示意图如图4所示。具体通过严密几何成像模型计算激光光斑在地球参考坐标系下的激光地面足印中心地理坐标：
公式(12)中，(x(t),y(t),z(t))t为成像时刻传感器投影中心在地球参考坐标系下的坐标，λ为积分时刻成像比例尺的分母(在本发明实施例中采用粗略的数字地形数据dem辅助获得，为现有技术)，
为了得到更准确的激光地面足印中心的地理坐标，在利用激光光斑影像强度数据和卫星轨道姿态信息，得到激光地面足印中心的地理坐标及形状信息的步骤之后，还包括计算由于大气折射造成的激光地面足印中心的偏移量，根据该偏移量，校正激光地面足印中心的地理坐标。具体地，激光在穿过大气层的过程中，由于受大气折射的影响，最终的地面足印中心的地理坐标位置与原始指向直线传播的位置会存在一定偏差，如图5所示。如果按激光指向角直线传播，最终地面足印点应该在p0点，而因大气折射导致激光地面足印中心的地理坐标的实际位置在p点。针对大气折射对光学卫星影像的几何位置偏差，采用视线跟踪几何算法进行计算，结果表明偏移量随着卫星观测角的增大而非线性迅速增大。偏移量的计算过程如下：
结合图5所示，根据折射定律，可得：
由于大气折射的原因，p相对于p0移动了距离pp0，
根据该偏移量pp0，校正激光地面足印中心的地理坐标。
得到激光地面足印中心的地理坐标后，上述步骤s104涉及到利用激光雷达点云数据坐标及强度信息分别得到不规则三角网模型和激光所在谱段的地表反射率，结合激光雷达方程，得到地表回波响应模型。在一个具体实施方式中，利用激光雷达点云数据坐标及强度信息，得到不规则三角网模型和地表反射率。具体地，基于激光地面足印中心的地理位置及形状，从1064nm波长的机载激光雷达lidar点云中提取激光足印所覆盖的点云数据信息，该激光雷达点云数据信息包括三维空间坐标以及强度信息pi(x,y,z,i),其中x,y,z为坐标信息，i为强度信息。利用激光雷达点云数据的强度信息经归一化后作为地表反射率，同时利用激光雷达点云数据中三维空间坐标来构建tin(triangleirregularnet)不规则三角网模型。激光雷达点云数据的强度信息经归一化后作为地表反射率的过程如下：
公式(15)中，imax为激光雷达点云数据的最大强度值，ii为激光雷达点云数据的强度值，ii'为归一化后的激光雷达点云数据的强度值，该值作为地表反射率。
然后利用激光雷达方程，计算不规则三角网模型的每个三角面元的回波光子数：
公式(16)中，ni为所述每个三角面元的回波光子数；ei为激光脉冲能量，从激光光斑影像中获取；h为普朗克常量；ν为激光频率；ar为望远镜接收孔径面积；ri为激光传播距离；τsys为光学系统透过率；τq为探测器量子效率；τatm为单程大气衰减系数；ρi为三角面元的地表反射率；θi为地表反射面法向量与望远镜视场方向的夹角。
设三角面元的三个顶点的三维坐标为a(x1,y1,z1)，b(x2,y2,z2)，c(x3,y3,z3)，则abc平面的法向量为:
公式(17)中，
该法向量与天顶方向的夹角为：
θi＝θs+θ(19)
其中θ为望远镜视场方向偏离天底方向的夹角，一般为0度。
在上述步骤之后，按预设高程间隔将地形均匀分层，根据每一层的面元数量和每个三角面元的回波光子数，得到每层地形的能量响应值，将每层地形的能量响应值按时间序列排列，得到地表回波响应模型。具体地，例如按一定高程间隔将地形分为m层，设每一层有k个面元，统计每一层地形的所有回波光子数得到该层的能量响应值，将所有地形层按时间序列排列起来就构成了地表回波响应模型n。
例如，激光足印范围内地形的高程范围为[hmin,hmax]，则m＝(hmax-hmin)/dh，一般取dh＝0.15m，对应时间间隔为1ns。第j层对应的高程为[hj,hj+dh]，对应的回波波形中的时间为[tj,tj+dh/0.15]。高程范围在[hj,hj+dh]内的所有面元对应的回波能量之和，即为该层的能量：
公式(20)中，n(tj)表示每层地形的能量响应值；j表示第j层地形；tj表示第j层地形对应的回波波形的时间；k表示每层地形的三角面元的总数。
地表回波相应模型为：
公式(21)中，n(t)表示时间t时地表回波响应模型；
然后将地表回波响应模型与步骤s101中提取的实际发射波形进行卷积，即可仿真获得卫星激光测高仪的回波波形数据：
e(t)＝t(t)*n(t)(22)
公式(22)中，e(t)表示时间t时卫星激光测高仪的回波波形数据；t(t)表示时间t时实际发射波形参数；*表示卷积。
本具体实施方式中还提供了一个回波波形仿真的具体流程图，如图6所示，从激光测高仪测量的辅助数据中剔除无效数据，然后从有效数据中获取激光发射波形、lpa强度数据和足印几何特征参数，根据原始机载lidar数据和激光足印位置数据，对其进行坐标变换，提取激光光斑范围内的lidar点云数据，剔除不完整的lidar点云数据，从完整的lidar数据中提取三角网内插dsm地形，从1064nm激光点云强度信息中获取地表发射率，然后根据激光发射波形、三角网内插dsm地形、地表反射率得到回波仿真模型，获取仿真回波波形，具体的详细计算机及获取过程同上述具体实施方式中的描述，在此不再赘述。综上所述，通过该过程获得更加精确的回波波形仿真数据，该回波波形数据是大光斑激光测高仪获取的核心数据，也是开展激光测高数据应用的基础，具有较好实际的应用价值。
本发明还实施例提供一种计算机设备，图7是本发明实施例提供的计算机设备的硬件结构示意图，如图7所示，该计算机设备包括一个或多个处理器71以及存储器72，图7中以一个处理器71为例。
该计算机设备还可以包括：输入装置73及输出装置74。
处理器71、存储器72、输入装置73及输出装置74可以通过总线或者其他方式连接，图7中以通过总线连接为例。
处理器71可以为中央处理器(centralprocessingunit，cpu)。处理器71还可以为其他通用处理器、数字信号处理器(digitalsignalprocessor，dsp)、专用集成电路(applicationspecificintegratedcircuit，asic)、现场可编程门阵列(field-programmablegatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等芯片，或者上述各类芯片的组合。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
存储器72作为一种非暂态计算机可读存储介质，可用于存储非暂态软件程序、非暂态计算机可执行程序以及模块，如本发明实施例中的卫星激光测高仪回波波形模拟仿真方法对应的程序指令/模块。处理器71通过运行存储在存储器72中的非暂态软件程序、指令以及模块，从而执行服务器的各种功能应用以及数据处理，即实现上述方法实施例所述的卫星激光测高仪回波波形模拟仿真方法。
存储器72可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储根据卫星激光测高仪回波波形模拟仿真装置的使用所创建的数据等。此外，存储器72可以包括高速随机存取存储器，还可以包括非暂态存储器，例如至少一个磁盘存储器件、闪存器件、或其他非暂态固态存储器件。在一些实施例中，存储器72可选包括相对于处理器71远程设置的存储器，这些远程存储器可以通过网络连接至卫星激光测高仪回波波形模拟仿真装置。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。
输入装置73可接收用户输入的查询请求(或其他数字或字符信息)，以及产生与卫星激光测高仪回波波形模拟仿真装置的用户设置以及功能控制有关的键信号输入。输出装置74可包括显示屏等显示设备，用以输出计算结果。
在本发明实施例中，可通过传感器等设备获取传动带的当前固有频率、张紧力等数据。
所述一个或者多个模块存储在所述存储器72中，当被所述一个或者多个处理器71执行时，执行如图1-图6所示的方法。
上述产品可执行本发明实施例所提供的方法，具备执行方法相应的功能模块和有益效果。未在本发明实施例中详尽描述的技术细节，具体可参见如图1-图6所示的实施例中的相关描述。
本发明实施例还提供一种非暂态计算机存储介质，所述计算机存储介质存储有计算机可执行指令，该计算机可执行指令可执行上述任意方法实施例中的卫星激光测高仪回波波形模拟仿真方法。其中，所述存储介质可为磁碟、光盘、只读存储记忆体(read-onlymemory，rom)、随机存储记忆体(randomaccessmemory，ram)、快闪存储器(flashmemory)、硬盘(harddiskdrive，缩写：hdd)或固态硬盘(solid-statedrive，ssd)等；所述存储介质还可以包括上述种类的存储器的组合。
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。