改进局部稀疏系数的轻小型光子计数激光雷达去噪方法

栾奎峰; 张昆宁; 邱振戈; 王洁; 王振华; 薛烨; 朱卫东; 林丹丹; 赵雪燕

doi:10.11834/jrs.20221854

模型与方法 | 浏览量 : 830 下载量: 678 CSCD: 0 更多指标

R-PDF
PDF
导出
分享给微信好友或者朋友圈

分享
收藏
专辑

改进局部稀疏系数的轻小型光子计数激光雷达去噪方法
Denoising method for light weight photon counting LiDAR based on an improved local sparse coefficient
- 角色： First author 第一作者
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
  上海河口海洋测绘工程技术研究中心, 上海 201306
- 邮箱： kfluan@shou.edu.cn
- 简介：栾奎峰，研究方向为新型激光成像质量控制理论和应用。E-mail： kfluan@shou.edu.cn
栾奎峰
1 2 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
张昆宁
1 ，
- 角色： Corresponding author 通讯作者
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
  上海河口海洋测绘工程技术研究中心, 上海 201306
- 邮箱： zgqiu@shou.edu.cn
- 简介：邱振戈，研究方向为海洋测绘与遥感。E-mail：zgqiu@shou.edu.cn
邱振戈
1 2 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
王洁
1 ，
- 机构：
  上海海洋大学信息学院, 上海 201306
王振华
3 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
薛烨
1 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
  上海河口海洋测绘工程技术研究中心, 上海 201306
朱卫东
1 2 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
林丹丹
1 ，
- 机构：
  上海海洋大学海洋科学学院, 上海 201306
赵雪燕
1
2023年27卷第2期页码：520-532
纸质出版日期： 2023-02-07 ，
DOI： 10.11834/jrs.20221854

扫描看全文

引用

阅读全文PDF
栾奎峰，张昆宁，邱振戈，王洁，王振华，薛烨，朱卫东，林丹丹，赵雪燕.2023.改进局部稀疏系数的轻小型光子计数激光雷达去噪方法.遥感学报，27（2）： 520-532

Luan K F，Zhang K N，Qiu Z G，Wang J，Wang Z H，Xue Y，Zhu W D，Ling D D and Zhao X Y. 2023. Denoising method for light weight photon counting LiDAR based on an improved local sparse coefficient. National Remote Sensing Bulletin， 27（2）：520-532
栾奎峰，张昆宁，邱振戈，王洁，王振华，薛烨，朱卫东，林丹丹，赵雪燕.2023.改进局部稀疏系数的轻小型光子计数激光雷达去噪方法.遥感学报，27（2）： 520-532 DOI： 10.11834/jrs.20221854.

Luan K F，Zhang K N，Qiu Z G，Wang J，Wang Z H，Xue Y，Zhu W D，Ling D D and Zhao X Y. 2023. Denoising method for light weight photon counting LiDAR based on an improved local sparse coefficient. National Remote Sensing Bulletin， 27（2）：520-532 DOI： 10.11834/jrs.20221854.

改进局部稀疏系数的轻小型光子计数激光雷达去噪方法

栾奎峰 ; 张昆宁 ; 邱振戈 ; 王洁 ; 王振华 ; 薛烨 ; 朱卫东 ; 林丹丹 ; 赵雪燕 ;

PDF下载

摘要

无人机载光子计数激光雷达测深系统具有高探测灵敏度、高密度、小光斑的特点，是海岛礁和浅海水深快速测量的重要技术手段。然而，高探测灵敏度也导致了获取的光子点云数据具有背景噪声大、信噪比与地物类型的相关性强以及光子的密度分布差异大，已有的去噪算法不能很好地适用。本文提出一种原始光子观测数据的去噪方法，首先基于直方图统计的方法计算原始光子观测数据的有效信号区间，其次利用网格统计法对区间内数据进行粗去噪，最后改进局部稀疏系数方法，采用水平椭圆搜索计算格网内每个光子数据的局部稀疏系数值，基于最大类间方差法确定噪声光子和信号光子的分隔阈值，实现原始光子观测数据的精去噪。本文选取了海南省的加井岛及临近的浅海地形为研究区，获取了无人机载光子计数激光雷达光子数据，验证本文提出的去噪算法。结果表明：该方法在高信噪比海岛植被覆盖区域和砂质潮间带区域的F1-Score均值达94.64%和98.96%，在低信噪比的近岸较浅和较深水体区域的F1-Score均值达93.04%和90.74%，总体F1-Score为94.34%，能有效剔除绝大部分噪声点，且对不同信噪比的海岛植被、沙地和不同深度的水下地形具有较强的适应能力。此外，本文还选取南海地区的珊瑚岛星载ICESat-2光子数据集，初步验证了本文所提出去噪算法在星载光子点云数据上的可用性。

译

Abstract

The photon counting LiDAR bathymetry system carried by UAVs is an important method for island reef mapping and shallow water bathymetry due to the characteristics of high detection sensitivity and high density. However, the high detection sensitivity also leads to the acquired photonic point cloud data with large background noise, a strong correlation between the signal-to-noise ratio and the type of ground objects, and large differences in the density distribution of photons, and the existing denoising algorithms cannot be well applied.

译

In this paper, a denoising method for raw photon observation data is proposed. First, the effective signal interval of the raw photon observation data is calculated based on the histogram statistics method, and then the data in the interval are coarsely denoised by the grid statistics method. Finally, the local sparse coefficient method is improved, the horizontal ellipse search is used to calculate the local sparse coefficient value of each photon data in the grid, and the method of maximum interclass variance is introduced to determine the separation threshold of noise photons and signal photons, which improves the original photon observation data. Denoising accuracy. Jiajing Island and the adjacent shallow sea terrain in Hainan Province are selected as the research area to verify the denoising algorithm proposed.

译

The results show that the average F1-score in the high signal-to-noise ratio areas, such as the island vegetation coverage area and the sandy intertidal zone, reaches 94.64% and 98.96%, respectively, and the average F1-score in the low signal-to-noise ratio area, such as the shallower and deeper water bodies near the coast, can also reach 93.04% and 90.74%, respectively. The overall F1-score is 94.34%, which can effectively remove most of the noise points and has strong adaptability to island vegetation, sandy land and underwater terrain of different depths with different signal-to-noise ratios.

译

In addition, this paper also selects the spaceborne ICESat-2 photon dataset of coral islands in the South China Sea, which further verifies the availability and applicability of the denoising algorithm proposed in this paper on spaceborne photonic point cloud data.

译

关键词

遥感; 光子计数; 无人机; 激光雷达; 局部稀疏系数; 最大类间方差法; 去噪

译

Keywords

remote sensing; photon counting; UAV; LiDAR; local sparse coefficient; OTSU; noise removal

译

1　引言

激光雷达是一种能直接获取目标高精度三维信息的主动式探测系统，在地形和水深测量（

刘帅等，2021；Wang等，2015）、建筑物的提取（Yang等，2014）和森林生物量的估计（李增元等，2016）等领域得到广泛的应用。传统线性探测体制的激光雷达由于对回波信号强度的要求较高，使得系统的体积、功耗增大，难以轻小型化，不适于近岸水下目标信息的快速获取（张河辉等，2019）。光子计数探测体制的激光雷达采用单光子探测器以及时间相关单光子计数TCSPC（Time-Correlated Single-Photon Counting）技术，能够识别和响应激光回波信号中一个或多个光子信号，实现对强度微弱且持续时间很短的信号的记录，减少激光发射脉冲所需要的能量，使得激光雷达系统的重量、体积和功耗大大降低（段雪洁，2017）。这些特点使得光子计数激光雷达可搭载在无人机载的飞行平台上，实现探测目标三维地理信息的快速获取。

译

由于光子计数激光雷达的高探测灵敏度，使得回波信号中包含了大量太阳辐射噪声、大气散射噪声、水体中后向散射噪声等背景辐射噪声，以及仪器所产生的暗计数噪声等，噪声来源情况复杂并与外部观测环境密切相关。同时由于水体对激光的强衰减作用，导致水下地形数据出现信号弱于噪声的低信噪比现象。虽然制造工艺和致冷方式的进步使仪器所产生的暗计数噪声得到有效控制，但如何有效抑制背景辐射噪声仍是光子计数激光雷达数据处理中的关键问题（

Fouche，2003；Shim和Lee，2014）。

译

国内外学者针对光子计数探测机制的激光雷达去噪方法开展了相关的研究。针对机载多波束测高实验激光雷达系统MABEL（Multiple Altimeter Beam Experimental Lidar）获取的光子点云数据，

夏少波等（2014）提出一种基于局部距离统计的方法，选用MABEL系统夜间时间段的低太阳光噪声点云数据进行去噪处理，整体的去噪精度达到97.8%。Zhang和Kerekes（2015）提出了一种基于密度聚类的噪声滤波方法，将密度聚类方法中的圆形搜索区域修改为椭圆，有效提升了在平坦地形区域的去噪精度。而相较于机载的MABEL系统，星载ICESat-2上搭载的先进地形激光高度计系统ATLAS（Advanced Topographic Laser Altimeter System）由于其飞行平台更高，虽然获取的数据覆盖范围更广，但返回的能量更弱，使得光子点云密度更低。Zhu等（2021）提出一种基于密度聚类的去噪方法OPTICS（Ordering Points To Identify the Clustering Structure），并利用水平椭圆搜索区域对其进行了改进，通过对获取的ATLAS光子点云数据去噪的结果进行定性和定量分析，验证该方法在植被区域提取信号光子的有效性。Chen等（2021）提出一种基于自适应椭圆滤波的光子计数激光雷达测深方法AVEBM（Adaptive Variable Ellipse Filtering Bathymetric Method），用于ATLAS系统获取的海岛近岸水体光子点云数据，能够精确检测水面和底部的信号光子。Xie等（2022）为了能够精确提取植被冠层的信号光子，提出了一种基于密度的自适应方法（DBAM）的地面和冠层光子检测方法，通过ICESat-2数据验证了方法的可行性以及提取信号光子的精确度。

译

现有的光子点云数据处理方法主要针对星载及有人机载的推扫式、高空低密度点云数据进行的滤波处理，但无人机载光子计数激光雷达系统飞行高度更低，返回的能量更强，足印更小，点密度差异也更大，再加上其采样的频率更高，使得接收到的信号和噪声的水平更高。

Wang等（2017）提出一种自适应椭球体搜索（AES）的去噪方法，应用于机载高分辨率量子激光雷达系统HRQLS（High Resolution Quantum Lidar System）获取的包含高背景噪声的三维光子点云数据中，提取出建筑物及植被区域有效信号光子点云。Shen等（2021）设计了一套高速机载单光子激光雷达系统，并提出一种基于网格密度统计的粗去噪方法，从高背景噪声的光子点云数据中提取出有效的地物信号光子。但上述两种方法对距离地物信号光子较近的噪声光子的剔除仍有进一步的改进空间。同时，无人机载轻小型光子计数激光雷达应用于包含植被、裸地、岩石、水面和水底等信噪比差异较大海岛及岸线区域观测，现有针对高程剖面光子点云数据噪声剔除的方法不能很好地适用于无人机载的轻小型光子计数激光雷达系统获取的原始观测数据。

译

本文提出一种改进局部稀疏系数的无人机载光子计数激光雷达原始观测数据去噪方法，实现复杂信噪比条件下噪声光子的精确去除，为光子计数激光雷达系统的海岛礁和浅海水深高精度测量提供可靠数据。

译

2　研究区和数据

2.1　研究区

为对不同地形、覆盖类型的海岛及浅海地形下的光子计数激光雷达数据产品质量进行控制，结合项目的可获取数据情况，本文选取海南省加井岛无人机载光子激光雷达数据开展实验。同时，为初步验证本文提出的光子点云数据去噪方法的适用性，选取了中国南海珊瑚岛区域的星载光子计数激光雷达ICESat-2数据开展验证实验，该岛位于中国西沙群岛永乐群岛中永乐环礁的西北部。加井岛的数据于2020年8月8日，在海南省万宁市加井岛及周边海域开展的无人机载轻小型光子计数激光雷达飞行试验获取。加井岛位于万宁市东部的石梅湾内，实验区域如图1所示，加井岛海岛面积约为0.18 km²，由两个低矮的山峰组成，岸线以陡峭山崖为主，西部为砂质岸线，水下是由多种类型珊瑚礁构成的复杂海底地貌。珊瑚岛是西沙群岛的一部分，主要由小岛、浅滩和珊瑚礁等组成。

译

图1 实验区域

Fig. 1 Experimental area

下载: 原图 | 高精图 | 低精图

2.2　实验数据

2.2.1　无人机载数据

无人机载轻小型光子计数激光雷达系统由中国科学院上海技术物理研究所研制，相关参数如表1所示，系统采用光子探测技术，激光波长为532 nm，具有一定的水体穿透能力。扫描模式为圆锥扫描，使每一束激光偏离天顶的角度相近，以保证每一束激光都有相近的水面入射角。飞行平台为四旋翼无人机，实验时平均飞行高度为120 m，平均飞行速度为4 m/s。实验获取了海岛及其附近海域内包括海岛植被、砂质潮间带和近岸2种水体区域在内的4种不同类型的光子数据。

译

表1 相关参数表

Table 1 Table of relevant parameters

参数	数值
单脉冲能量/nJ	300
波长/nm	532
重复频率/kHz	500
时间测量分辨率/ps	64
脉冲宽度/ps	300
反射镜转速/PRM	600

下载: 导出CSV

下载: 下载表格图片

系统获取的原始观测数据包含激光脉冲信息以及激光脉冲对应的时间计数信息。系统的激光重频为500 kHz，时间测量分辨率为64 ps，对应1 s的测量时间约有500000条激光脉冲，距离分辨率为0.96 cm。以激光雷达系统反射镜旋转一圈（0.1 s）获取的植被类型地物的原始观测数据为例，将每个激光回波脉冲接收到的一个或多个光子事件时间计数、脉冲序号信息分别解算成对应的距离（h）及时间信息（t），在二维平面上得到每个光子事件的时间与距离的对应关系图（图2）。

译

图2 植被覆盖区原始测距数据

Fig. 2 Raw ranging data of vegetation coverage area

下载: 原图 | 高精图 | 低精图

从图2中能直观的看到光子原始测距信号中存在大量的噪声，且水平方向上的信号光子点密度高于垂直方向上。在纵轴斜距方向上约0—30 m范围内，测距数据主要为激光回波在扫描反射镜上产生的窗口回波噪声，这部分噪声需要进行手动的剔除。在约80—100 m范围内，测距数据包括地物信号和噪声，其中密集的测距数据为信号光子部分，其余范围内为背景辐射、系统暗计数等所产生的噪声。

译

采用人工判读的方法对提取的有效信号范围内海岛及其附近水域的4种类型原始测距数据的信噪比SNR（Signal to Noise Ratio）进行了统计。每种类型选取10组时长0.1 s的数据，分别计算10组数据的平均值作为该地物类型光子数据的SNR值，其统计结果如表2所示。

译

表 2 不同地物类型的SNR

Table 2 SNR under different feature types

类型	平均SNR
海岛植被覆盖区域	11.62
砂质潮间带区域	6.19
近岸较浅水体区域	2.17
近岸较深水体区域	0.72

下载: 导出CSV

下载: 下载表格图片

从表2中可看出，在海岛上的植被覆盖区域与砂质潮间带区域的SNR值较高，分别为11.62和6.19。砂质潮间带区域（图3（b））主要以砂石等细小的颗粒物组成，对阳光有着更强的反射率，因此产生的大气背景辐射也就越多，使得激光器接收到太阳光背景辐射的噪声增多。海岛植被覆盖区域（图2（a））由于植被本身形成的不规则面多，加大了阳光反射的离散程度，使得对太阳光的反射率较低，因此产生的大气背景辐射相对较少，激光器接收到的太阳光背景辐射噪声也就减少（

Wu和Liu，2002），因此海岛植被覆盖区域的SNR高于砂质潮间带区域。

译

下载: | 高精图 | 低精图

图3 植被覆盖区与砂质潮间带区域

Fig. 3 Vegetation coverage area and sandy intertidal zone

在海域水体区域，近岸较浅水体区域与近岸较深水体区域的SNR均较低，分别为2.17和0.72，主要受水中各种微小悬浮物的影响以及水体本身对太阳光的反射、后向散射等现象的影响，使得水体中的噪声数量增多，导致水体的SNR值相对较低，并且随着深度增加（图4（b）），激光在海水中传输会发生严重的衰减（

Liu等，2020），导致水底反射的回波信号数量减少，使得近岸较深水体区域的SNR值低于近岸较浅水体区域（图4（a）），增加了该区域有效光子信号提取的难度。

译

下载: | 高精图 | 低精图

图4 近岸水体区域

Fig. 4 Nearshore water area

2.2.2　星载ICESat-2 数据

星载光子计数激光雷达ICESat-2由美国国家航天局NASA于2018年9月发射，其上搭载的先进地形激光测高仪ATLAS（Advanced Topographic Laser Altimeter System）采用多波束微脉冲光子计数的形式获取全球范围内地物表面的高程数据（

Markus 等，2017），相关参数见表3。

译

表3 ATLAS相关参数表

Table 3 Table of ATLAS relevant parameters

参数	数值
单脉冲能量/mJ	0.12—1.2
波长/nm	532
重复频率/kHz	10
脉冲宽度/ns	1
光斑直径/m	17
光斑间隔/m	0.7

下载: 导出CSV

下载: 下载表格图片

NASA根据数据的处理方式和处理程度的不同，将ICESat-2数据共计分成21级产品，分别命名为ATL01~21。选取位于研究区的ATL03数据集，ATL03数据为全球地理定位光子数据集，是ICESat-2数据的二级产品，数据中包含了激光器所接收到的每个光子的高程、纬度、精度和时间信息（

Neumann 等，2019）。本文在南海珊瑚岛研究区使用的ATL03数据集是获取于2019年2月22日的ATL03_2019022213519_08570207_005_01.‍h5，其中所用到的数据编号为GT3L。

译

通过经纬度将珊瑚岛研究区的范围截取出来，高程剖面点云结果如图5所示。横轴为沿轨距离、纵轴为光子的高程。图5中密集的光子为信号光子部分，且有明显的水面和水底形态特征。在水面和水底的信号光子周围存在着离散分布的噪声光子，其中在水面上的噪声光子明显少于水中的噪声光子。

译

图5 珊瑚岛区域剖面点云

Fig. 5 The profile point cloud of Shanhu Island

下载: 原图 | 高精图 | 低精图

3　研究方法

本文提出的轻小型光子计数激光雷达原始观测数据去噪方法，包括光子数据有效信号范围的确定、基于网格统计法的粗去噪和基于改进局部稀疏系数LSC（Local Sparsity Coefficient）方法的精去噪等几个过程，具体流程如图6所示。

译

图6 算法流程图

Fig. 6 Algorithm flow chart

下载: 原图 | 高精图 | 低精图

首先采用直方图统计的方法确定光子数据的有效信号范围，其次利用网格统计法消除低于平均点密度的格网内光子数据，去除远离信号中心的噪声，完成粗去噪。然后利用改进LSC的精去噪方法，计算剩余格网内的每个光子点的LSC值，并引入最大类间方差法（OTSU）（

Otsu，1979）确定LSC分割阈值，最后将大于分割阈值的噪声点去除，从而实现原始观测数据的精去噪。

译

3.1　有效信号范围确定

因原始观测数据中有效光子信号范围只占纵轴斜距方向上的小部分，为了减少计算量，根据纵轴斜距方向上信号分布的差异性，即有效信号范围内的光子数远高于其他仅为噪声光子的区域，本文采用直方图统计的方法确定光子信号的有效范围，并提取有效范围内的光子数据。

译

将实验区域内获得的单个文件内光子原始测距数据在纵轴方向上分割成一系列的空间垂直段，则分割单元的数量 $H_{m}$ 为：

$H_{m} = \frac{m a x (H_{i}) - m i n (H_{i})}{Δ H} (i = 1, 2, \dots, D)$

(1)

式中， $H_{i}$ 为数据集中每个点的距离值，D为数据集内点的个数， $Δ H$ 为垂直方向上分割单元的距离值， $Δ H$ 通常取经验值。首先统计每个垂向分割单元内点的数量 $p_{m}$ ，并计算所有垂向分割单元的平均光子数m，沿纵轴方向上生成距离—离散点数量统计直方图。然后利用高斯曲线拟合的方法确定直方图中的峰值位置 $H_{m a x}$ 和对应的距离值，并以距离值为中心分别向上向下选取最后一个 $p_{m}$ 大于m的节段，以此为边界得到有效光子信号的区间范围。

3.2　基于网格统计法的粗去噪

经上述方法处理后得到的光子数据包含了信号光子和噪声光子，本文采用基于网格统计的粗去噪方法，分别在横轴及纵轴方向上进行等时间（ $Δ t$ ）、等距离（ $Δ h$ ）间隔的划分，网格数量M为

$\{\begin{array}{l} \begin{array}{l} t_{m} = \frac{m a x (t_{i}) - m i n (t_{i})}{Δ t} \\ h_{m} = \frac{m a x (h_{i}) - m i n (h_{i})}{Δ h} \end{array} \\ M = t_{m} \times h_{m} \end{array}$

(2)

式中， $t_{i}$ 为光子的时间信息， $h_{i}$ 为光子的距离信息， $t_{m}$ 与 $h_{m}$ 分别为横纵轴上划分间隔的数量， $Δ t$ 与 $Δ h$ 采取试错法进行确定。方法首先统计每个网格单元中光子数据的数量n，并计算网格单元的平均光子数 $h_{a v e_m}$ 。噪声光子随机分布在整个网格空间中，并且包含信号光子和噪声光子网格内的光子数量明高于只包含噪声光子的网格（

Popescu等，2018；Ma等，2018）。基于此思想，本文将n小于

$h_{a v e_m}$ 的网格作为噪声网格并剔除网格内的所有光子数据，大于

$h_{a v e_m}$ 的则作为有效信号网格保留，实现远离信号光子中心噪声点的粗剔除。

3.3　改进LSC方法精去噪

粗去噪处理后保留的网格内除包含信号光子的同时，邻近信号的噪声光子也被保留了下来，但信号光子和噪声光子数据仍具有显著的密度差异性，可以使用基于密度的LSC算法对网格内的噪声光子进行处理。LSC算法基本思想是通过计算待检测点在其圆形邻域内点的距离与数量的比值关系得到LSC值，LSC值越大，则该点为噪声光子的可能性就越大，反之该点为信号光子的可能性就越大（

Agyemang，2002）。

译

经粗去噪后保留的网格内信号光子在水平方向的聚合度高于垂直方向，噪声光子离散随机分布于信号光子周围，原LSC方法中综合度量了待检测光子在圆形搜索邻域内光子的数量和局部距离值，但并未考虑信号光子与噪声光子在空间分布上的形态差异，无法准确判别垂直测距方向上距离信号光子较近的噪声光子。因此本文提出了改进水平椭圆搜索的LSC精去噪方法，从而提升了噪声光子在垂直测距方向上的剔除精度。在2.2.1节建立的坐标系中，对于给定任意待检测光子 $x (t_{x}, h_{x})$ 和其他任意光子 $y (t_{y}, h_{y})$ ，椭圆搜索区域由以下等式定义：

$d i s t (x, y) = \sqrt[]{\frac{(t_{x} - t_{y})^{2}}{a^{2}} + \frac{(h_{x} - h_{y})^{2}}{b^{2}}}$

(3)

式中，a和b分别表示椭圆长轴和短轴，椭圆搜索形状（扁率）由长轴和短轴的比值决定。

译

椭圆搜索的长短半轴比值受实际获取数据的地形特征和地面反射率影响较大，且海岛礁和浅海水域的噪声光子数据空间分布形态复杂，无法采用自适应的方法来确定长短半轴的比值。一般情况下，当信号光子在水平方向的聚合度远大于噪声光子时，信噪比较高，此时椭圆长短半轴的比值可取较小值，以便较好的去除噪声光子。当信号光子和噪声光子在水平方向的聚合度差异较小时，信噪比较低，此时椭圆长短半轴的比值可取较大值，以减少信号光子的误剔除。本文根据实际获取的数据集中信号光子和噪声光子的聚合程度和空间分布形态情况，取a∶b=10∶1。LSC方法中的初始阈值k决定着搜索区域大小，并影响着最终LSC值的计算结果。在SNR较高的植被覆盖区域与砂质潮间带区域，当初始阈值k偏低时，有利于去除离信号光子较近的噪声光子，在SNR较低的近岸水体区，当初始阈值k偏高时，有利于水下信号光子的提取。

译

利用构建的椭圆搜索区域，计算粗去噪后保留下的网格内每个光子的局部稀疏率Lsr（Local sparsity rate）和裁剪因子pf（Pruning factor of p），将Lsr小于pf的光子剔除，余下的光子放入候选集中，最终计算得到候选集中每个光子的LSC值（

Agyemang，2004）。用LSC值表征待检测光子与周围光子的聚合程度，该值越小表明待检测光子与周围光子的聚合程度高，成为信号光子的可能性则越大，原LSC方法中需要通过经验法取值来确定分割信号光子和噪声光子的LSC值。但在包含海岛礁和浅海水域的复杂信噪比区域经验法取值存在不准确性和适用性较差的问题，本文采用基于最大类间方差法确定噪声光子和信号光子的分割阈值，采用理论计算得到最大类间方差值，并使用其对应的LSC值为分割阈值，从而提升分割阈值的准确性，基本公式如下：

译

$\{\begin{array}{l} ω_{1} (m) = \frac{m}{N} \\ ω_{2} (m) = \frac{N - m}{N} \\ μ_{1} (m) = \frac{\sum_{i = 1}^{N} L S C (i)}{N - m} \\ μ_{2} (m) = \frac{\sum_{i = m + 1}^{N} L S C (i)}{N - m} \\ μ (m) = ω_{1} (m) \cdot μ_{1} (m) + ω_{2} (m) \cdot μ_{2} (m) \\ \begin{array}{l} g = ω_{1} (m) \cdot (μ (m) - μ_{1} {(m))}^{2} + \\ ω_{2} (m) \cdot (μ (m) - μ_{2} {(m))}^{2} \end{array} \end{array}$

(4)

式中，N是候选集中待检测光子的数量，m为候选集中所有的待检测光子按其LSC值从小到大排序后的光子序号，范围从1到N。假设前m个光子构成“信号”类，其余光子构成“噪声”类，属于信号类的光子数量占整个候选集光子数的比例记为 $ω_{1} (m)$ ，其平均的LSC值为 $μ_{1} (m)$ 。属于噪声类的光子数量占整个候选集光子数的比例记为 $ω_{2} (m)$ ，其平均的LSC值为 $μ_{2} (m)$ ，此时整个候选集的平均LSC值为 $μ (m)$ ，类间方差值记为g。循环计算得到N个类间方差值g，其中g最大时所对应的LSC值即为候选集中噪声光子与信号光子的分割阈值。最后将候选集中小于LSC分割阈值的光子标记为信号光子，反之标记为噪声光子。

4　结果与分析

4.1　无人机载光子数据去噪结果与分析

4.1.1　粗去噪结果与分析

实验首先采用直方图统计的方法确定了4种类型实验数据光子信号的有效范围，其中 $Δ H$ 的取值为0.5 m，部分提取结果如图7所示。可以看出，图中密集处为信号光子的部分，噪声光子离散分布于信号光子周围。在SNR较高的海岛植被覆盖区域与砂质潮间带区域，信号光子的分布形态清晰，但在SNR较低的近岸较浅水体区域和近岸较深水体区域，由于水体本身对太阳光的反射、后向散射等影响，使得水面的信号光子数偏少，且随着水深的增加激光在水体中的衰减也在加强，接收到的水底有效信号光子数大幅减少，使得该区域水面与水底信号光子的分布差异不明显。

下载: | 高精图 | 低精图

图7 实验区域光子数据有效信号范围

Fig. 7 Effective signal range of photon data in the experimental area

其次，使用网格统计法对提取有效信号范围后的光子原始观测数据进行粗去噪，横轴间隔 $Δ t$ 取值为5 ms，纵轴间隔 $Δ h$ 的取值为1 m，计算了海岛植被、砂质潮间带、近岸较浅水底与较深水体等4种类型实验数据的 $h_{a v e_m}$ ，将n小于 $h_{a v e_m}$ 的网格内光子数据剔除，得到了4种类型实验数据的粗去噪结果（图8），蓝色虚线框为网格统计法粗去噪后保留下的网格范围。可以发现粗去噪后的海岛植被覆盖区（图8（a））与砂质潮间带区（图8（b））远离信号中心的噪声光子得到了有效的剔除，且保留下的网格中均含有信号光子。但在SNR较低的水体区域（图8（c）与图8（d）），用网格统计的方法虽然能去除空气中的噪声光子，但水体中一些仅包含噪声光子的网格未能有效剔除。为了进一步分析粗去噪结果，将图8中的绿色框区域放大，结果如图9所示。可以看出，SNR较高的砂质潮间带区（图9（a））保留的网格内，在地物边缘附近位置还残存着离散的噪声光子，但相较于噪声光子，网格内的信号光子聚集性更明显。SNR较低的近岸较深水体区域（图9（b））不仅存在着类似于SNR较高区域的离散噪声光子，同时因受水中各种微小悬浮物及水体本身的后向散射等的影响，使得水体中网格内的噪声光子数量接近信号光子数量，造成粗去噪后在水体中仍保留了较多的噪声网格。同时激光能量在水体中随着水深的增加，海底地形的光子信号衰减严重，致使水底光子信号的SNR降低。

下载: | 高精图 | 低精图

图 8 4种不同实验区域粗去噪结果

Fig. 8 Rough denoising results of four different of experimental areas

下载: | 高精图 | 低精图

图9 粗去噪局部网格结果

Fig. 9 Coarse denoising local mesh result

4.1.2　精去噪结果与分析

对于粗去噪后保留网格内的光子数据，采用本文提出的改进LSC算法进行精去噪。初始阈值k的取值采用试错法确定，在SNR较高的植被与沙地区 $k$ 取值为5，在SNR较低的近岸较浅水域与近岸较深水域 $k$ 取值为10，精去噪结果如图10所示。可以看出SNR较高的海岛植被区（图10（a））和砂质潮间带区（10（b））信号光子附近离散分布的噪声光子得到了有效的剔除，SNR较低的近岸较浅水体区域（图10（c））和近岸较深水体区域（图10（d））的水体中仅包含噪声光子的网格内光子数据也得到了有效的剔除，同时很好的保留了水面与水底有效信号光子。

下载: | 高精图 | 低精图

图10 4种不同类型实验区域精去噪结果

Fig. 10 Fine denoising results of four different types of experimental areas

图11为图10中绿色框选的区域放大结果。精去噪后SNR较高的沙地区信号上下离散分布的噪声光子得到了有效的剔除，地物有效信号光子数据被完整保留。SNR较低的近岸较深水体区域（图11（b）），水体中残存的大量噪声点以及水面和水底地形信号附近离散的噪声光子均得到了有效剔除。

译

下载: | 高精图 | 低精图

图11 精去噪局部结果

Fig. 11 Fine denoising local results

4.2　ICESat-2去噪结果与分析

对于星载ICESat-2获取的海岛近岸水体光子剖面点云数据，首先采用基于网格统计的粗去噪方法将远离信号光子中心的噪声光子剔除，该试验区横纵轴网格的间隔大小分别设为100 m和5 m。统计每个网格内光子的数量，计算得到网格单元的平均光子数，并将小于平均光子数的网格作为噪声网格剔除，粗去噪结果如图12所示。从图中可看出，经网格统计法粗去噪后，水面以上以及水下远离信号光子中心的大部分离散的噪声光子得到有效剔除，保留下的网格内均包含有效的信号光子，但在水面和水底信号光子的附近，仍然存在一些离信号光子较近的噪声光子。

译

图12 珊瑚岛区域粗去噪结果

Fig. 12 Shanhu Island area rough denoising results

下载: 原图 | 高精图 | 低精图

采用改进LSC方法对保留下网格内的噪声光子进行进一步去噪处理，在该区域实验数据集中，初始阈值k设为5，椭圆长短半轴的比值为a∶b=5∶1。精去噪结果如图13所示，从图13中可看出，保留下的网格中内大部分距离信号光子较近的噪声光子得到有效剔除，同时较好地保留了水下地形的形态特征。

译

图13 珊瑚岛区域精去噪结果

Fig. 13 Shanhu Island area Fine denoising results

下载: 原图 | 高精图 | 低精图

4.3　3种不同去噪算法结果

为了评价本文提出的无人机载轻小型光子计数激光雷达原始测距数据去噪方法，采用加井岛的无人机载光子数据，对比经典的局部距离统计方法（

Zhang和Kerekes，2015）、LSC方法（Agyemang，2004）和改进LSC方法的去噪结果。采用人工目视标注的光子数据作为验证数据，利用混淆矩阵的F1-Score作为精度评价指标（Zhang和Kerekes，2015）。表3列出了不同类型下8组数据的去噪精度。其中，局部距离统计法、LSC方法和改进LSC方法的海岛植被覆盖区域F1-Score均值分别为91.86%、93.04%和94.64%，在砂质潮间带区域F1-Score均值分别为98.43%、97.86和98.96%。虽然在SNR较高区域每种方法的F1-Score均高于90%，但本文改进LSC方法F1-Score均高于局部距离统计法和LSC方法，表明改进LSC方法的去噪效果优于前面两种方法。在近岸较浅水体区域，水下地形测量中，3种方法的F1-Score均值分别为88.31%、89.96%和93.04%，在近岸较深水体区域3种方法的F1-Score均值分别为85.91%、87.19%和90.74%。改进的LSC去噪方法在低信噪比区域适用性明显优于其他两种方法。

译

表4 3种去噪方法的精度比较

Table 4 The evaluation results of algorithm accuracy

类型	序号	F1-Score/%
类型	序号	局部距离统计	LSC	改进LSC
植被区	D1	91.27	93.66	95.23
植被区	D2	92.44	92.42	94.04
沙地区	D3	98.25	97.56	98.75
沙地区	D4	98.61	98.17	99.18
近岸浅水区	D5	88.93	90.69	93.11
近岸浅水区	D6	87.68	89.22	92.97
近岸深水	D7	86.64	87.54	91.12
近岸深水	D8	85.19	86.84	90.36

下载: 导出CSV

下载: 下载表格图片

图14为加井岛近岸较深水体区域3种方法的去噪结果，从图14中可以看出局部距离统计法和LSC方法的去噪实验结果虽然有效剔除了水中的噪声光子，但由于该区域SNR较低，两种方法在去除水体中大量噪声光子的同时，也将水底部分有效的信号光子进行了误去除。改进的LSC去噪方法不仅将水体中大量噪声光子进行了有效的剔除，同时SNR较低的水下地形信号光子得到了较好的保留，从而提升了在低SNR的水下地形测深能力。

译

下载: | 高精图 | 低精图

图14 近岸深水区域3种方法去噪效果对比

Fig. 14 Comparison of denoising effects of three algorithms in nearshore deep water area

改进LSC方法、LSC方法和局部距离统计法的去噪方法在不同SNR下4种典型海岛及水下地形光子原始测距数据去噪后的F1-Score均值分别为94.34%、92.14%和91.13%。综上所述，改进LSC的去噪方法，能有效解决低空无人机载光子计数激光雷达系统的噪声去除问题。

译

5　结论

本文针对机载轻小型光子计数激光雷达原始观测数据的噪声去除问题，提出的一种基于改进LSC的去噪方法，并选取加井岛的海岛植被、砂质潮间带和近岸2种浅海水下地形等不同信噪比无人机载光子数据开展试验，F1-Score值均优于90%，为后续水深信号的提取提供可靠的数据基础。同时利用珊瑚岛星载光子计数激光雷达ICESat-2的剖面点云数据初步验证本文提出去噪方法的在不同平台数据上的适用性。

译

算法的去噪精度总体上已基本满足无人机载轻小型光子计数激光雷达的噪声处理要求，但光子数据受海岛测量中植被覆盖类型和潮间带类型、水深测量中的水质、深度和海底底质等多种因素影响，会造成本文提出的去噪方法中相关参数的自适应性不足，针对此类似情况拟在后续研究中添加更多的海岛及水下地形的数据类型，研究建立自适应性更强的参数确定方法。

译

参考文献（References）

Agyemang M. 2002. Local Sparsity Coefficient-Based Mining of Outliers. Windsor, Ontario: University of Windsor [百度学术]

Agyemang M. 2004. LSC-mine: algorithm for mining local outliers//IRMA International Conference. Hershey: IRM Press: 5-8 [百度学术]

Chen Y F, Le Y, Zhang D F, Wang Y, Qiu Z G and Wang L Z. 2021. A photon-counting LiDAR bathymetric method based on adaptive variable ellipse filtering. Remote Sensing of Environment, 256: 112326 [DOI: 10.1016/j.rse.2021.112326] [百度学术]

Duan X J. 2017. Imaging Technology Based on Time-Correlated Single Photon Counting. Xi’an: Xidian University [百度学术]

段雪洁. 2017. 基于时间相关单光子计数的成像技术. 西安: 西安电子科技大学 [百度学术]

Fouche D G. 2003. Detection and false-alarm probabilities for laser radars that use Geiger-mode detectors. Applied Optics, 42(27): 5388-5398 [DOI: 10.1364/ao.42.005388] [百度学术]

Li Z Y, Liu Q W and Pang Y. 2016. Review on forest parameters inversion using LiDAR. Journal of Remote Sensing, 20(5): 1138-1150 [百度学术]

李增元, 刘清旺, 庞勇. 2016. 激光雷达森林参数反演研究进展. 遥感学报, 20(5): 1138-1150 [DOI: 10.11834/jrs.20165130] [百度学术]

Liu H, Chen P, Mao Z H and Pan D L. 2020. Iterative retrieval method for ocean attenuation profiles measured by airborne lidar. Applied Optics, 59(10): C42-C51 [DOI: 10.1364/AQ.379406] [百度学术]

Liu S, Luan K F, Tan K and Zhang W G. 2021. Multi-type vegetation coverage tidal flat terrain filtering based on UAV LiDAR point cloud. Remote Sensing Technology and Application, 36(6): 1272-1283 [百度学术]

刘帅, 栾奎峰, 谭凯, 张卫国. 2021. 基于无人机LiDAR点云的多类型植被覆盖滩涂地形滤波. 遥感技术与应用, 36(6): 1272-1283 [DOI: 10.11873/j.issn.1004-0323.2021.6.1272] [百度学术]

Ma Y, Liu R, Li S, Zhang W H, Yang F L and Su D P. 2018. Detecting the ocean surface from the raw data of the MABEL photon-counting lidar. Optics Express, 26(19): 24752-24762 [DOI: 10.1364/OE.26.024752] [百度学术]

Markus T, Neumann T, Martino A, Abdalati W, Brunt K, Csatho B, Farrell S, Fricker H, Gardner A, Harding D, Jasinski M, Kwok R, Magruder L, Lubin D, Luthcke S, Morison J, Nelson R, Neuenschwander A, Palm S, Popescu S, Shum C K, Schutz B E, Smith B, Yang Y K and Zwally J. 2017. The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2): science requirements, concept, and implementation. Remote Sensing of Environment, 190: 260-273 [DOI: 10.1016/j.rse.2016.12.029] [百度学术]

Neumann T, Brenner A, Hancock D, Robbins J, Saba J, Harbeck K and Gibbons A. 2019. Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) Project: Algorithm Theoretical Basis Document (ATBD) for Global Geolocated Photons ATL03. National Aeronautics and Space Administration, Goddard Space Flight Center. Available online: https://icesat-2.gsfc.nasa.gov/sites/default/files/files/ATL03_05June2018.pdf (accessed on 21 May 2019). [百度学术]

Otsu N. 1979. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics, 9(1): 62-66 [DOI: 10.1109/tsmc.1979.4310076] [百度学术]

Popescu S C, Zhou T, Nelson R, Neuenschwander A, Sheridan R, Narine L and Walsh K M. 2018. Photon counting LiDAR: an adaptive ground and canopy height retrieval algorithm for ICESat-2 data. Remote Sensing of Environment, 208: 154-170 [DOI: 10.1016/j.rse.2018.02.019] [百度学术]

Shen G Y, Zheng T X, Li Z H, Wu E, Yang L, Tao Y L, Wang C H and Wu G. 2021. High-speed airborne single-photon LiDAR with GHz-gated single-photon detector at 1550 nm. Optics and Laser Technology, 141: 107109 [DOI: 10.1016/j.optlastec.2021.107109] [百度学术]

Shim H and Lee S. 2014. Hybrid exposure for depth imaging of a time-of-flight depth sensor. Optics Express, 22(11): 13393-13402 [DOI: 10.1364/oe.22.013393] [百度学术]

Wang C S, Li Q Q, Liu Y X, Wu G F, Liu P and Ding X L. 2015. A comparison of waveform processing algorithms for single-wavelength LiDAR bathymetry. ISPRS Journal of Photogrammetry and Remote Sensing, 101: 22-35 [DOI: 10.1016/j.isprsjprs.2014.11.005] [百度学术]

Wang X, Glennie C and Pan Z G. 2017. An adaptive ellipsoid searching filter for airborne single-photon lidar. IEEE Geoscience and Remote Sensing Letters, 14(8): 1258-1262 [DOI: 10.1109/lgrs.2017.2704917] [百度学术]

Wu Z S and Liu A N. 2002. Scattering of solar and atmospheric background radiation from a target. International Journal of Infrared and Millimeter Waves, 23(6): 907-917 [DOI: 10.1023/a:1015703418994] [百度学术]

Xia S B, Wang C, Xi X H, Luo S Z and Zeng H C. 2014. Point cloud filtering and tree height estimation using airborne experiment data of ICESat-2. Journal of Remote Sensing, 18(6): 1199-1207 [百度学术]

夏少波, 王成, 习晓环, 骆社周, 曾鸿程. 2014. ICESat-2机载试验点云滤波及植被高度反演. 遥感学报, 18(6): 1199-1207 [DOI: 10.11834/jrs.20144029] [百度学术]

Xie H, Ye D, Xu Q, Sun Y, Huang P Q, Tong X H, Guo Y L, Liu X S and Liu S J. 2022. A density-based adaptive ground and canopy detecting method for ICESat-2 photon-counting data. IEEE Transactions on Geoscience and Remote Sensing, 60: 4411813 [DOI: 10.1109/TGRS.2022.3176982] [百度学术]

Yang B S, Xu W X and Yao W. 2014. Extracting buildings from airborne laser scanning point clouds using a marked point process. GIScience and Remote Sensing, 51(5): 555-574 [DOI: 10.1080/15481603.2014.950117] [百度学术]

Zhang H H, Ding Y X and Huang G H. 2019. Photon counting laser bathymetry system. Infrared and Laser Engineering, 48(1): 0106002 [百度学术]

张河辉, 丁宇星, 黄庚华. 2019. 光子计数激光雷达测深系统. 红外与激光工程, 48(1): 0106002 [DOI: 10.3788/IRLA201948.0106002] [百度学术]

Zhang J S and Kerekes J. 2015. An adaptive density-based model for extracting surface returns from photon-counting laser altimeter data. IEEE Geoscience and Remote Sensing Letters, 12(4): 726-730 [DOI: 10.1109/LGRS.2014.2360367] [百度学术]

Zhu X X, Nie S, Wang C, Xi X H, Wang J S, Li D and Zhou H Y. 2021. A noise removal algorithm based on OPTICS for photon-counting LiDAR data. IEEE Geoscience and Remote Sensing Letters, 18(8): 1471-1475 [DOI: 10.1109/LGRS.2020.3003191] [百度学术]

摘要

Abstract

the high detection sensitivity also leads to the acquired photonic point cloud data with large background noise

a strong correlation between the signal-to-noise ratio and the type of ground objects

and large differences in the density distribution of photons

and the existing denoising algorithms cannot be well applied.

In this paper

a denoising method for raw photon observation data is proposed. First

the effective signal interval of the raw photon observation data is calculated based on the histogram statistics method

and then the data in the interval are coarsely denoised by the grid statistics method. Finally

the local sparse coefficient method is improved

the horizontal ellipse search is used to calculate the local sparse coefficient value of each photon data in the grid

and the method of maximum interclass variance is introduced to determine the separation threshold of noise photons and signal photons

which improves the original photon observation data. Denoising accuracy. Jiajing Island and the adjacent shallow sea terrain in Hainan Province are selected as the research area to verify the denoising algorithm proposed.

The results show that the average F1-score in the high signal-to-noise ratio areas

such as the island vegetation coverage area and the sandy intertidal zone

reaches 94.64% and 98.96%

respectively

and the average F1-score in the low signal-to-noise ratio area

such as the shallower and deeper water bodies near the coast

can also reach 93.04% and 90.74%

respectively. The overall F1-score is 94.34%

which can effectively remove most of the noise points and has strong adaptability to island vegetation

sandy land and underwater terrain of different depths with different signal-to-noise ratios.

In addition

this paper also selects the spaceborne ICESat-2 photon dataset of coral islands in the South China Sea

which further verifies the availability and applicability of the denoising algorithm proposed in this paper on spaceborne photonic point cloud data.

关键词

遥感光子计数无人机激光雷达局部稀疏系数最大类间方差法去噪

Keywords

remote sensingphoton countingUAVLiDARlocal sparse coefficientOTSUnoise removal

references

Agyemang M. 2002. Local Sparsity Coefficient-Based Mining of Outliers. Windsor, Ontario: University of Windsor

Agyemang M. 2004. LSC-mine: algorithm for mining local outliers//IRMA International Conference. Hershey: IRM Press: 5-8

Duan X J. 2017. Imaging Technology Based on Time-Correlated Single Photon Counting. Xi’an: Xidian University

段雪洁. 2017. 基于时间相关单光子计数的成像技术. 西安: 西安电子科技大学

Fouche D G. 2003. Detection and false-alarm probabilities for laser radars that use Geiger-mode detectors. Applied Optics, 42(27): 5388-5398 [DOI: 10.1364/ao.42.005388http://dx.doi.org/10.1364/ao.42.005388]

Li Z Y, Liu Q W and Pang Y. 2016. Review on forest parameters inversion using LiDAR. Journal of Remote Sensing, 20(5): 1138-1150

李增元, 刘清旺, 庞勇. 2016. 激光雷达森林参数反演研究进展. 遥感学报, 20(5): 1138-1150 [DOI: 10.11834/jrs.20165130http://dx.doi.org/10.11834/jrs.20165130]

Liu H, Chen P, Mao Z H and Pan D L. 2020. Iterative retrieval method for ocean attenuation profiles measured by airborne lidar. Applied Optics, 59(10): C42-C51 [DOI: 10.1364/AQ.379406http://dx.doi.org/10.1364/AQ.379406]

Liu S, Luan K F, Tan K and Zhang W G. 2021. Multi-type vegetation coverage tidal flat terrain filtering based on UAV LiDAR point cloud. Remote Sensing Technology and Application, 36(6): 1272-1283

刘帅, 栾奎峰, 谭凯, 张卫国. 2021. 基于无人机LiDAR点云的多类型植被覆盖滩涂地形滤波. 遥感技术与应用, 36(6): 1272-1283 [DOI: 10.11873/j.issn.1004-0323.2021.6.1272http://dx.doi.org/10.11873/j.issn.1004-0323.2021.6.1272]

Otsu N. 1979. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics, 9(1): 62-66 [DOI: 10.1109/tsmc.1979.4310076http://dx.doi.org/10.1109/tsmc.1979.4310076]

Shim H and Lee S. 2014. Hybrid exposure for depth imaging of a time-of-flight depth sensor. Optics Express, 22(11): 13393-13402 [DOI: 10.1364/oe.22.013393http://dx.doi.org/10.1364/oe.22.013393]

Xia S B, Wang C, Xi X H, Luo S Z and Zeng H C. 2014. Point cloud filtering and tree height estimation using airborne experiment data of ICESat-2. Journal of Remote Sensing, 18(6): 1199-1207

夏少波, 王成, 习晓环, 骆社周, 曾鸿程. 2014. ICESat-2机载试验点云滤波及植被高度反演. 遥感学报, 18(6): 1199-1207 [DOI: 10.11834/jrs.20144029http://dx.doi.org/10.11834/jrs.20144029]

Zhang H H, Ding Y X and Huang G H. 2019. Photon counting laser bathymetry system. Infrared and Laser Engineering, 48(1): 0106002

张河辉, 丁宇星, 黄庚华. 2019. 光子计数激光雷达测深系统. 红外与激光工程, 48(1): 0106002 [DOI: 10.3788/IRLA201948.0106002http://dx.doi.org/10.3788/IRLA201948.0106002]

文章被引用时，请邮件提醒。

提交

无人机激光雷达人工林林分高估测模型分析

无人机多光谱和LiDAR的红树林精细识别与生物量估算