星载单光子激光雷达的叶绿素浓度垂直剖面反演方法

1.本发明公开了星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，属于激光遥感技术领域。


背景技术：

2.传统的星载线性体制激光雷达用于水体叶绿素探测时，主要处理海洋目标反射、散射的激光雷达回波波形数据；相对于传统线性体制激光雷达，单光子激光雷达探测器灵敏度提升了2-3个数量级，使得其垂直分辨率以及探测效率都有明显提升，通过卫星平台获取水下微弱的后向散射信号成为可能，进而能够获取海水参数、叶绿素参数的垂直剖面。但新一代单光子体制激光雷达不具备波形记录功能，只能记录单个光子信号的有无，不能记录信号的大小，所有被单光子探测器响应的入射光子不论是信号还是噪声，都对应一个空间三维光子点云。此外，单光子探测器所特有的余脉冲效应对辐射传输的非线性影响也必须被考虑和校正。


技术实现要素：

3.本发明公开了星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，以解决现有技术中，单光子激光雷达不能记录信号大小、余脉冲效应对结果产生影响的问题。
4.星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，包括：s1：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；s2：利用rl反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；s3：构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；s4：构建光谱经验解析模型，计算水体叶绿素的垂直剖面。
5.优选地，所述s1包括：s1.1.利用空间点云密度分类方法，提取密度更高的水面反射光子点云；s1.2.通过水面光子点云，计算水面的波浪高度信息，用于确定水面反射信号光子在垂直方向的上下边界；s1.3.以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；s1.4.计算单次激光脉冲获取的平均水面光子点云数量以及在深度方向等间隔采样的水体后向散射信号光子点云序列。
6.优选地，所述s1.1包括：将单光子激光雷达获取的光子点云数据进行空间网格化，每个空间网格大小在沿卫星飞行方向的长度间隔为δl，在高度方向的间隔为δh，每个空间网格内的光子点云密度n(i，j)定义为沿卫星飞行方向第i个且沿高度方向第j个网格空间范围内的光子点云数量，计算所有空间网格点云密度的平均值μ和标准差σ，并以此为依
据计算出水面反射信号光子点云的鉴别阈值th，满足：（1）（2）ea是比例系数，取值为2~4，n
l
为沿卫星飞行方向的网格数量n
l
=[max(l)-min(l)]/δl，[ ]为向上取整运算符，l为光子点云沿卫星飞行方向的坐标；nh为在高度方向的网格数量nh=[max(h)-min(h)]/δh，h为光子点云沿高度方向的坐标，若空间网格密度n(i，j)大于鉴别阈值th，则当前网格空间范围内的光子点云都被分类为水面反射的点云光子。
[0007]
优选地，所述s1.2包括：移除靠近水面的激光雷达水体后向散射光子点云，选择水面有效波高swh的下限作为水体光子点云的起算深度，具体为：在沿卫星飞行方向每个长度为δl的分段中，计算s1.1获取的水面反射信号光子的平均高度h
pho
与标准差σ
sub
，在每个δl分段中，在垂直方向以平均高度h
pho
为中心，采用水面有效波高swh=4σ
sub
作为波浪高度的垂直范围阈值，确定当前分段水面光子对应高度的上下边界[h
min
，h
max
]，筛选位于水面上下边界之间的光子点云，用于构建在沿卫星飞行方向每个长度为δl分段内的水面反射信号光子点云集ns，并构建水体后向散射信号光子集nu，如式（3）所示，第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：（3）[lj，hj]是全部接收光子点云中第j个光子的坐标，lj对应第j个光子沿卫星飞行方向的距离，hj对应高度；l(i)是第i个沿卫星飞行方向分段的起始位置。
[0008]
优选地，所述s1.3包括：在沿卫星飞行方向累计连续k个分段的水体后向散射信号光子集，计算平均每次激光脉冲获取光子点云在深度方向等间隔采样的水体后向散射信号点云数量序列nu(z)=[nu(z1)，nu(z2)，
…
，nu(zn)]，z代表深度序列，z1，z2，
…
，zn代表不同深度，深度方向统计间隔，nu(z1)表示在深度z1处的平均每次激光脉冲后向散射对应的光子点云数量，以相同的沿卫星飞行方向累计方式，计算平均每次激光脉冲在水面获取的光子点云数量ns，ns对应平均水面位置，深度为0，即ns=nu(0)，合并表示，，此时z1=0。
[0009]
优选地，所述s2包括：单光子激光雷达观测到的激光脉冲点云序列y(z)包括水面反射平均信号点云数量ns与水体后向散射平均信号点云数量序列nu(z)，是每次激光脉冲真实获取的水面、水体光子点云序列与单光子激光雷达系统的脉冲响应函数h(z)的卷积，表示为，*表示卷积运算，基于rl反卷积算法通过迭代求取真实获取的水面、水体光子点云序列x(z)，它的第i次迭代计算的表
达为：（4）x
i+1
(z)为第i次迭代的真实光子点云序列的估计值，当残差小于设定阈值或者迭代次数达到上限时迭代终止；h(z)为单光子激光雷达探测器的脉冲响应函数，包括收发光学系统、单光子探测器件对发射和接收激光脉冲时间分布以及探测器件温度和电压的相关影响，由激光雷达系统自身硬件决定，对于给定的单光子激光雷达是已知或可算的；当已知脉冲响应函数h(z)时，通过rl反卷积校正余脉冲效应，滤除主脉冲之后非接收光子引起的小幅度余脉冲，经过反卷积后不同深度平均每次激光脉冲真实获取的光子点云序列x(z)可以表示为，对应水体后向散射光子点云数量序列nu(z)经过反卷积后得到真实水体后向散射光子点云数量序列nu′
(z)、水面反射光子点云数量ns进行反卷积后得到真实水面反射光子点云数量ns′
。
[0010]
优选地，所述s3包括：星载单光子激光雷达水面反射信号光子的概率统计模型为：（5）η为单光子激光雷达系统的综合效率，e
t
为单次脉冲激光发射能量，ρs为水面反射系数，对于532nm的绿激光为0.02，ar为激光雷达望远镜接收孔径的有效面积，ta为激光脉冲通过大气层时的单程衰减系数，h
p
为普朗克常数约等于6.63
×
10-34j·
s，v为激光频率，r为星载单光子激光雷达飞行高度，s2为水面均方斜率，其与海面上方风速u
10
有关，表示为：（6）单光子激光雷达水体后向散射信号光子的概率统计模型为：（7）exp为指数函数；β(π，z)为深度为z，角度为π时体散射系数；kd代表水体漫射衰减系数；a为激光雷达系统常数，由单光子激光雷达系统参数和测量时对应的环境参数综合决定；nw为水体折射率，对于532nm的绿激光为1.3334；tw为激光脉冲经过水面进入水体时的衰减系数，对于532nm的绿激光为0.98，式（7）中，水深z的影响可以忽略；结合式（5）和式（7），单光子激光雷达水面反射信号光子的概率统计模型改写为；nw、ρs、tw在给定条件下为常数值，水面均方斜率s2可以通过式（6）由水面上方的风速u
10
计算得到，通过（7）式和水面反射激光脉冲平均接收光子数量n
′s计算激光雷达系统常数a。
[0011]
优选地，所述s4包括：求解nu′
(z)随深度z变化的斜率得到kd，将nu′
(z)，激光雷达系数常数a、漫射衰减系数kd带入式（7）得到不同深度的β(π，z)。
[0012]
优选地，所述s4包括：假设水体后向散射各向同性，不同深度z的水体后向散射系数bb表示为bb(z)=2πβ(π，z)，散射系数b表示为b=bb/b，b为水体的后向散射比，表征了水体中后向散射在中总体散射的比例，对于清洁的大洋水体，后向散射比b为0.044，对于近海岸水体，后向散射比b为0.013。
[0013]
优选地，所述s4包括：海水散射系数b(λ)写成包含纯水的散射bw(λ)和海水中水体颗粒物的散射b
p
(λ)之和的形式，随波长λ的变化，水体的散射系数bw(λ)关系表示为，532nm时bw≈0.0016，对于不同浓度的叶绿素海水，水体颗粒物的散射b
p
(λ)可以用来计算叶绿素浓度c
chl
：（8）
ϕ
(λ)是与波长相关的常数，，根据在不同深度z处得到的b
p
(λ，z)，代入式（8）计算得到在单光子激光雷达卫星当前地理位置处、不同深度z对应的叶绿素剖面c
chl
(z)。
[0014]
与现有技术相比，本发明在给定星载单光子激光雷达的系统硬件参数、测量光子点云数据与海面风速数据的前提下，利用水面反射的光子点云概率统计模型、水体（水下）散射光子点云概率统计模型、光谱解析经验模型，能够准确反演待测水域的水体叶绿素垂直剖面结果。
附图说明
[0015]
图1是本发明的技术流程图；图2是典型星载单光子激光雷达在开阔大洋区域、夜间测量获取的光子点云分布图；图3是icesat-2反演的叶绿素浓度结果与浮标实测数据在垂直方向的第一剖面结果对比图；图4是icesat-2反演的叶绿素浓度结果与浮标实测数据在垂直方向的第二剖面结果对比图。
具体实施方式
[0016]
下面结合具体实施方式对本发明作进一步详细说明：本发明的步骤如图1，以美国icesat-2(ice，cloud，and land elevation satellite-2)单光子激光雷达所获取的光子点云数据及其系统硬件参数作为样例。icesat-2是全球首颗搭载单光子体制激光雷达的卫星，其携带名为atlas(advanced topographic laser altimeter system)的单光子激光雷达，其在垂直于卫星轨道方向分别配置了三束强波束和弱波束可以构建六条地形剖面，垂直精度高达几十厘米。同时，单光子激光雷达点云的足印大小、沿轨方向（沿卫星飞行方向）和垂轨方向（垂直于卫星飞行方向）的点云密度都大幅度提高，具有较高的水平分辨率和垂直分辨率，这极大的提升了卫星激光雷达的性能。icesat-2/atlas提供的atl03光子点云数据产品，包含所有测量光子的纬度、经度、wgs84椭球体高度、光子捕获时刻等信息。实施例的样例数据选择icesat-2/atlas
单光子激光雷达在2019年3月27日夜间飞越印度洋海域所获取的光子点云数据。实施例中的海面上方风速数据选择美国国家环境预报中心(national centers for environmental prediction，ncep)全球范围再分析数据集（公开）提供的海面上方10m风速数据。
[0017]
星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，包括：s1：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；s2：利用rl反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；s3：构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；s4：构建光谱经验解析模型，计算水体叶绿素的垂直剖面。
[0018]
s1中水面和水体信号光子的分离提取主要包含三个步骤，首先利用空间点云密度分类方法，提取密度更高的水面反射光子点云（水面反射信号更强）；然后，通过水面反射光子点云，计算水面有效波高信息，用于确定水面反射信号光子在垂直方向的上下边界；进而，以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；最后，计算单次激光脉冲获取的平均水面光子点云数量、在深度方向等间隔采样的水体后向散射信号光子点云序列，所述s1包括：s1.1.利用空间点云密度分类方法，提取密度更高的水面反射光子点云；s1.2.通过水面光子点云，计算水面的波浪高度信息，用于确定水面反射信号光子在垂直方向的上下边界；s1.3.以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；s1.4.计算单次激光脉冲获取的平均水面光子点云数量以及在深度方向等间隔采样的水体后向散射信号光子点云序列。
[0019]
所述s1.1包括：在夜间浅水区域，单光子激光雷达接收到信号光子主要由四部分组成：（1）大气后向散射信号：单光子激光雷达发射的激光脉冲经过大气后向散射进入激光雷达接收系统，被探测器响应；（2）水面反射信号：单光子激光雷达发射的激光脉冲经过大气衰减后到达水面，部分光子在水面反射后再次经过大气传输后进入接收系统并被探测器响应；（3）水体后向散射信号：部分信号进入水体内部，经水体中粒子的散射后再次经过水面和大气进入接收系统并被探测器响应（4）水底反射信号：当水深较浅时，部分发射的激光脉冲信号可以到达水底，经过水底目标反射，再次经过水体、水面以及大气后，进入接收系统并被探测器响应。为了有效分离提取出水体后向散射的光子信号点云，本发明首先对原始光子点云数据进行筛选，选择在远海开阔水域夜间测量得到的数据，排除水底反射激光脉冲的光子点云影响。将单光子激光雷达获取的光子点云数据进行空间网格化，每个空间网格大小在沿卫星飞行方向的长度间隔为δl，在高度方向的间隔为δh，每个空间网格内的光子点云密度n(i，j)定义为沿卫星飞行方向第i个且沿高度方向第j个网格空间范围内的光子点云数量，计算所有空间网格点云密度的平均值μ和标准差σ，并以此为依据计算出水面反射信号光子点云的鉴别阈值th，满足：
（1）（2）ea是比例系数，取值为2~4，n
l
为沿卫星飞行方向的网格数量n
l
=[max(l)-min(l)]/δl，[ ]为向上取整运算符，l为光子点云沿卫星飞行方向的坐标；nh为在高度方向的网格数量nh=[max(h)-min(h)]/δh，h为光子点云沿高度方向的坐标，若空间网格密度n(i，j)大于鉴别阈值th，则当前网格空间范围内的光子点云都被分类为水面反射的点云光子。
[0020]
所述s1.2包括：移除靠近水面的激光雷达水体后向散射光子点云，选择水面有效波高swh的下限作为水体光子点云的起算深度，具体为：由于水面可能存在强镜面反射，所以可能会出现探测器过饱和现象，而步骤4中的漫射衰减系数是基于水体（水下）光子的指数衰减特性求解的，所以靠近水面的激光雷达水体后向散射光子点云光子通常被移除，选择水面有效波高significant wave height(swh)的下限作为水体（水下）光子点云的起算深度（如图2中黑色曲线所示），图2中黑色曲线为水体后向散射信号深度起算线，横轴是沿卫星飞行方向的纬度坐标，纵轴是高度方向，使用wgs84椭球高程基准，单位都为米。图中的网格是对光子点云进行的空间网格划分，空间网格大小在沿卫星飞行方向长度间隔为δl，在高度方向间隔为δh。
[0021]
在沿卫星飞行方向每个长度为δl的分段中，计算s1.1获取的水面反射信号光子的平均高度h
pho
与标准差σ
sub
，在每个δl分段中，在垂直方向以平均高度h
pho
为中心，采用水面有效波高swh=4σ
sub
作为波浪高度的垂直范围阈值，确定当前分段水面光子对应高度的上下边界[h
min
，h
max
]，筛选位于水面上下边界之间的光子点云，用于构建在沿卫星飞行方向每个长度为δl分段内的水面反射信号光子点云集ns，并构建水体后向散射信号光子集nu，如式（3）所示，第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：（3）[lj，hj]是全部接收光子点云中第j个光子的坐标，lj对应第j个光子沿卫星飞行方向的距离，hj对应高度；l(i)是第i个沿卫星飞行方向分段的起始位置。
[0022]
所述s1.3包括：由于水体后向散射激光脉冲的信号较弱（即：光子点云数量较少），且单次或连续几次激光脉冲测量具有一定的随机性，因此在沿卫星飞行方向累计连续k个分段的水体后向散射信号光子集，计算平均每次激光脉冲获取光子点云在深度方向等间隔采样的水体后向散射信号点云数量序列nu(z)=[nu(z1)，nu(z2)，
…
，nu(zn)]，z代表深度序列，z1，z2，
…
，zn代表不同深度，深度方向统计间隔，nu(z1)表示在深度z1处的平均每次激光脉冲后向散射对应的光子点云数量，以相同的沿卫星飞行方向累计方式，计算平均每次激光脉冲在水面获取的光子点云数量ns，ns对应平均水面位置，深度为0，即ns=nu(0)，合
并表示，此时z1=0。
[0023]
所述s2包括：步骤2中的反卷积过程是为了剔除耦合在水面反射、水体后向散射激光脉冲形成光子点云之中的余脉冲影响。单光子激光雷达观测到的激光脉冲点云序列y(z)包括水面反射平均信号点云数量ns与水体后向散射平均信号点云数量序列nu(z)，是每次激光脉冲真实获取的水面、水体光子点云序列与单光子激光雷达系统的脉冲响应函数h(z)的卷积，表示为，*表示卷积运算，基于rl反卷积算法通过迭代求取真实获取的水面、水体光子点云序列x(z)，它的第i次迭代计算的表达为： （4）x
i+1
(z)为第i次迭代的真实光子点云序列的估计值，当残差小于设定阈值或者迭代次数达到上限时迭代终止；h(z)为单光子激光雷达探测器的脉冲响应函数，包括收发光学系统、单光子探测器件对发射和接收激光脉冲时间分布以及探测器件温度和电压的相关影响，由激光雷达系统自身硬件决定，对于给定的单光子激光雷达是已知或可算的；在本实验例中，有两种方法可以获取icesat-2/atlas的脉冲响应函数：（1）通过icesat-2/atlas系统内部tep(transmitter echo path)模块所截取部分激光脉冲能量对应的接收光子信号点云计算。（2）选用icesat-2/atlas夜间的激光脉冲入射至平面目标（沙漠区域）的接收光子点云信号计算。当已知脉冲响应函数h(z)时，通过rl反卷积校正余脉冲效应，滤除主脉冲之后非接收光子引起的小幅度余脉冲，经过反卷积后不同深度平均每次激光脉冲真实获取的光子点云序列x(z)可以表示为，对应水体后向散射光子点云数量序列nu(z)经过反卷积后得到真实水体后向散射光子点云数量序列nu′
(z)、水面反射光子点云数量ns进行反卷积后得到真实水面反射光子点云数量ns′
。
[0024]
所述s3包括：星载单光子激光雷达水面反射信号光子的概率统计模型为：（5）η为单光子激光雷达系统的综合效率，e
t
为单次脉冲激光发射能量，ρs为水面反射系数，对于532nm的绿激光为0.02，ar为激光雷达望远镜接收孔径的有效面积，ta为激光脉冲通过大气层时的单程衰减系数，h
p
为普朗克常数约等于6.63
×
10-34j·
s，v为激光频率，r为星载单光子激光雷达飞行高度，s2为水面均方斜率，其与海面上方风速u
10
有关，表示为：（6）单光子激光雷达水体后向散射信号光子的概率统计模型为：
（7）exp为指数函数；β(π，z)为深度为z，角度为π时体散射系数；kd代表水体漫射衰减系数；a为激光雷达系统常数，由单光子激光雷达系统参数和测量时对应的环境参数综合决定；nw为水体折射率，对于532nm的绿激光为1.3334；tw为激光脉冲经过水面进入水体时的衰减系数，对于532nm的绿激光为0.98，式（7）中，水深z的影响可以忽略；结合式（5）和式（7），单光子激光雷达水面反射信号光子的概率统计模型改写为；nw、ρs、tw在给定条件下为常数值，水面均方斜率s2可以通过式（6）由水面上方的风速u
10
计算得到，通过（7）式和水面反射激光脉冲平均接收光子数量n
′s计算激光雷达系统常数a。
[0025]
所述s4包括：首先需要利用斜率法求解水体漫射衰减系数kd，在公式（7）水体后向散射信号概率统计模型方程中有两个未知变量，深度为z，角度为π时（后向）的体散射系数β(π，z)和水体的漫射衰减系数kd，水体信号呈指数衰减，假设漫射衰减系数kd不随水深变化，可以通过求解nu′
(z)随深度z变化的斜率计算得到。求解nu′
(z)随深度z变化的斜率得到kd，将nu′
(z)，激光雷达系数常数a、漫射衰减系数kd带入式（7）得到不同深度的β(π，z)。
[0026]
所述s4包括：假设水体后向散射各向同性，不同深度z的水体后向散射系数bb表示为bb(z)=2πβ(π，z)，散射系数b表示为b=bb/b，b为水体的后向散射比，表征了水体中后向散射在中总体散射的比例，对于清洁的大洋水体，后向散射比b为0.044，对于近海岸水体，后向散射比b为0.013。
[0027]
所述s4包括：海水散射系数b(λ)写成包含纯水的散射bw(λ)和海水中水体颗粒物的散射b
p
(λ)之和的形式，随波长λ的变化，水体的散射系数bw(λ)关系表示为，532nm时bw≈0.0016，对于不同浓度的叶绿素海水，水体颗粒物的散射b
p
(λ)可以用来计算叶绿素浓度c
chl
：（8）
ϕ
(λ)是与波长相关的常数，，根据在不同深度z处得到的b
p
(λ，z)，代入式（8）计算得到在单光子激光雷达卫星当前地理位置处、不同深度z对应的叶绿素剖面c
chl
(z)。
[0028]
根据前述步骤，在实施例中使用icesat-2星载单光子激光雷达最终解算的印度洋待测水域叶绿素浓度垂直剖面，选择了bgc-argo编号为2902120的两条现场叶绿素浓度剖面测量数据，如图3和图4，两条剖面结果在不同深度对应的rmse(root mean squared error均方根误差)分别为0.0583mg/m3、0.0399mg/m3，绝对值比例误差mape(mean absolute percentage error)分别为8.35%、13.81%，这也证明了本发明方法的可行性和准确性。因此，本发明可以通过星载单光子激光雷达飞越不同水域所获取的光子点云数据，快速准确地计算当地叶绿素浓度的垂直剖面结果，这将有助于定量化研究海洋次表层叶绿素浓度的时空变化以及浮游植物生物量的研究。
[0029]
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。