探索 LiDAR 栅格数据中的不确定性(python)

主要是针对阴天和晴天获取的雷达数据的不确定性进行查看，预期是这两次收集之间没有任何变化。

内容包括：

• 1 读取数据
• 2 区域整体不确定性
• 3 植被区不确定性
• 4 非植被区不确定性

1 读取数据

chm1_fname = os.path.join(chm_dir,'2016_PRIN_1_CHM.tif')
dsm1_fname = os.path.join(dem_dir,'2016_PRIN_1_DSM.tif')
dtm1_fname = os.path.join(dem_dir,'2016_PRIN_1_DTM.tif')

chm2_fname = os.path.join(chm_dir,'2016_PRIN_2_CHM.tif')
dsm2_fname = os.path.join(dem_dir,'2016_PRIN_2_DSM.tif')
dtm2_fname = os.path.join(dem_dir,'2016_PRIN_2_DTM.tif')

chm1_dataset = rio.open(chm1_fname)
dsm1_dataset = rio.open(dsm1_fname)
dtm1_dataset = rio.open(dtm1_fname)

chm2_dataset = rio.open(chm2_fname)
dsm2_dataset = rio.open(dsm2_fname)
dtm2_dataset = rio.open(dtm2_fname)

2 区域整体不确定性

查看这两天数据的差异，直接两个 DSM 数据相减：

dsm1_data = dsm1_dataset.read(1)
dsm2_data = dsm2_dataset.read(1)
diff_dsm_array = np.subtract(dsm1_data,dsm2_data)
diff_dsm_array_mean = np.mean(diff_dsm_array)
diff_dsm_array_std = np.std(diff_dsm_array)
print('Mean difference in DSMs: ',round(diff_dsm_array_mean,3),'m')
print('Standard deviation of difference in DSMs: ',round(diff_dsm_array_std,3),'m')

# Mean difference in DSMs:  0.019 m
# Standard deviation of difference in DSMs:  0.743 m

均值为零 (0.019 米) 附近，表明这两天数据几乎没有系统偏差。但数据的标准差相当高，为 0.743 米。通常，期望激光雷达数据误差低于 0.15 米。

看一下 DSM 差异的直方图：

直方图有长长的尾部，掩盖了中心附近的分布。为了限制直方图的 x 轴限制，将使用刚刚计算出的均值和标准差。由于数据似乎呈正态分布，可以通过包括均值上下 2 个标准差来将直方图限制在 95% 的数据范围内。

直方图显示了 DSM 差异的广泛变化，其中 95% 的极限值在约 ±1.5 米。

查看 DSM 差异的空间分布，并将色条的数值约束在 95% 的观测值范围内：

从图上看，误差似乎存在空间模型，接下来查看 DTM 差异的统计信息：

dtm1_data = dtm1_dataset.read(1)
dtm2_data = dtm2_dataset.read(1)
diff_dtm_array = np.subtract(dtm1_data,dtm2_data)

diff_dtm_array_mean = np.mean(diff_dtm_array)
diff_dtm_array_std = np.std(diff_dtm_array)
print('Mean difference in DTMs: ',round(diff_dtm_array_mean,3),'m')
print('Standard deviation of difference in DTMs: ',round(diff_dtm_array_std,3),'m')
# Mean difference in DTMs:  0.014 m
# Standard deviation of difference in DTMs:  0.102 m

差异的整体幅度小于 DSM，但误差的相同空间模式是明显的。

3 植被区不确定性

绘制冠层高度模型 (CHM)。在 CHM 中，树冠高度被表示出来，所有地面像素的海拔高度为零。使用一个颜色条，其中地面显示为浅绿色，最高植被显示为深绿色：

从 CHM 来看，误差模式的空间分布遵循植被的位置。

提取 DSM 差异中对应植被位置的像素，计算其均值和标准差，并绘制相关直方图。在显示直方图之前，将从 DSM 差异移除非数值像素以及从 CHM 中移除非零像素。

chm1_data = chm1_dataset.read(1)
# 植被数字表面模型差异
diff_dsm_array_veg_mean = np.nanmean(diff_dsm_array[chm1_data!=0.0])
diff_dsm_array_veg_std = np.nanstd(diff_dsm_array[chm1_data!=0.0])
print('Mean difference in DSMs on veg points: ',round(diff_dsm_array_veg_mean,3),'m')
print('Standard deviations of difference in DSMs on veg points: ',round(diff_dsm_array_veg_std,3),'m')
# Mean difference in DSMs on veg points:  0.072 m
# Standard deviations of difference in DSMs on veg points:  1.405 m

plt.figure();
# 移除DSM差异的非数值像素
diff_dsm_array_nodata_removed = diff_dsm_array[~np.isnan(diff_dsm_array)]
chm_dsm_nodata_removed = chm1_data[~np.isnan(diff_dsm_array)]

dsm_diff_veg_vmin = diff_dsm_array_veg_mean-2*diff_dsm_array_veg_std
dsm_diff_veg_vmax = diff_dsm_array_veg_mean+2*diff_dsm_array_veg_std

# 移除CHM的零值(地面部分)
plt.hist(diff_dsm_array_nodata_removed[chm_dsm_nodata_removed!=0.0],100,range=[dsm_diff_veg_vmin, dsm_diff_veg_vmax])
plt.title('Histogram of PRIN DSM Difference in Vegetated Areas')
plt.xlabel('Height Difference(m)'); plt.ylabel('Frequency');

结果表明平均差异接近零，但标准差高达 1.405 米！由于数字表面模型(DSM)代表树冠顶部，这表明冠层高度的不确定性水平。

植被区 DTM 差异分布：

diff_dtm_array_veg_mean = np.nanmean(diff_dtm_array[chm1_data!=0.0])
diff_dtm_array_veg_std = np.nanstd(diff_dtm_array[chm1_data!=0.0])
print('Mean difference in DTMs on vegetated pixels: ',round(diff_dtm_array_veg_mean,3),'m')
print('Standard deviations of difference in DTMs on vegetated pixels: ',round(diff_dtm_array_veg_std,3),'m')
# Mean difference in DTMs on vegetated pixels:  0.023 m
# Standard deviations of difference in DTMs on vegetated pixels:  0.163 m

平均差几乎为零 (0.023 米)，变化小于 DSM 变化 (0.163 米)。

尽管变化有所减少，但仍然比预期的大。这是因为植被下方可能没有太多激光能量到达地面，而到达地面的点可能返回的信号较弱。点的稀疏性导致在较大距离上进行表面插值，这可能会忽略地形的变化。由于 LIDAR 点的分布每天不同，这导致了不同的地形表示和地表的不确定性。这表明在有植被存在时，LiDAR DTM 的精度会降低。

4 非植被区不确定性

仅查看地面点 (CHM = 0) 上的 DTM 差异。

diff_dtm_array_ground_mean = np.nanmean(diff_dtm_array[chm1_data==0.0])
diff_dtm_array_ground_std = np.nanstd(diff_dtm_array[chm1_data==0.0])
print('Mean difference in DTMs on ground points: ',round(diff_dtm_array_ground_mean,3),'m')
print('Standard deviations of difference in DTMs on ground points: ',round(diff_dtm_array_ground_std,3),'m')

# Mean difference in DTMs on ground points:  0.011 m
# Standard deviations of difference in DTMs on ground points:  0.069 m

平均差仅为 0.011 米，变化为 0.068 米。

表明 LiDAR 系统的不确定性 (约 0.15 米) 仅在裸露、开阔、硬质表面时才有效。当存在植被时，不能期望 LiDAR 的精度如此之高。量化树冠顶部尤其困难，并且可能导致任何给定像素的不确定性超过 1 米。