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

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

Info

Publication number
CN114509409B
CN114509409B CN202210412903.5A CN202210412903A CN114509409B CN 114509409 B CN114509409 B CN 114509409B CN 202210412903 A CN202210412903 A CN 202210412903A CN 114509409 B CN114509409 B CN 114509409B
Authority
CN
China
Prior art keywords
photon
water surface
laser radar
point cloud
satellite
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210412903.5A
Other languages
English (en)
Other versions
CN114509409A (zh
Inventor
马跃
郑慧莹
张智宇
宿殿鹏
阳凡林
黄珏
刘振
卜宪海
崔晓东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong University of Science and Technology
Original Assignee
Shandong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong University of Science and Technology filed Critical Shandong University of Science and Technology
Priority to CN202210412903.5A priority Critical patent/CN114509409B/zh
Publication of CN114509409A publication Critical patent/CN114509409A/zh
Application granted granted Critical
Publication of CN114509409B publication Critical patent/CN114509409B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • G01N21/49—Scattering, i.e. diffuse reflection within a body or fluid
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • G01N2021/4704—Angular selective
    • G01N2021/4709—Backscatter

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Chemical & Material Sciences (AREA)
  • Biochemistry (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Algebra (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Biology (AREA)
  • Remote Sensing (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Optical Radar Systems And Details Thereof (AREA)

Abstract

本发明公开了星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，属于激光遥感技术领域，用于获取海水参数、叶绿素参数的垂直剖面，包括：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；利用RL反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；构建光谱经验解析模型，计算水体叶绿素的垂直剖面。与现有技术相比，本发明能够准确反演待测水域的水体叶绿素垂直剖面结果，对于监测水体污染、藻类爆发、海洋碳循环、以及全球气候变化等多个领域具有重要意义。

Description

星载单光子激光雷达的叶绿素浓度垂直剖面反演方法
技术领域
本发明公开了星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，属于激光遥感技术领域。
背景技术
传统的星载线性体制激光雷达用于水体叶绿素探测时，主要处理海洋目标反射、散射的激光雷达回波波形数据；相对于传统线性体制激光雷达，单光子激光雷达探测器灵敏度提升了2-3个数量级，使得其垂直分辨率以及探测效率都有明显提升，通过卫星平台获取水下微弱的后向散射信号成为可能，进而能够获取海水参数、叶绿素参数的垂直剖面。但新一代单光子体制激光雷达不具备波形记录功能，只能记录单个光子信号的有无，不能记录信号的大小，所有被单光子探测器响应的入射光子不论是信号还是噪声，都对应一个空间三维光子点云。此外，单光子探测器所特有的余脉冲效应对辐射传输的非线性影响也必须被考虑和校正。
发明内容
本发明公开了星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，以解决现有技术中，单光子激光雷达不能记录信号大小、余脉冲效应对结果产生影响的问题。
星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，包括：
S1：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；
S2：利用RL反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；
S3：构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；
S4：构建光谱经验解析模型，计算水体叶绿素的垂直剖面。
优选地，所述S1包括：
S1.1.利用空间点云密度分类方法，提取密度更高的水面反射光子点云；
S1.2.通过水面光子点云，计算水面的波浪高度信息，用于确定水面反射信号光子在垂直方向的上下边界；
S1.3.以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；
S1.4.计算单次激光脉冲获取的平均水面光子点云数量以及在深度方向等间隔采样的水体后向散射信号光子点云序列。
优选地，所述S1.1包括：将单光子激光雷达获取的光子点云数据进行空间网格化，每个空间网格大小在沿卫星飞行方向的长度间隔为Δl，在高度方向的间隔为Δh，每个空间网格内的光子点云密度N(i，j)定义为沿卫星飞行方向第i个且沿高度方向第j个网格空间范围内的光子点云数量，计算所有空间网格点云密度的平均值μ和标准差σ，并以此为依据计算出水面反射信号光子点云的鉴别阈值TH，满足：
Figure 370465DEST_PATH_IMAGE001
（1）
Figure 615501DEST_PATH_IMAGE002
（2）
e a 是比例系数，取值为2~4，n l 为沿卫星飞行方向的网格数量n l =[max(l)-min(l)]/Δl，[ ]为向上取整运算符，l为光子点云沿卫星飞行方向的坐标；n h 为在高度方向的网格数量n h =[max(h)-min(h)]/Δh，h为光子点云沿高度方向的坐标，若空间网格密度N(i，j)大于鉴别阈值TH，则当前网格空间范围内的光子点云都被分类为水面反射的点云光子。
优选地，所述S1.2包括：移除靠近水面的激光雷达水体后向散射光子点云，选择水面有效波高SWH的下限作为水体光子点云的起算深度，具体为：
在沿卫星飞行方向每个长度为Δl的分段中，计算S1.1获取的水面反射信号光子的平均高度h pho 与标准差σ sub ，在每个Δl分段中，在垂直方向以平均高度h pho 为中心，采用水面有效波高SWH=4σ sub 作为波浪高度的垂直范围阈值，确定当前分段水面光子对应高度的上下边界[H min ，H max ]，筛选位于水面上下边界之间的光子点云，用于构建在沿卫星飞行方向每个长度为Δl分段内的水面反射信号光子点云集n s ，并构建水体后向散射信号光子集n u ，如式（3）所示，第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：
Figure 208288DEST_PATH_IMAGE003
Figure 665814DEST_PATH_IMAGE004
Figure 558814DEST_PATH_IMAGE005
Figure 709173DEST_PATH_IMAGE006
（3）
[l j ，h j ]是全部接收光子点云中第j个光子的坐标，l j 对应第j个光子沿卫星飞行方向的距离，h j 对应高度；L(i)是第i个沿卫星飞行方向分段的起始位置。
优选地，所述S1.3包括：
在沿卫星飞行方向累计连续k个分段的水体后向散射信号光子集，计算平均每次 激光脉冲获取光子点云在深度方向等间隔采样的水体后向散射信号点云数量序列N u (z)= [N u (z 1 )，N u (z 2 )，…，N u (z n )]，z代表深度序列，z 1 ，z 2 ，…，z n 代表不同深度，深度方向统计间隔
Figure 179469DEST_PATH_IMAGE007
，N u (z 1 )表示在深度z 1 处的平均每次激光脉冲后向散射对应的光子点云 数量，以相同的沿卫星飞行方向累计方式，计算平均每次激光脉冲在水面获取的光子点云 数量N s ，N s 对应平均水面位置，深度为0，即N s =N u (0)，合并表示，
Figure 456997DEST_PATH_IMAGE008
，此时z 1 =0。
优选地，所述S2包括：单光子激光雷达观测到的激光脉冲点云序列y(z)包括水面 反射平均信号点云数量N s 与水体后向散射平均信号点云数量序列N u (z)，是每次激光脉冲真 实获取的水面、水体光子点云序列
Figure 188193DEST_PATH_IMAGE009
与单光子激 光雷达系统的脉冲响应函数H(z)的卷积，表示为
Figure 260185DEST_PATH_IMAGE010
，*表示卷积运算，基于RL反 卷积算法通过迭代求取真实获取的水面、水体光子点云序列x(z)，它的第i次迭代计算的表 达为：
Figure 76832DEST_PATH_IMAGE011
（4）
x i+1 (z)为第i次迭代的真实光子点云序列的估计值，当残差
Figure 892472DEST_PATH_IMAGE012
小 于设定阈值或者迭代次数达到上限时迭代终止；H(z)为单光子激光雷达探测器的脉冲响应 函数，包括收发光学系统、单光子探测器件对发射和接收激光脉冲时间分布以及探测器件 温度和电压的相关影响，由激光雷达系统自身硬件决定，对于给定的单光子激光雷达是已 知或可算的；
当已知脉冲响应函数H(z)时，通过RL反卷积校正余脉冲效应，滤除主脉冲之后非 接收光子引起的小幅度余脉冲，经过反卷积后不同深度平均每次激光脉冲真实获取的光子 点云序列x(z)可以表示为
Figure 619120DEST_PATH_IMAGE013
，对应水体后向散射光子点云数量序列N u (z)经过反卷积后得到真实水体后向散射光子点云数量序列N u ′(z)、水面反射光子点云数 量N s 进行反卷积后得到真实水面反射光子点云数量N s ′。
优选地，所述S3包括：星载单光子激光雷达水面反射信号光子的概率统计模型为：
Figure 376860DEST_PATH_IMAGE014
（5）
η为单光子激光雷达系统的综合效率，E t 为单次脉冲激光发射能量，ρ s 为水面反射系数，对于532nm的绿激光为0.02，A r 为激光雷达望远镜接收孔径的有效面积，T a 为激光脉冲通过大气层时的单程衰减系数，h p 为普朗克常数约等于6.63×10-34J·s，v为激光频率，R为星载单光子激光雷达飞行高度，s 2 为水面均方斜率，其与海面上方风速U 10 有关，表示为：
Figure 697114DEST_PATH_IMAGE015
（6）
单光子激光雷达水体后向散射信号光子的概率统计模型为：
Figure 300134DEST_PATH_IMAGE016
（7）
exp为指数函数；β(π，z)为深度为z，角度为π时体散射系数；K d 代表水体漫射衰减系数；A为激光雷达系统常数，由单光子激光雷达系统参数和测量时对应的环境参数综合决定；n w 为水体折射率，对于532nm的绿激光为1.3334；T w 为激光脉冲经过水面进入水体时的衰减系数，对于532nm的绿激光为0.98，式（7）中，水深z的影响可以忽略；
结合式（5）和式（7），单光子激光雷达水面反射信号光子的概率统计模型改写为
Figure 146867DEST_PATH_IMAGE017
；n w 、ρ s 、T w 在给定条件下为常数值，水面均方斜率s 2 可以通过式（6）由 水面上方的风速U 10 计算得到，通过（7）式和水面反射激光脉冲平均接收光子数量N′ s 计算激 光雷达系统常数A。
优选地，所述S4包括：求解N u ′(z)随深度z变化的斜率得到K d ，将N u ′(z)，激光雷达系数常数A、漫射衰减系数K d 带入式（7）得到不同深度的β(π，z)。
优选地，所述S4包括：假设水体后向散射各向同性，不同深度z的水体后向散射系数b b 表示为b b (z)=2πβ(π，z)，散射系数b表示为b=b b /B，B为水体的后向散射比，表征了水体中后向散射在中总体散射的比例，对于清洁的大洋水体，后向散射比B为0.044，对于近海岸水体，后向散射比B为0.013。
优选地，所述S4包括：
海水散射系数b(λ)写成包含纯水的散射b w (λ)和海水中水体颗粒物的散射b p (λ)之 和的形式，随波长λ的变化，水体的散射系数b w (λ)关系表示为
Figure 560662DEST_PATH_IMAGE018
， 532nm时b w ≈0.0016，对于不同浓度的叶绿素海水，水体颗粒物的散射b p (λ)可以用来计算叶 绿素浓度C chl ：
Figure 351900DEST_PATH_IMAGE019
（8）
ϕ(λ)是与波长相关的常数，
Figure 509343DEST_PATH_IMAGE020
，根据在不同深度z处得 到的b p (λ，z)，代入式（8）计算得到在单光子激光雷达卫星当前地理位置处、不同深度z对应 的叶绿素剖面C chl (z)。
与现有技术相比，本发明在给定星载单光子激光雷达的系统硬件参数、测量光子点云数据与海面风速数据的前提下，利用水面反射的光子点云概率统计模型、水体（水下）散射光子点云概率统计模型、光谱解析经验模型，能够准确反演待测水域的水体叶绿素垂直剖面结果。
附图说明
图1是本发明的技术流程图；
图2是典型星载单光子激光雷达在开阔大洋区域、夜间测量获取的光子点云分布图；
图3是ICESat-2反演的叶绿素浓度结果与浮标实测数据在垂直方向的第一剖面结果对比图；
图4是ICESat-2反演的叶绿素浓度结果与浮标实测数据在垂直方向的第二剖面结果对比图。
具体实施方式
下面结合具体实施方式对本发明作进一步详细说明：
本发明的步骤如图1，以美国ICESat-2(Ice，Cloud，and land ElevationSatellite-2)单光子激光雷达所获取的光子点云数据及其系统硬件参数作为样例。ICESat-2是全球首颗搭载单光子体制激光雷达的卫星，其携带名为ATLAS(AdvancedTopographic Laser Altimeter System)的单光子激光雷达，其在垂直于卫星轨道方向分别配置了三束强波束和弱波束可以构建六条地形剖面，垂直精度高达几十厘米。同时，单光子激光雷达点云的足印大小、沿轨方向（沿卫星飞行方向）和垂轨方向（垂直于卫星飞行方向）的点云密度都大幅度提高，具有较高的水平分辨率和垂直分辨率，这极大的提升了卫星激光雷达的性能。ICESat-2/ATLAS提供的ATL03光子点云数据产品，包含所有测量光子的纬度、经度、WGS84椭球体高度、光子捕获时刻等信息。实施例的样例数据选择ICESat-2/ATLAS单光子激光雷达在2019年3月27日夜间飞越印度洋海域所获取的光子点云数据。实施例中的海面上方风速数据选择美国国家环境预报中心(National Centers for EnvironmentalPrediction，NCEP)全球范围再分析数据集（公开）提供的海面上方10m风速数据。
星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，包括：
S1：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；
S2：利用RL反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；
S3：构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；
S4：构建光谱经验解析模型，计算水体叶绿素的垂直剖面。
S1中水面和水体信号光子的分离提取主要包含三个步骤，首先利用空间点云密度分类方法，提取密度更高的水面反射光子点云（水面反射信号更强）；然后，通过水面反射光子点云，计算水面有效波高信息，用于确定水面反射信号光子在垂直方向的上下边界；进而，以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；最后，计算单次激光脉冲获取的平均水面光子点云数量、在深度方向等间隔采样的水体后向散射信号光子点云序列，所述S1包括：
S1.1.利用空间点云密度分类方法，提取密度更高的水面反射光子点云；
S1.2.通过水面光子点云，计算水面的波浪高度信息，用于确定水面反射信号光子在垂直方向的上下边界；
S1.3.以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；
S1.4.计算单次激光脉冲获取的平均水面光子点云数量以及在深度方向等间隔采样的水体后向散射信号光子点云序列。
所述S1.1包括：在夜间浅水区域，单光子激光雷达接收到信号光子主要由四部分组成：（1）大气后向散射信号：单光子激光雷达发射的激光脉冲经过大气后向散射进入激光雷达接收系统，被探测器响应；（2）水面反射信号：单光子激光雷达发射的激光脉冲经过大气衰减后到达水面，部分光子在水面反射后再次经过大气传输后进入接收系统并被探测器响应；（3）水体后向散射信号：部分信号进入水体内部，经水体中粒子的散射后再次经过水面和大气进入接收系统并被探测器响应（4）水底反射信号：当水深较浅时，部分发射的激光脉冲信号可以到达水底，经过水底目标反射，再次经过水体、水面以及大气后，进入接收系统并被探测器响应。为了有效分离提取出水体后向散射的光子信号点云，本发明首先对原始光子点云数据进行筛选，选择在远海开阔水域夜间测量得到的数据，排除水底反射激光脉冲的光子点云影响。将单光子激光雷达获取的光子点云数据进行空间网格化，每个空间网格大小在沿卫星飞行方向的长度间隔为Δl，在高度方向的间隔为Δh，每个空间网格内的光子点云密度N(i，j)定义为沿卫星飞行方向第i个且沿高度方向第j个网格空间范围内的光子点云数量，计算所有空间网格点云密度的平均值μ和标准差σ，并以此为依据计算出水面反射信号光子点云的鉴别阈值TH，满足：
Figure 804058DEST_PATH_IMAGE021
（1）
Figure 654334DEST_PATH_IMAGE022
（2）
e a 是比例系数，取值为2~4，n l 为沿卫星飞行方向的网格数量n l =[max(l)-min(l)]/Δl，[ ]为向上取整运算符，l为光子点云沿卫星飞行方向的坐标；n h 为在高度方向的网格数量n h =[max(h)-min(h)]/Δh，h为光子点云沿高度方向的坐标，若空间网格密度N(i，j)大于鉴别阈值TH，则当前网格空间范围内的光子点云都被分类为水面反射的点云光子。
所述S1.2包括：移除靠近水面的激光雷达水体后向散射光子点云，选择水面有效波高SWH的下限作为水体光子点云的起算深度，具体为：
由于水面可能存在强镜面反射，所以可能会出现探测器过饱和现象，而步骤4中的漫射衰减系数是基于水体（水下）光子的指数衰减特性求解的，所以靠近水面的激光雷达水体后向散射光子点云光子通常被移除，选择水面有效波高significant wave height(SWH)的下限作为水体（水下）光子点云的起算深度（如图2中黑色曲线所示），图2中黑色曲线为水体后向散射信号深度起算线，横轴是沿卫星飞行方向的纬度坐标，纵轴是高度方向，使用WGS84椭球高程基准，单位都为米。图中的网格是对光子点云进行的空间网格划分，空间网格大小在沿卫星飞行方向长度间隔为Δl，在高度方向间隔为Δh。
在沿卫星飞行方向每个长度为Δl的分段中，计算S1.1获取的水面反射信号光子的平均高度h pho 与标准差σ sub ，在每个Δl分段中，在垂直方向以平均高度h pho 为中心，采用水面有效波高SWH=4σ sub 作为波浪高度的垂直范围阈值，确定当前分段水面光子对应高度的上下边界[H min ，H max ]，筛选位于水面上下边界之间的光子点云，用于构建在沿卫星飞行方向每个长度为Δl分段内的水面反射信号光子点云集n s ，并构建水体后向散射信号光子集n u ，如式（3）所示，第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：
Figure 932869DEST_PATH_IMAGE023
Figure 628423DEST_PATH_IMAGE024
Figure 184169DEST_PATH_IMAGE005
Figure 189035DEST_PATH_IMAGE025
（3）
[l j ，h j ]是全部接收光子点云中第j个光子的坐标，l j 对应第j个光子沿卫星飞行方向的距离，h j 对应高度；L(i)是第i个沿卫星飞行方向分段的起始位置。
所述S1.3包括：
由于水体后向散射激光脉冲的信号较弱（即：光子点云数量较少），且单次或连续 几次激光脉冲测量具有一定的随机性，因此在沿卫星飞行方向累计连续k个分段的水体后 向散射信号光子集，计算平均每次激光脉冲获取光子点云在深度方向等间隔采样的水体后 向散射信号点云数量序列N u (z)=[N u (z 1 )，N u (z 2 )，…，N u (z n )]，z代表深度序列，z 1 ，z 2 ，…，z n 代表不同深度，深度方向统计间隔
Figure 705598DEST_PATH_IMAGE026
，N u (z 1 )表示在深度z 1 处的平均每次 激光脉冲后向散射对应的光子点云数量，以相同的沿卫星飞行方向累计方式，计算平均每 次激光脉冲在水面获取的光子点云数量N s ，N s 对应平均水面位置，深度为0，即N s =N u (0)，合 并表示
Figure 454111DEST_PATH_IMAGE027
，此时z 1 =0。
所述S2包括：步骤2中的反卷积过程是为了剔除耦合在水面反射、水体后向散射激 光脉冲形成光子点云之中的余脉冲影响。单光子激光雷达观测到的激光脉冲点云序列y(z) 包括水面反射平均信号点云数量N s 与水体后向散射平均信号点云数量序列N u (z)，是每次激 光脉冲真实获取的水面、水体光子点云序列
Figure 739730DEST_PATH_IMAGE009
与 单光子激光雷达系统的脉冲响应函数H(z)的卷积，表示为
Figure 181075DEST_PATH_IMAGE028
，*表示卷积运算， 基于RL反卷积算法通过迭代求取真实获取的水面、水体光子点云序列x(z)，它的第i次迭代 计算的表达为：
Figure 309568DEST_PATH_IMAGE029
（4）
x i+1 (z)为第i次迭代的真实光子点云序列的估计值，当残差
Figure 612505DEST_PATH_IMAGE030
小 于设定阈值或者迭代次数达到上限时迭代终止；H(z)为单光子激光雷达探测器的脉冲响应 函数，包括收发光学系统、单光子探测器件对发射和接收激光脉冲时间分布以及探测器件 温度和电压的相关影响，由激光雷达系统自身硬件决定，对于给定的单光子激光雷达是已 知或可算的；在本实验例中，有两种方法可以获取ICESat-2/ATLAS的脉冲响应函数：（1）通 过ICESat-2/ATLAS系统内部TEP(Transmitter Echo Path)模块所截取部分激光脉冲能量 对应的接收光子信号点云计算。（2）选用ICESat-2/ATLAS夜间的激光脉冲入射至平面目标 （沙漠区域）的接收光子点云信号计算。当已知脉冲响应函数H(z)时，通过RL反卷积校正余 脉冲效应，滤除主脉冲之后非接收光子引起的小幅度余脉冲，经过反卷积后不同深度平均 每次激光脉冲真实获取的光子点云序列x(z)可以表示为
Figure 1898DEST_PATH_IMAGE031
，对应水体后向散射光子点云数量序列N u (z)经过反卷积 后得到真实水体后向散射光子点云数量序列N u ′(z)、水面反射光子点云数量N s 进行反卷积 后得到真实水面反射光子点云数量N s ′。
所述S3包括：星载单光子激光雷达水面反射信号光子的概率统计模型为：
Figure 364877DEST_PATH_IMAGE032
（5）
η为单光子激光雷达系统的综合效率，E t 为单次脉冲激光发射能量，ρ s 为水面反射系数，对于532nm的绿激光为0.02，A r 为激光雷达望远镜接收孔径的有效面积，T a 为激光脉冲通过大气层时的单程衰减系数，h p 为普朗克常数约等于6.63×10-34J·s，v为激光频率，R为星载单光子激光雷达飞行高度，s 2 为水面均方斜率，其与海面上方风速U 10 有关，表示为：
Figure 839721DEST_PATH_IMAGE033
（6）
单光子激光雷达水体后向散射信号光子的概率统计模型为：
Figure 70982DEST_PATH_IMAGE034
（7）
exp为指数函数；β(π，z)为深度为z，角度为π时体散射系数；K d 代表水体漫射衰减系数；A为激光雷达系统常数，由单光子激光雷达系统参数和测量时对应的环境参数综合决定；n w 为水体折射率，对于532nm的绿激光为1.3334；T w 为激光脉冲经过水面进入水体时的衰减系数，对于532nm的绿激光为0.98，式（7）中，水深z的影响可以忽略；
结合式（5）和式（7），单光子激光雷达水面反射信号光子的概率统计模型改写为
Figure 65614DEST_PATH_IMAGE035
；n w 、ρ s 、T w 在给定条件下为常数值，水面均方斜率s 2 可以通过式（6）由水 面上方的风速U 10 计算得到，通过（7）式和水面反射激光脉冲平均接收光子数量N′ s 计算激光 雷达系统常数A。
所述S4包括：首先需要利用斜率法求解水体漫射衰减系数K d ，在公式（7）水体后向散射信号概率统计模型方程中有两个未知变量，深度为z，角度为π时（后向）的体散射系数β(π，z)和水体的漫射衰减系数K d ，水体信号呈指数衰减，假设漫射衰减系数K d 不随水深变化，可以通过求解N u ′(z)随深度z变化的斜率计算得到。求解N u ′(z)随深度z变化的斜率得到K d ，将N u ′(z)，激光雷达系数常数A、漫射衰减系数K d 带入式（7）得到不同深度的β(π，z)。
所述S4包括：假设水体后向散射各向同性，不同深度z的水体后向散射系数b b 表示为b b (z)=2πβ(π，z)，散射系数b表示为b=b b /B，B为水体的后向散射比，表征了水体中后向散射在中总体散射的比例，对于清洁的大洋水体，后向散射比B为0.044，对于近海岸水体，后向散射比B为0.013。
所述S4包括：
海水散射系数b(λ)写成包含纯水的散射b w (λ)和海水中水体颗粒物的散射b p (λ)之 和的形式，随波长λ的变化，水体的散射系数b w (λ)关系表示为
Figure 583183DEST_PATH_IMAGE036
， 532nm时b w ≈0.0016，对于不同浓度的叶绿素海水，水体颗粒物的散射b p (λ)可以用来计算叶 绿素浓度C chl ：
Figure 561635DEST_PATH_IMAGE037
（8）
ϕ(λ)是与波长相关的常数，
Figure 190062DEST_PATH_IMAGE038
，根据在不同深度z处 得到的b p (λ，z)，代入式（8）计算得到在单光子激光雷达卫星当前地理位置处、不同深度z对 应的叶绿素剖面C chl (z)。
根据前述步骤，在实施例中使用ICESat-2星载单光子激光雷达最终解算的印度洋待测水域叶绿素浓度垂直剖面，选择了BGC-ARGO编号为2902120的两条现场叶绿素浓度剖面测量数据，如图3和图4，两条剖面结果在不同深度对应的RMSE(Root Mean SquaredError均方根误差)分别为0.0583mg/m3、0.0399mg/m3，绝对值比例误差MAPE(Mean AbsolutePercentage Error)分别为8.35%、13.81%，这也证明了本发明方法的可行性和准确性。因此，本发明可以通过星载单光子激光雷达飞越不同水域所获取的光子点云数据，快速准确地计算当地叶绿素浓度的垂直剖面结果，这将有助于定量化研究海洋次表层叶绿素浓度的时空变化以及浮游植物生物量的研究。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.星载单光子激光雷达的叶绿素浓度垂直剖面反演方法，其特征在于，包括：
S1：沿卫星飞行方向，将星载单光子激光雷达获取的水面信号光子点云与水体后向散射信号光子点云进行分离提取；
S2：利用RL反卷积算法校正单光子激光雷达余脉冲效应对光子信号的影响；
S3：构建星载单光子激光雷达水面反射信号光子与水体后向散射信号光子的概率统计模型，利用风速数据计算激光雷达系统常数；
S4：构建光谱经验解析模型，计算水体叶绿素的垂直剖面；
所述S1包括：
S1.1.利用空间点云密度分类方法，提取密度更高的水面反射光子点云；
S1.2.通过水面光子点云，计算水面的波浪高度信息，用于确定水面反射信号光子在垂直方向的上下边界；
S1.3.以水面上下边界为高度阈值，在沿卫星飞行方向分段构建水面反射信号光子点云集和水体后向散射信号光子点云集；
S1.4.计算单次激光脉冲获取的平均水面光子点云数量以及在深度方向等间隔采样的水体后向散射信号光子点云序列；
所述S1.1包括：将单光子激光雷达获取的光子点云数据进行空间网格化，每个空间网格大小在沿卫星飞行方向的长度间隔为Δl，在高度方向的间隔为Δh，每个空间网格内的光子点云密度N(i，j)定义为沿卫星飞行方向第i个且沿高度方向第j个网格空间范围内的光子点云数量，计算所有空间网格点云密度的平均值μ和标准差σ，并以此为依据计算出水面反射信号光子点云的鉴别阈值TH，满足：
Figure DEST_PATH_IMAGE002
（1）
Figure DEST_PATH_IMAGE004
（2）
e a 是比例系数，取值为2~4，n l 为沿卫星飞行方向的网格数量n l =[max(l)-min(l)]/Δl，[ ]为向上取整运算符，l为光子点云沿卫星飞行方向的坐标；n h 为在高度方向的网格数量n h =[max(h)-min(h)]/Δh，h为光子点云沿高度方向的坐标，若空间网格密度N(i，j)大于鉴别阈值TH，则当前网格空间范围内的光子点云都被分类为水面反射的点云光子；
所述S1.2包括：移除靠近水面的激光雷达水体后向散射光子点云，选择水面有效波高SWH的下限作为水体光子点云的起算深度，具体为：
在沿卫星飞行方向每个长度为Δl的分段中，计算S1.1获取的水面反射信号光子的平均高度h pho 与标准差σ sub ，在每个Δl分段中，在垂直方向以平均高度h pho 为中心，采用水面有效波高SWH=4σ sub 作为波浪高度的垂直范围阈值，确定当前分段水面光子对应高度的上下边界[H min ，H max ]，筛选位于水面上下边界之间的光子点云，用于构建在沿卫星飞行方向每个长度为Δl分段内的水面反射信号光子点云集n s ，并构建水体后向散射信号光子集n u ，如式（3）所示，第i个分段内的水面反射信号光子集与水体后向散射信号光子集满足：
Figure DEST_PATH_IMAGE006
Figure DEST_PATH_IMAGE008
Figure DEST_PATH_IMAGE010
Figure DEST_PATH_IMAGE012
（3）
[l j ，h j ]是全部接收光子点云中第j个光子的坐标，l j 对应第j个光子沿卫星飞行方向的距离，h j 对应高度；L(i)是第i个沿卫星飞行方向分段的起始位置；
所述S1.3包括：
在沿卫星飞行方向累计连续k个分段的水体后向散射信号光子集，计算平均每次激光脉冲获取光子点云在深度方向等间隔采样的水体后向散射信号点云数量序列N u (z)=[N u (z 1 )，N u (z 2 )，…，N u (z n )]，z代表深度序列，z 1 ，z 2 ，…，z n 代表不同深度，深度方向统计间隔
Figure DEST_PATH_IMAGE014
，N u (z 1 )表示在深度z 1 处的平均每次激光脉冲后向散射对应的光子点云数量，以相同的沿卫星飞行方向累计方式，计算平均每次激光脉冲在水面获取的光子点云数量N s ，N s 对应平均水面位置，深度为0，即N s =N u (0)，合并表示
Figure DEST_PATH_IMAGE016
，此时z 1 =0；
所述S2包括：单光子激光雷达观测到的激光脉冲点云序列y(z)包括水面反射平均信号点云数量N s 与水体后向散射平均信号点云数量序列N u (z)，是每次激光脉冲真实获取的水面、水体光子点云序列
Figure DEST_PATH_IMAGE018
与单光子激光雷达系统的脉冲响应函数H(z)的卷积，表示为
Figure DEST_PATH_IMAGE020
，*表示卷积运算，基于RL反卷积算法通过迭代求取真实获取的水面、水体光子点云序列x(z)，它的第i次迭代计算的表达为：
Figure DEST_PATH_IMAGE022
（4）
x i+1 (z)为第i次迭代的真实光子点云序列的估计值，当残差
Figure DEST_PATH_IMAGE024
小于设定阈值或者迭代次数达到上限时迭代终止；H(z)为单光子激光雷达探测器的脉冲响应函数，包括收发光学系统、单光子探测器件对发射和接收激光脉冲时间分布以及探测器件温度和电压的相关影响，由激光雷达系统自身硬件决定，对于给定的单光子激光雷达是已知或可算的；
当已知脉冲响应函数H(z)时，通过RL反卷积校正余脉冲效应，滤除主脉冲之后非接收光子引起的小幅度余脉冲，经过反卷积后不同深度平均每次激光脉冲真实获取的光子点云序列x(z)可以表示为
Figure DEST_PATH_IMAGE026
，对应水体后向散射光子点云数量序列N u (z)经过反卷积后得到真实水体后向散射光子点云数量序列N u ′(z)、水面反射光子点云数量N s 进行反卷积后得到真实水面反射光子点云数量N s ′；
所述S3包括：星载单光子激光雷达水面反射信号光子的概率统计模型为：
Figure DEST_PATH_IMAGE028
（5）
η为单光子激光雷达系统的综合效率，E t 为单次脉冲激光发射能量，ρ s 为水面反射系数，对于532nm的绿激光为0.02，A r 为激光雷达望远镜接收孔径的有效面积，T a 为激光脉冲通过大气层时的单程衰减系数，h p 为普朗克常数约等于6.63×10-34J·s，v为激光频率，R为星载单光子激光雷达飞行高度，s 2 为水面均方斜率，其与海面上方风速U 10 有关，表示为：
Figure DEST_PATH_IMAGE030
（6）
单光子激光雷达水体后向散射信号光子的概率统计模型为：
Figure DEST_PATH_IMAGE032
（7）
exp为指数函数；β(π，z)为深度为z，角度为π时体散射系数；K d 代表水体漫射衰减系数；A为激光雷达系统常数，由单光子激光雷达系统参数和测量时对应的环境参数综合决定；n w 为水体折射率，对于532nm的绿激光为1.3334；T w 为激光脉冲经过水面进入水体时的衰减系数，对于532nm的绿激光为0.98，式（7）中，水深z的影响可以忽略；
结合式（5）和式（7），单光子激光雷达水面反射信号光子的概率统计模型改写为
Figure DEST_PATH_IMAGE034
；n w 、ρ s 、T w 在给定条件下为常数值，水面均方斜率s 2 可以通过式（6）由水面上方的风速U 10 计算得到，通过（7）式和水面反射激光脉冲平均接收光子数量N′ s 计算激光雷达系统常数A；
所述S4包括：求解N u ′(z)随深度z变化的斜率得到K d ，将N u ′(z)，激光雷达系数常数A、漫射衰减系数K d 带入式（7）得到不同深度的β(π，z)；
所述S4包括：假设水体后向散射各向同性，不同深度z的水体后向散射系数b b 表示为b b (z)=2πβ(π，z)，散射系数b表示为b=b b /B，B为水体的后向散射比，表征了水体中后向散射在总体散射的比例，对于清洁的大洋水体，后向散射比B为0.044，对于近海岸水体，后向散射比B为0.013；
所述S4包括：
海水散射系数b(λ)写成包含纯水的散射b w (λ)和海水中水体颗粒物的散射b p (λ)之和的形式，随波长λ的变化，水体的散射系数b w (λ)关系表示为
Figure DEST_PATH_IMAGE036
，532nm时b w ≈0.0016，对于不同浓度的叶绿素海水，水体颗粒物的散射b p (λ)可以用来计算叶绿素浓度C chl ：
Figure DEST_PATH_IMAGE038
（8）
ϕ(λ)是与波长相关的常数，
Figure DEST_PATH_IMAGE040
，根据在不同深度z处得到的b p (λ，z)，代入式（8）计算得到在单光子激光雷达卫星当前地理位置处、不同深度z对应的叶绿素剖面C chl (z)。
CN202210412903.5A 2022-04-20 2022-04-20 星载单光子激光雷达的叶绿素浓度垂直剖面反演方法 Active CN114509409B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210412903.5A CN114509409B (zh) 2022-04-20 2022-04-20 星载单光子激光雷达的叶绿素浓度垂直剖面反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210412903.5A CN114509409B (zh) 2022-04-20 2022-04-20 星载单光子激光雷达的叶绿素浓度垂直剖面反演方法

Publications (2)

Publication Number Publication Date
CN114509409A CN114509409A (zh) 2022-05-17
CN114509409B true CN114509409B (zh) 2022-08-02

Family

ID=81555217

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210412903.5A Active CN114509409B (zh) 2022-04-20 2022-04-20 星载单光子激光雷达的叶绿素浓度垂直剖面反演方法

Country Status (1)

Country Link
CN (1) CN114509409B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115266509B (zh) * 2022-09-26 2023-02-24 水利部交通运输部国家能源局南京水利科学研究院 一种基于激光雷达的水下垂向悬浮物浓度探测方法及系统
CN115855882B (zh) * 2023-03-03 2023-04-28 山东科技大学 利用星载激光雷达背景噪声反演水体遥感反射率的方法
CN117930203B (zh) * 2024-03-22 2024-06-04 山东科技大学 星载光子激光雷达的冰雪反射信号在轨辐射校正方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111965608A (zh) * 2020-07-16 2020-11-20 自然资源部第二海洋研究所 一种基于水体叶绿素浓度的星载海洋激光雷达探测能力评估方法
CN114355367A (zh) * 2022-01-10 2022-04-15 中国人民解放军61540部队 一种基于星载单光子激光雷达数据测量浅海水深的方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2001256611A1 (en) * 2000-05-16 2001-11-26 Jeacle Limited Photometric analysis of natural waters
CN109406457B (zh) * 2018-11-08 2021-04-06 北京航空航天大学 一种基于半分析模型的沉水植被光谱水体影响校正方法
CN110673108B (zh) * 2019-09-25 2021-07-16 自然资源部第二海洋研究所 一种基于迭代Klett的机载海洋激光雷达信号处理方法
CN111239713B (zh) * 2020-03-18 2022-03-04 武汉大学 一种星载单光子激光雷达的最大测量深度评估方法
US11815625B2 (en) * 2020-05-26 2023-11-14 China University Of Geosciences, Wuhan Methods and devices for correcting underwater photon displacement and for depth sounding with single-photon Lidar
CN114089366A (zh) * 2021-11-11 2022-02-25 山东科技大学 一种星载单光子激光雷达的水体光学参数反演方法
CN114167437B (zh) * 2021-11-22 2023-03-14 桂林理工大学 测水激光雷达多通道设计方法
CN114295585B (zh) * 2022-01-04 2024-03-22 浙江大学 一种基于解析模型的多视场海洋激光雷达数据正则化反演方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111965608A (zh) * 2020-07-16 2020-11-20 自然资源部第二海洋研究所 一种基于水体叶绿素浓度的星载海洋激光雷达探测能力评估方法
CN114355367A (zh) * 2022-01-10 2022-04-15 中国人民解放军61540部队 一种基于星载单光子激光雷达数据测量浅海水深的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Multivariate reconstruction of missing data in sea surface temperature, chlorophyll, and wind satellite fields;A. Alvera-Azcárate等;《JOURNAL OF GEOPHYSICAL RESEARCH》;20070315;1-11 *

Also Published As

Publication number Publication date
CN114509409A (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN114509409B (zh) 星载单光子激光雷达的叶绿素浓度垂直剖面反演方法
Smith et al. Land ice height-retrieval algorithm for NASA's ICESat-2 photon-counting laser altimeter
Abdallah et al. Wa-LiD: A new LiDAR simulator for waters
Ma et al. Estimating water levels and volumes of lakes dated back to the 1980s using Landsat imagery and photon-counting lidar datasets
Hsu et al. A semi-empirical scheme for bathymetric mapping in shallow water by ICESat-2 and Sentinel-2: A case study in the South China Sea
Zhang et al. An adaptive density-based model for extracting surface returns from photon-counting laser altimeter data
Tulldahl et al. Classification of aquatic macrovegetation and substrates with airborne lidar
Abdallah et al. Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters
Brunt et al. Determination of local slope on the Greenland Ice Sheet using a multibeam photon-counting Lidar in preparation for the ICESat-2 Mission
Zhang et al. A maximum bathymetric depth model to simulate satellite photon-counting lidar performance
Yueh et al. QuikSCAT geophysical model function for tropical cyclones and application to Hurricane Floyd
Moudrý et al. Effects of environmental conditions on ICESat-2 terrain and canopy heights retrievals in Central European mountains
CN114089366A (zh) 一种星载单光子激光雷达的水体光学参数反演方法
Steinvall et al. Experimental evaluation of an airborne depth-sounding lidar
CN111239713A (zh) 一种星载单光子激光雷达的最大测量深度评估方法
CN116817869B (zh) 一种利用激光雷达数据的海底光子信号确定方法
Farrell et al. Sea-ice freeboard retrieval using digital photon-counting laser altimetry
CN115508854A (zh) 一种利用星载单光子激光雷达海面点云反演水深的方法
Li et al. An improved method for automatic determination of the planetary boundary layer height based on lidar data
Zhang et al. An automatic algorithm to extract nearshore bathymetric photons using pre-pruning quadtree isolation for ICESat-2 data
CN114235173A (zh) 一种光子计数星载海洋激光雷达探测仿真方法
Wang et al. Determining planetary boundary layer height by micro-pulse lidar with validation by UAV measurements
Smith et al. Understanding biases in ICESat-2 data due to subsurface scattering using Airborne Topographic Mapper waveform data
Saputra et al. Effect of Turbidity, Temperature and Salinity of Waters on Depth Data from Airborne LiDAR Bathymetry
Streßer et al. Surface Wave and Roller Dissipation Observed With Shore‐Based Doppler Marine Radar

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant