前言
对森林高度进行持续的大尺度监测，对于估算森林碳排放、分析森林退化情况以及量化森林恢复措施的效果具有重要意义。星载激光雷达GEDI（Global Ecosystem Dynamics Investigation）通过发射激光脉冲，记录植被和地表反射回波的能量。利用森林冠层与地表回波的时间差，可以计算出森林冠层的高度。本文将介绍如何下载和读取GEDI_01B数据，并讲解使用GEDI数据提取冠层高度的相关理论知识。
数据下载与读取
如图所示，GEDI提供了不同级别的产品数据，其中01_B是经过地理校正的全波形数据，通常用于更灵活的处理和应用。此外，GEDI提供的树高和生物量产品数据也具有很高的应用价值。
1、GEDI01_B数据下载
下载地址：https://search.earthdata.nasa.gov/
先注册账号
在搜索窗口，搜索GEDI，并填上相应的筛选条件，最后进行下载
2、数据读取
在Python中使用h5py或pyGEDI这两个库可以实现01_B数据的读取，但是h5py更应用性更高，pyGEDI可以读取，并展示高程与回波能量的图，但是我怎么也没找到地理信息数据的读取与导出函数，单纯用pyGEDI应用很受限。
2.1 pyGEDI
在使用pyGEDI需要安装gdal库，但是即使安装好，你运行也会显示没有gdal模块，这时你需要在安装pyGEDI包的位置打开int文件。打开文件后在开头加入代码：from osgeo import gdal，就可以了。因为这个文件中是直接import gdal，这是无法引入的。
from pyGEDI import *
import matplotlib.pyplot as plt
fileh5_1B=r'C:\Users\YSY\Desktop\pyGEDI-master\pyGEDI-master\notebook\data\GEDI01_B_2019108080338_O01964_T05337_02_003_01_sub.h5'
h5_1B=getH5(fileh5_1B)
print(h5_1B)
shot_number=19640521100108408
rxwaveform,elevation = waveForm(shot_number,h5_1B)
beam=getBeam(shot_number,h5_1B)
print(beam)
fig = plt.subplots(figsize=(7,7))
plt.plot(rxwaveform,elevation, 'green')
plt.xlabel("Waveform Amplitude")
plt.ylabel("Elevation (m)")
plt.show()
2.2 h5py
import h5py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# d读取h5格式的波形数据
fileh5_1B = r'C:\Users\YSY\Desktop\pyGEDI-master\pyGEDI-master\notebook\data\GEDI01_B_2019108080338_O01964_T05337_02_003_01_sub.h5'
with h5py.File(fileh5_1B, 'r') as h5_file:
# 选择激光束
beam_group = h5_file['BEAM0101']
# 选择光斑
shot_number = 19640521100108408
shots = beam_group['shot_number'][:]# 可以提前打印看光斑号
# 找到先对应的光斑
shot_index = np.where(shots == shot_number)[0][0]
# 提取地理位置和回波能量
rx_sample_count = beam_group['rx_sample_count'][shot_index]
rx_start_index = beam_group['rx_sample_start_index'][shot_index]
rx_waveform = beam_group['rxwaveform'][rx_start_index:rx_start_index + rx_sample_count]
latitudes = beam_group['geolocation']['latitude_bin0'][shot_index]
longitudes = beam_group['geolocation']['longitude_bin0'][shot_index]
elevation_bin0 = beam_group['geolocation']['elevation_bin0'][shot_index]
elevation_lastbin = beam_group['geolocation']['elevation_lastbin'][shot_index]
# Generate elevation range
elevation_range = np.linspace(elevation_lastbin, elevation_bin0, rx_sample_count)
# 将地理位置和回波能量导出
data = {
'Latitude': [latitudes] * len(rx_waveform),
'Longitude': [longitudes] * len(rx_waveform),
'Elevation (m)': elevation_range,
'Waveform Amplitude': rx_waveform
}
df = pd.DataFrame(data)
df.to_csv('gedi_waveform_data.csv', index=False)
# Plot waveform amplitude vs. elevation
plt.figure(figsize=(7, 7))
plt.plot(rx_waveform, elevation_range, color='green')
plt.xlabel("Waveform Amplitude")
plt.ylabel("Elevation (m)")
plt.title(f"Waveform Amplitude vs Elevation for Shot {shot_number}")
plt.gca().invert_yaxis() # Invert y-axis for correct elevation view
plt.show()
这几天我翻阅资料有以下几点需要在注意：
- 01_B数据需要根据“stable_return_flag”和“degrade”的值进行质量过滤，通常会删除这些值为0的行，以确保数据质量；
- GEDI光斑的地理位置信息通常具有15至20米的误差，因此在应用中需考虑这种空间偏差；
- BEAM0101、BEAM0110、BEAM1000、BEAM1011是全功率脉冲（Full power beam），具有更强的穿透能力，由两台激光器生成，并快速偏转一次形成。而BEAM0000、BEAM0001、BEAM0010、BEAM0011是覆盖光束（Coverage beam），由一台激光器分成两束，并以1.5毫弧度（mrad）的角度偏转，生成四条光束。
结果如下图所示，高度变化相对正常，每个bin之间的差距约为15厘米，这与GEDI记录时间间隔为1ns相对应。
3、冠层高度提取
在我撰写的第一篇推文《利用GEDI数据反演森林冠层高度》中，已经介绍了如何在GEE平台上使用GEDI数据与被动光学影像（如Sentinel-2）融合，实现大范围高分辨率的森林冠层高度制图。那时我还不熟悉Markdown的使用，代码是以图片形式上传的，也缺少理论的详细解释，当时心情比较紧张。提取冠层高度时，使用的是GEDI_02A数据，关键在于我们应使用相对高度中的哪一个作为森林冠层高度值。根据《Mapping global forest canopy height through integration of GEDI and Landsat data》一文，参考基于机载点云计算的冠层高度，使用RH95高度作为冠层高度会更加准确。
那么，为什么不直接使用提取的RH100，而是要将高度乘以一个百分比呢？
我认为主要有两个原因：第一是脉冲宽度引起的误差，即使激光打到平地，回波能量的波形仍会显示大约四米的高度；第二是坡度引起的误差，坡度会导致土壤层的回波宽度增加，从而影响测量精度。
因此，为了减少这些误差，我们通常会将波形高度乘以一个小于1的比值（如0.95，RH95），作为冠层高度。注意一点使用的波形数据对应位置的森林林木状态应该未落叶，否则高度会被低估。
