一种星载单光子激光雷达的水体光学参数反演方法

1.本发明属于激光遥感技术领域，尤其涉及一种星载单光子激光雷达的水体光学参数反演方法。


背景技术：

2.针对海洋次表层的研究，船载方式、浮标方式的原位测量可以提供高精度的结果，但是其耗费较大，探测效率低，很难大面积覆盖实施；星载被动光学遥感方式受制于太阳辐射和卫星轨道限制，较难在高纬度地区进行工作，并且无法提供有效的垂直剖面分布。主动探测型激光雷达遥感得益于主动发射激光脉冲，可以进行全天时工作，星载激光雷达又具有覆盖范围广等优点，逐渐成为了研究的可靠工具。传统线性体制激光雷达的极限探测能力与探测效率有限，而雪崩光电二极管(geiger-mode avalanche photodiode,gm-apd)和光电倍增管(photomultiplier tube,pmt)等单光子探测器能对单个光子的能量进行响应，比传统线性体制探测器的探测灵敏度提高了约1000倍，这极大的提升了激光雷达的性能。
3.在海洋次表层的光学参数反演层面，利用传统线性体制激光雷达获取海洋次表面光学参数已有较完善的理论、方法和实验，即：利用传统线性体制激光雷达所获取的全波形数据，通过建立波形参数与水体参数之间的物理或经验模型，提取和分析所测量波形的主要参数，反演水体衰减系数、散射系数。但线性体制激光雷达的水体光学参数反演方法并不适用于单光子体制，这是因为：单光子激光雷达并不具备波形记录功能，每次探测的结果为一个离散的光子事件，只能记录信号的有无，不能识别信号的大小；并且，由于单光子探测器灵敏度很高，受背景光影响较大，原始数据中含有大量的噪声光子，而在水体光学参数反演时仅需要水下后向散射信号光子；此外，单光子探测器受余脉冲效应的影响，会带来与传统线性体制探测器所不具有的额外噪声，需要在实际应用中进行修正。


技术实现要素：

4.本发明的目的在于提供一种星载单光子激光雷达的水体光学参数反演方法，以解决现有技术中，线性体制激光雷达的水体光学参数反演方法并不适用于单光子体制的问题。
5.一种星载单光子激光雷达的水体光学参数反演方法，包括：
6.s1.获取星载单光子激光雷达在水面上方飞行并进行测量时的有效信号；
7.s2.利用反卷积算法消除单光子激光雷达的余脉冲效应；
8.s3.基于水体后向散射的理论模型，通过最小二乘拟合获取水体后向体散射系数与水体漫射衰减系数；
9.s4.计算水体后向散射系数。
10.优选地，所述有效信号为水体后向散射点云，即水体后向散射激光脉冲形成的点云，通过空间点云密度分类方法，提取空间点云密度更高的水面反射点云，所述水面反射点云即水面反射激光脉冲形成的点云，比水面反射点云高度值更高的点云对应大气散射激光
脉冲形成的点云，比水面反射点云高度值更低的点云对应水体后向散射点云。
11.优选地，步骤s1包括以下子步骤：
12.s1.1.计算单位空间内的点云密度；
13.s1.2.计算空间网格内的平均值μ和标准差σ；
14.s1.3.计算水面反射点云与大气和水体后向散射点云的鉴别阈值th；
15.s1.4.计算平均每次激光脉冲获取点云在深度方向等间隔累计的后向散射激光脉冲点云序列y(z)＝[y(z1),y(z2),
…
,y(zn)]，z代表深度序列，z1,z2,
…
,zn代表不同深度，深度方向统计间隔δz＝z
2-z1＝z
3-z2＝z
4-z3。
[0016]
优选地，步骤s1.1包括：每个空间网格大小在沿卫星飞行方向长度间隔为δl，在高度方向间隔为δh，在研究区域星载单光子激光雷达所获取的数据沿卫星飞行方向的总长度为l，高度方向总长度为h，则在沿卫星飞行方向的网格数量n
l
和在高度方向的网格数量nh分别可以表示为：n
l
＝ceil(l/δl)，nh＝ceil(h/δh)，ceil表示向上取整函数；将研究区域内的数据进行空间网格划分，在沿卫星飞行方向的第i个且沿高度方向第j个网格的空间坐标范围f(i,j)满足：
[0017][0018]
l
start
和l
end
分别是研究区域内沿卫星飞行方向的起始坐标和终止坐标，h
min
和h
max
分别是研究区域内高度方向记录范围，即激光雷达的距离窗的最小值和最大值；统计每个空间网格坐标范围f(i,j)内的光子点云数量，即单位空间内的点云密度，记为n(i,j)。
[0019]
优选地，步骤s1.2中，
[0020]
优选地，步骤s1.3中，th＝μ+ea×
σ，ea是比例系数，在研究区域内，所有空间网格密度n(i,j)大于鉴别阈值th的空间范围内点云都分类为水面反射点云。
[0021]
优选地，步骤s2中，y(z)是真实水体后向散射信号序列x(z)＝[x(z1),x(z2),
…
,x(zn)]与单光子激光雷达系统的脉冲响应函数h(z)的卷积，即y(z)＝h(z)x(z)，z代表深度序列，y(z)＝h(z)x(z)等同于：
[0022][0023]
z1,z2,
…
,zn代表不同深度。
[0024]
优选地，单光子激光雷达系统的脉冲响应函数h(z)可以通过其激光脉冲入射至平面目标的接收信号计算获取，或者在单光子激光雷达系统内部截取部分激光脉冲能量对应的接收信号计算获取，当已知脉冲响应函数h(z)时，通过反卷积x(z)＝h-1
(z)y(z)获取消除余脉冲效应后的真实水下后向散射信号序列x(z)，其中h-1
(z)表示h(z)的逆矩阵。
[0025]
优选地，水体后向散射信号序列x(z)与水体参数后向散射理论函数可以表示为：
[0026]
[0027]
exp为指数函数，sec为正割三角函数；β(π)代表角度为π时(后向)的体散射系数，kd代表水体漫射衰减系数，为待计算的两项参数；e
t
为单脉冲激光发射能量，η为单光子激光雷达系统的综合效率，h
p
为普朗克常数约等于6.63
×
10-34j
·
s，v为激光波长效率，r为星载单光子激光雷达飞行高度，ar为激光雷达望远镜接收孔径有效面积，以上参数对于给定激光雷达系统都是已知参数；nw为水体折射率，θ为发射激光脉冲进入水体后的折射角，ta为激光脉冲通过大气层时的衰减系数，tw为激光脉冲通过水面进入水体时的衰减系数，以上参数通常在给定条件下为常数值或可以通过激光雷达系统参数计算获取；
[0028]
上式中，代求参数为2个，后向散射信号序列x(z)中的数据通常大于10个，存在数据冗余问题，为了取得更稳定的结果，抑制数据中存在的随机误差，采用最小二乘拟合的方法计算后向体散射系数β(π)，漫射衰减系数kd。
[0029]
优选地，步骤s4中的水体后向散射系数为更常用的水体光学参数，利用步骤3中得到的后向体散射系数β(π)，计算水体后向散射系数bb，假设后向散射各向同性条件下，后向散射系数bb定义为入射光线后半球空间的体散射系数之和：bb＝2πβ(π)。
[0030]
相对比现有技术，本发明建立一种通过理论建模反演水体光学参数的方法，只需给定单光子激光雷达的系统硬件参数与测量数据，结合信号统计模型、水体参数后向散射模型，能够准确反演待测水域次表层的水体光学参数，这对水体污染、水体富营养化、水质下降等进行有效监测具有重要的作用。
附图说明
[0031]
图1是本发明方法的技术流程图。
[0032]
图2是星载单光子激光雷达在开阔大洋区域夜间测量获取的光子点云分布图，横轴是沿卫星飞行轨道方向的相对累计距离，纵轴是高度方向，使用椭球高程基准，单位都为米，图中的网格是对光子点云进行的空间网格划分，空间网格大小在沿卫星飞行方向长度间隔为δl，在高度方向间隔为δh，右图是左图的局部放大展示。
[0033]
图3是实施例中利用icesat-2/atlas在撒哈拉沙漠区域夜间采集的点云数据计算获取的水深值zi与对应的脉冲响应函数h(zi)之间的曲线关系。
具体实施方式
[0034]
为使本发明的目的、技术方案和优点更加清楚，下面对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0035]
本发明以美国icesat-2(ice,cloud,and land elevation satellite-2)搭载全球首颗星载单光子激光雷达atlas(advance topographic laser altimeter system)的系统参数和所获取数据作为样例。icesat-2/atlas具有6个激光波束，激光脉冲重复频率10khz，沿轨激光采样间隔为0.7米，光斑直径为12米，多波束采样模式使得空间点云更为密集，对于粗糙和不均匀的地形将会提供更高的真实度。icesat-2/atlas提供的atl03数据产品，给出了其沿卫星飞行轨迹所获取的点云数据，点云数据中包含了icesat-2/atlas所获取的信号点和噪声点，每个点云具有唯一的三维坐标(经度、纬度和高度)，坐标系统基于
wgs84椭球。实施例的样例数据选择icesat-2/atlas单光子激光雷达在2018年12月13日晚上20:38:55(当地时间)飞越印度洋海域所获取的点云数据。
[0036]
一种星载单光子激光雷达的水体光学参数反演方法，包括：
[0037]
s1.获取星载单光子激光雷达在水面上方飞行并进行测量时的有效信号；
[0038]
s2.利用反卷积算法消除单光子激光雷达的余脉冲效应；
[0039]
s3.基于水体后向散射的理论模型，通过最小二乘拟合获取水体后向体散射系数与水体漫射衰减系数；
[0040]
s4.计算水体后向散射系数。
[0041]
优选地，所述有效信号为水体后向散射点云，即水体后向散射激光脉冲形成的点云，通过空间点云密度分类方法，提取空间点云密度更高的水面反射点云，所述水面反射点云即水面反射激光脉冲形成的点云，比水面反射点云高度值更高的点云对应大气散射激光脉冲形成的点云，比水面反射点云高度值更低的点云对应水体后向散射点云。
[0042]
优选地，步骤s1包括以下子步骤：
[0043]
s1.1.计算单位空间内的点云密度；
[0044]
s1.2.计算空间网格内的平均值μ和标准差σ；
[0045]
s1.3.计算水面反射点云与大气和水体后向散射点云的鉴别阈值th；
[0046]
s1.4.计算平均每次激光脉冲获取点云在深度方向等间隔累计的后向散射激光脉冲点云序列y(z)＝[y(z1),y(z2),
…
,y(zn)]，z代表深度序列，z1,z2,
…
,zn代表不同深度，深度方向统计间隔δz＝z
2-z1＝z
3-z2＝z
4-z3。
[0047]
单光子激光雷达在进行目标探测时，不能像传统线性体制探测器进行波形记录，但是可以记录每个光子点云的三维坐标，图2展示了星载单光子激光雷达在开阔大洋区域夜间测量获取的光子点云分布图；其中，横轴是沿卫星飞行轨道方向的相对累计距离，纵轴是高度方向，使用椭球高程基准，单位都为米。在进行水体目标探测时，星载单光子激光雷达接收到的光子信号主要由大气散射激光脉冲形成的点云、水面反射激光脉冲形成的点云、水体后向散射激光脉冲形成的点云、水底反射激光脉冲形成的点云、大气散射太阳背景辐射形成的点云、水面反射太阳背景辐射形成的点云、水体后向散射太阳背景辐射形成的点云、水底反射太阳背景辐射形成的点云、探测器暗计数形成的点云、探测器余脉冲形成的点云。通过单光子激光雷达获取水体光学参数，只有水体后向散射激光脉冲形成的点云属于信号，所有与太阳背景辐射相关的点云都属于噪声，大气散射、水面反射与水底反射激光脉冲形成的点云都属于噪声，探测器暗计数、激光雷达系统余脉冲形成点云也属于噪声。其中，探测器暗计数影响很小，可以忽略；激光雷达系统余脉冲影响耦合在水体后向散射激光脉冲形成点云之中，需要首先提取水体后向散射激光脉冲形成的点云。
[0048]
针对上述问题，为了更好的提取有效水体后向散射信号点云，本发明采用以下方法对数据进行处理和筛选：(1)选用在夜间测量得到的数据，排除所有与太阳背景辐射相关的影响；(2)选用远离岛礁、海岸带，在开阔海面测量得到的数据，排除水底反射激光脉冲的影响；(3)通过空间点云密度分类方法，提取空间点云密度更高的水面反射激光脉冲形成的点云，比这些水面点云高度值更高的点云对应大气散射激光脉冲形成的点云，比这些水面点云高度值更低的点云对应水体后向散射激光脉冲形成的点云。
[0049]
优选地，步骤s1.1包括：每个空间网格大小在沿卫星飞行方向长度间隔为δl，在
高度方向间隔为δh，在研究区域星载单光子激光雷达所获取的数据沿卫星飞行方向的总长度为l，高度方向总长度为h，则在沿卫星飞行方向的网格数量n
l
和在高度方向的网格数量nh分别可以表示为：n
l
＝ceil(l/δl)，nh＝ceil(h/δh)，ceil表示向上取整函数；将研究区域内的数据进行空间网格划分，在沿卫星飞行方向的第i个且沿高度方向第j个网格的空间坐标范围f(i,j)满足：
[0050][0051]
l
start
和l
end
分别是研究区域内沿卫星飞行方向的起始坐标和终止坐标，h
min
和h
max
分别是研究区域内高度方向记录范围，即激光雷达的距离窗的最小值和最大值；统计每个空间网格坐标范围f(i,j)内的光子点云数量，即单位空间内的点云密度，记为n(i,j)。
[0052]
优选地，步骤s1.2中，
[0053]
优选地，步骤s1.3中，由于水面反射系数远大于大气和水体的后向散射系数，因此空间点云密度较高网格内的点云对应水面反射点云，以此为依据计算出水面点云与大气和水体后向散射点云的鉴别阈值th＝μ+ea×
σ，ea是比例系数，在研究区域内，所有空间网格密度n(i, j)大于鉴别阈值th的空间范围内点云都分类为水面反射点云。
[0054]
优选地，步骤s1.4中，由于水体后向散射激光脉冲形成的点云信号较弱(点云数量较少)，且单次测量具有一定的随机性，因此在沿轨(卫星飞行)方向累计连续多个脉冲后，计算平均每次激光脉冲获取点云在深度方向等间隔累计的后向散射激光脉冲点云序列y(z)＝[y(z1),y(z2),
…
,y(zn)]，z代表深度序列，z1,z2,
…
,zn代表不同深度，深度方向统计间隔δz＝z
2-z1＝z
3-z2＝z
4-z3。在沿卫星飞行方向，不同位置的深度零点等于当前位置处分类为水面点云的空间网格坐标下限值。
[0055]
优选地，步骤s2中，y(z)是真实水体后向散射信号序列x(z)＝[x(z1),x(z2),
…
,x(zn)]与单光子激光雷达系统的脉冲响应函数h(z)的卷积，即y(z)＝h(z)x(z)，z代表深度序列，y(z)＝h(z)x(z)等同于：
[0056][0057]
z1,z2,
…
,zn代表不同深度。
[0058]
优选地，单光子激光雷达系统的脉冲响应函数h(z)可以通过其激光脉冲入射至平面目标的接收信号计算获取，或者在单光子激光雷达系统内部截取部分激光脉冲能量对应的接收信号计算获取。在本实施例中，有两种方法可以获取icesat-2/atlas的脉冲响应函数h(z)。首先，可以通过icesat-2/atlas系统内部tep(transmitter echo path)模块，其在内部截取部分激光脉冲能量，形成对应的接收信号，获取h(z)。此外，可以选用icesat-2/atlas在撒哈拉沙漠区域获得的夜间测量数据，这是因为：(1)夜间测量得到的数据排除所有与太阳背景辐射相关的影响；(2)icesat-2/atlas发射532nm波长的激光脉冲无法穿透沙漠表面；(3)撒哈拉沙漠表面覆盖物(例如植被)很少，且表面在12m激光光斑大小的尺寸上
较平坦。将icesat-2/atlas在撒哈拉沙漠区域夜间获得的点云数据沿卫星飞行轨迹进行累计，再沿深度方向进行分层，每层δz＝0.15m。其中，第一层z1对应点云数量最多的高度，其深度z1＝0，该层点云数量为h
′
(z1)；依次统计第i个分层的点云数量作为h
′
(zi)，对每一个分层点云数量进行归一化处理得到不同深度的h(zi)，即h(zi)＝h
′
(zi)/h
′
(z1)。这以百分比的形式定义了真实后向散射信号序列x(z)对单光子激光雷达所观测到的后向散射激光脉冲点云序列y(z)的贡献，图3显示了实施例中水深值zi与对应的脉冲响应函数h(zi)之间的曲线关系。当已知脉冲响应函数h(z)时，通过反卷积x(z)＝h-1
(z)y(z)获取消除余脉冲效应后的真实水下后向散射信号序列x(z)，其中h-1
(z)表示h(z)的逆矩阵。
[0059]
优选地，水体后向散射信号序列x(z)与水体参数后向散射理论函数可以表示为：
[0060][0061]
exp为指数函数，sec为正割三角函数；β(π)代表角度为π时(后向)的体散射系数，kd代表水体漫射衰减系数，为待计算的两项参数；e
t
为单脉冲激光发射能量，η为单光子激光雷达系统的综合效率，h
p
为普朗克常数约等于6.63
×
10-34j
·
s，v为激光波长效率，r为星载单光子激光雷达飞行高度，ar为激光雷达望远镜接收孔径有效面积，以上参数对于给定激光雷达系统都是已知参数；单脉冲发射能量e
t
＝93μj，单光子激光雷达系统的综合效率η＝0.0615，激光波长效率v＝3e8/532e-9
，激光雷达望远镜接收孔径有效面积ar＝0.41m2，星载单光子激光雷达的飞行高度r＝500km。nw为水体折射率，通常取1.334；θ为发射激光脉冲进入水体后的折射角，icesat-2/atlas为近天底入射，θ≈0
°
；ta为激光脉冲通过大气层时的衰减系数，icesat-2/atlas可以实时获取，在实施例中ta＝0.88；tw为激光脉冲通过水面进入水体的衰减系数，对于icesat-2/atlas的532nm激光波长，tw≈0.98，以上参数通常在给定条件下为常数值或可以通过激光雷达系统参数计算获取；
[0062]
上式中，代求参数为2个，后向散射信号序列x(z)中的数据通常大于10个，存在数据冗余问题，为了取得更稳定的结果，抑制数据中存在的随机误差，采用最小二乘拟合的方法计算后向体散射系数β(π)，漫射衰减系数kd。后向散射信号序列x(z)中数据通常大于10个(即n值)，在实施例中，水深范围10m，水深方向间隔δz＝0.15m，n＝67，存在数据冗余问题，为了取得更稳定的结果，抑制数据中存在的随机误差，采用最小二乘拟合的方法计算后向体散射系数β(π)，漫射衰减系数kd。进行最小二乘拟合时，需要给予后向体散射系数β(π)和漫射衰减系数kd的拟合初始值，可以选择典型水体光学参数作为初始值，532nm波长的光束对应不同类型水域的水体吸收系数a、散射系数b和后向散射系数bb表1所示。估算漫射衰减系数kd：kd＝a+4.18bb[1-0.52exp(-10.8a)]。
[0063]
优选地，步骤s4中的水体后向散射系数为更常用的水体光学参数，利用步骤3中得到的后向体散射系数β(π)，计算水体后向散射系数bb，假设后向散射各向同性条件下，后向散射系数bb定义为入射光线后半球空间的体散射系数之和：bb＝2πβ(π)。
[0064]
表1三类水体的典型光学参数数值(包括水体吸收系数、散射系数、后向散射系数参考值)
[0065] 清洁大洋水体近海岸水体港口浑浊水体吸收系数a(m-1
)0.1140.1790.366散射系数b(m-1
)0.0370.3981.824
后向散射系数bb(m-1
)0.001630.000520.0365
[0066]
根据前述步骤，在实施例中，最终解算印度洋待测水域的后向散射系数bb＝0.00247m-1
，漫射衰减系数kd＝0.1047m-1
，实际上modis在该水域所获取的后向散射系数bb＝0.00252m-1
，漫射衰减系数kd＝0.1036m-1
，与本发明所计算的结果非常一致，也证明了本发明方法的可行性和准确性。因此，本发明可以通过星载单光子激光雷达飞越海洋时所获取的点云数据，快速准确地计算当地水体的主要光学参数，作为进一步分析水体中各种组成成分和环境参数的基础。
[0067]
以上实施例仅用于说明本发明的技术方案，而非对其限制，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换，而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。