本发明涉及建筑高度估算领域，具体涉及一种顾及不同地形场景下的城市建筑高度提取方法。
背景技术：
1、城市用地仅占地球陆地表面积的3％左右，却造成了全球能源消耗的60-80％和碳排放的75％。若要厘清城市化进程对城市环境的影响，城市三维结构信息至关重要。由于数据获取等限制，现有城市化研究多探讨的是城市水平扩张问题。相比之下，揭示城市垂直扩张分布的研究较少，可见大尺度高分辨率的城市建筑高度数据是实现城市可持续发展等重大需求较紧缺的关键参数。另外，已有研究表明诸多城市物理量与建筑高度存在显著函数关系，例如城市人口分布、物资存量分配、城市热岛效应等。在此背景下，可将城市建筑高度估算作为当下研究的重要课题。
2、迅速发展的遥感技术为城市建筑高度制图提供了丰富的数据。从获取城市建筑三维结构信息的精细程度，可将其划分为两方面讨论：区域尺度的城市建筑高度反演、精细尺度的城市建筑高度反演。在区域和全球尺度上，sar估计不同空间分辨率建筑高度的可行性和有效性。然而，对于异质性较高的城市区域，不同建筑物的sar后向散射信号难以分离、互相干扰，导致建筑高度估计存在明显的不确定性。精细尺度下的城市建筑高度反演，目前可划分为以下四类：
3、1)基于高分辨率光学影像捕捉建筑物阴影长度来估算。liasis and stavrou的研究表明，基于高分辨率光学图像，可根据建筑物的阴影长度及其与太阳和传感器之间的几何关系来推算建筑高度。其算法如下：a.通过空间、光谱分析来自定义滤波器，增强建筑物阴影区域；b.基于活动轮廓模型(active contour models)以图像分割方式提取阴影，辅以形态学操作剔除阴影杂质；c.结合建筑物阴影长度、预定义或估计的太阳高度角定义几何关系推算建筑高度。但该研究是假设研究区内建筑物方位角一致条件下，仅限于小范围区域。zhao et al.认为当建筑物方位角小于180°，部分建筑物阴影将存在叠加的情况。于是在搭建建筑物阴影长度、建筑高度之间的定量关系时，纳入了根据建筑物方位角计算的校定参数。综上来看，基于阴影的方法很依赖于阴影分割的准确性，但复杂城市景观中阴影不可避免受不同地物遮挡，进而出现建筑物阴影分布不连续的情况，导致城市建筑高度被低估。
4、2)基于街景图像(svis)来估算：google maps可高效率、低成本的获取城市街道场景图像(street view images,svis)。yan and huang基于svis通过消失点/灭点(vanishing point)检测、图像线段检测与语义分割等过程分割出图像中所有的垂直线段，最后筛选出最佳的垂直线段来估算建筑高度。但在复杂街道场景下，建筑物之间、建筑物与其他地物之间存在遮挡，导致识别的垂直线段很难匹配到实际对应的建筑物上。与上不同，zhao et al.提出了新的建筑高度估算策略，同时采用了从沿街道方向的街道场景图像和正视建筑物的图像，并且改识别建筑物垂直线段为识别建筑物角点、屋顶线，联合相机位置来估算城市建筑高度。此方法在复杂街道场景下估算建筑高度稳定性更高，但问题是同一建筑物的街景图像的可用性有限。
5、3)基于高分辨率立体图像来估算：随着立体影像(zy-3/gf-7)和深度神经网络的迅速发展，很多学者开始从立体图像中生成高质量的视差数据来捕捉建筑物几何特征。caoand huang引入了可同时获取多光谱和多视角图像的zy-3数据，并提出了多光谱、多视角、多任务的深度网络(m3net)用于建筑高度估计。其中，多任务网络集成了多视角图像计算建筑高度的网络、多光谱图像计算建筑高度的网络和识别建筑物脚印的随机森林模型。最后，级联各任务的输出结果来预测最终的建筑高度。chen et al.基于gf-2立体图像结合深度网络stereonet和attention u-net，分别提取了城市dsm和建筑物屋顶数据，最后以形态学操作方法去除dsm中的地形信息来得到表征建筑高度的非地形高度，即归一化dsm(normalized dsm,ndsm)。明显的不足是上述方法未考虑城市坡度，虽然huang et al.基于形态学操作理论提出了dsm地形校正方法，但在地形复杂的城市场景下仍有待进一步探索。总体上立体图像在城市建筑高度估算上，精度非常可靠，但其数据获取困难、且受云量干扰严重等问题严重阻碍了其大尺度应用。
6、4)地基、车载和机载等激光雷达数据虽然具备了获取城市建筑物高精准三维结构信息的能力，但成本过高，不便于推广至大尺度应用。
7、综上可见，区域尺度的城市建筑高度反演研究，其结果数据无法映射每个建筑高度信息且缺乏高精度的建筑高度样本；然而精细尺度的城市建筑高度反演研究，则是不同数据源存在诸多不足，很难推广至区域或全球尺度计算。在未来实现“即见城市又见建筑”尺度的城市建筑高度反演，目前已有研究从gedi自带高度数据出发，进行初步探索。目标是从密度更高、精度更高的icesat-2/atlas光子点云中获取精准的城市建筑高度样本，以此支撑未来“即见城市又见建筑”尺度的城市建筑高度反演。目前，lao et al.已基于icesat-2/atlas探究了相应的城市建筑高度提取算法，但在城市光子点云去噪方面仍存在不足、以及并未考虑地形起伏对建筑高度估算的影响。
技术实现思路
1、本发明的目的在于：针对目前存在的问题，提供了一种顾及不同地形场景下的城市建筑高度提取方法，基于icesat-2/atl03光子点云数据提出一种参数自适应的城市光子点云分类算法，可在不同地形场景下获取更精准的城市建筑高度数据。
2、本发明的技术方案如下：
3、一种顾及不同地形场景下的城市建筑高度提取方法，包括以下步骤：
4、从原始光子中识别建筑光子和地面光子，基于局部参数统计原理提取光子局部特征acdi来自适应的识别建筑光子及地面光子；
5、从建筑光子和地面光子中识别地面光子，从城市地面光子在水平和垂直两个方向上的分布特征，来准确识别地面光子；
6、识别建筑光子及建筑高度估算，在获取地面曲线后，即获取地面曲线3m以上的信号光子为建筑光子，并对建筑光子进行后处理；
7、结合地面曲线和建筑光子数据可估算准确的建筑高度参数。
8、进一步的，所述从原始光子中识别建筑光子和地面光子，基于局部参数统计原理提取光子局部特征acdi来自适应的识别建筑光子及地面光子具体包括以下步骤：
9、icesat-2/alt03光子点云数据读取；
10、基于地形自适应的高程缓冲带实现粗去噪；
11、光子点云等距分段，且各分段作边缘缓冲和镜像处理；
12、基于局部参数统计提取建筑光子及地面光子；
13、精确识别建筑光子和地面光子。
14、进一步的，所述基于局部参数统计提取建筑光子及地面光子包括以下步骤：
15、步骤一、设置数据处理尺度及规则化：以1km为数据处理的基本尺度；通过acdi识别建筑光子和地面光子时，需凭借移动窗口法在相对沿轨方向作等距分段，并对每个分段进行缓冲和镜像处理；
16、步骤二、自适应局部参数统计及acdi计算：在每个分段构建高程频率直方图并获取最高直方柱，并以ransac对相应的光子数据作线性拟合；
17、获取各点至拟合线的平均距离，定义其值为和的长度，定义τminor的长度为短轴的十倍；
18、将光子点云数据置于x-y-z三维空间内，椭球体τminor依次绕z轴、x轴旋转变换，两次变换的坐标计算如公式如下所示：
19、
20、其中，(x，y，z)是单光子的三维坐标，其分别表示dist_ph_cross、dist_ph_along和h_ph三个位置参数；而(x′，y′，z′)是单光子(x，y，z)绕z轴旋转α°后坐标；(x″，y″，z″)则是单光子(x′，y′，z′)绕z轴旋转β°后坐标；
21、获取单光子p在三维空间内的最大邻域及垂直邻域，获取垂直邻域的光子数ncvertical和最大邻域的光子数nchorizonta并结合以下公式计算单光子p的vhdrp指标：
22、
23、其中，nc(θ0°)，...，nc(θ90°)，...，nc(θ180°)是不同角度θ下的单光子的邻域点数；ncvertical是垂直邻域的光子数；nchorizontal是最大邻域的光子数，nr是ncvertical和nchorizontal中最大值与最小值的比值；
24、获取最大邻域内所有的邻域光子，并通过ransac对其线性拟合，并计算每个邻域光子至拟合线的距离接着即可结合以下公式计算单光子p的phalp指标：
25、
26、其中，是最大邻域中邻域光子qi和邻域中心光子p的欧式距离；distmean则是所有邻域光子至邻域中心光子的平均距离；
27、综合vhdrp和phalp两个光子邻域特征值，根据以下公式计算出单光子p的acdip：
28、acdip＝phalp+vhdrp；
29、基于该指标从总体光子中识别出建筑光子和地面光子；
30、步骤三、自适应确定去噪阈值：对acdip阈值的训练样本范围定在[-1，2]之间，以0.05为跨度依次进行点云分类，在不同阈值i下将获取信号集和噪声集，经pca计算后获取前三主成分的方差；接着通过以下公式计算信号集和噪声集各主成分方差的标准差和当和的绝对差值最大时对应的i即为acdi最佳分割阈值，
31、
32、
33、其中，var(pcai)是第i主成分的方差；varmean是主成分分析的前三主成分方差的平均值；std是前三主成分方差的标准差；和是当acdi分类阈值为i时，被分类为信号光子和噪声光子的数据，通过主成分分析后前三主成分方差的标准差；另外，当和差值达到最大时，以此时的i为acdi的最佳分割阈值。
34、进一步的，所述粗去噪包括以下步骤：
35、以1km为icesat-2剖面点云数据处理尺度，根据50m跨度的移动窗口在沿轨方向建立等距分段，在各分段内，在高程范围内划分10个等距分箱并作高程频率统计；获取频数最高的分箱，并以在其上下方各取一分箱的位置为缓冲区，除此之外的分箱将视作噪声被剔除。
36、进一步的，所述从建筑光子和地面光子中识别地面光子具体包括以下步骤：
37、移动窗口下通过百分位统计法识别初始地面光子；
38、基于地面光子的水平分布特征对地面光子初过滤；
39、区域生长识别所有地面光子，同时以mdnn突变时作为停止区域生长标志；
40、基于地面光子的垂直分布特征对地面光子终过滤，得到精确的地面光子。
41、进一步的，所述从建筑光子和地面光子中识别地面光子具体包括以下步骤：
42、初始地面光子提取和去噪：沿初始地面光子水平分布方向进行迭代ransac线性拟合，并以拟合点至拟合线平均距离、拟合线斜率阈值作为分割条件，将初始地面光子聚类为若干线状点簇；过滤线状点簇，包括全局过滤，即根据线状点簇的拟合斜率进行过滤；以及局部过滤；将每个线性点簇当作一个泵，泵可将水流从地势高的地方输送至地势低的地方，也可在一定增压能力下将水流从地势低的地方输送至地势高的地方，泵的增压能力通过地势高程差来衡量；那么可在一定高程差阈值下去除伪地面光子；以连续分布的三个线性点簇为步长，获取中间点簇相对两侧点簇的高度差，相对两侧起伏较大的点簇被滤除，该过程迭代进行直至相邻线性点簇的相对位置趋于平稳；
43、区域生长完善地面光子识别：对去噪后的初始地面光子构建kd-tree索引，依次获取单光子的最大邻域方向，结合椭圆结构对单光子进行区域生长，以完善地面光子识别；通过以下公式计算mdnn，并设置mdnn阈值，在区域生长中若mdnn出现突变或零变化则停止区域生长，
44、
45、其中，n表示椭圆生长域内光子数据；表示第i组最邻近点对pi和qi间的欧式距离；
46、地面光子二次去噪：需要对已识别的地面光子进行分段，在适宜窗口下进行该去噪处理；在每个窗口内，根据识别的地面光子构建高程频率直方图，将所有的直方柱聚类分组，实现地面光子的点簇分类；定义每个地面点簇为待判断的目标点簇，并根据目标点簇的横坐标范围获取原始数据中的对应光子，构建该部分光子数据的高程频率直方图并对直方柱进行聚类分组，将定位目标点簇在分组中的位置；根据目标点簇在原始数据高程频率直方图中的相对位置来判断目标点簇是否为地面光子；
47、地面曲线获取：由每个窗口中获取最终的地面光子数据，并联合傅里叶变换和lowess拟合地面曲线，结合b样条插值算法将地面曲线插值至原始光子分辨率。
48、进一步的，所述结合地面曲线和建筑光子数据可估算准确的建筑高度参数具体包括以下步骤：
49、联合lowess拟合和b样条插值获取信号光子分辨率的地面曲线，得到精准的地面曲线；
50、获取地面曲线3m以上的信号光子为建筑光子，对建筑光子后处理：以区域生长完善建筑光子，以knn剔除离群伪建筑光子，得到精准的建筑光子；
51、结合精准的地面曲线和建筑光子数据可估算准确的建筑高度参数。
52、进一步的，所述对建筑光子进行后处理包括以下步骤：
53、获取建筑光子和地面光子集合、地面曲线，定义地面曲线3m以上的光子为初始建筑光子；
54、以mdnn作为终止区域生长的区域生长法来完善建筑光子；
55、通过knn算法来检测离群光子，并采用otsu算法确定阈值进行剔除，knn算法的最邻近点数为5个。
56、进一步的，所述估算准确的建筑高度参数包括以下步骤：
57、获取建筑光子和地面曲线后，利用b样条插值法将地面曲线插值至建筑光子分辨率；将每个建筑光子将匹配对应位置上的地面光子，计算两者的高程差为建筑高度；根据实景三维模型量算的高精度建筑高度数据为参照，然后分别计算均方根误差、平均绝对误差和决定系数三个指标来评估不同场景下的建筑高度估算精度。
58、与现有的技术相比本发明的有益效果是：
59、1、一种顾及不同地形场景下的城市建筑高度提取方法，所提出方法框架可准确识别不同城市地形场景下的建筑光子和地面光子。通过量化不同场景下建筑高度的估算精度，与参考建筑高度具备较高一致性；平坦地形下，rmse在0.230m-0.315m之间，mae在0.161m-0.243m之间；起伏地形下rmse在0.912m-1.424m之间，mae在0.647m-0.871m之间。与同行研究相比，该方法在不同场景下表现出更良好的适应性；
60、2、一种顾及不同地形场景下的城市建筑高度提取方法，提出一种基于单光子局部特征acdi的自适应信号光子提取方法，并且对不同的城市场景作了验证(包括日间-起伏地形、日间-平坦地形、夜间-起伏地形、夜间-平坦地形)，相比lao et al.的研究，所提出的算法能够更准确的捕捉icesat-2/atlas沿轨方向的信号光子(即建筑光子和地面光子)数据。其中，acdi是在自适应椭球体结构下，获取单光子邻域参数phal和vhdr来计算，最后在pca辅助下自适应获取分割信号光子和其他光子(即噪声光子和其他地物光子)的阈值，与otsu算法对比发现，该方法的阈值分割算法效果更佳；
61、3、一种顾及不同地形场景下的城市建筑高度提取方法，针对起伏地形下，城市地形识别不稳定的问题。根据地面光子的水平分布特征和垂直分布特征提出一种简洁、高效的地面识别方法。相比同行lao et al.的研究，该方法能更稳定的解决城市地面识别问题。其方法划分三个过程：a.通过移动窗口获取局部高程最低点为初始地面光子；b.借鉴li etal.研究中根据水泵布设准则剔除伪地面光子的思路，对初始地面光子作了粗去噪；c.将未剔除干净的伪地面光子(源自建筑光子)通过区域生长从“离群点”延伸至“带状建筑光子”，接着基于高程频率直方图，从光子垂直分布方向划分出地面光子簇和伪地面光子簇(即“带状建筑光子”)。结果表明，该方法显然是通用于不同地形场景下的地面光子识别；
62、4、一种顾及不同地形场景下的城市建筑高度提取方法，在获取准确的地面光子后，联合傅里叶变换和lowess算法拟合了平稳的地形曲线，接着提取地面曲线3m以上的信号光子为建筑光子。同时，为减少建筑光子的缺失，以区域生长和knn方法分别对建筑光子进行补齐及剔除离群噪声光子。
1.一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述从原始光子中识别建筑光子和地面光子，基于局部参数统计原理提取光子局部特征acdi来自适应的识别建筑光子及地面光子具体包括以下步骤：
3.根据权利要求2所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述粗去噪包括以下步骤：
4.根据权利要求1或2所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述基于局部参数统计提取建筑光子及地面光子包括以下步骤：
5.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述从建筑光子和地面光子中识别地面光子具体包括以下步骤：
6.根据权利要求1或5所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述从建筑光子和地面光子中识别地面光子具体包括以下步骤：
7.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述结合地面曲线和建筑光子数据可估算准确的建筑高度参数具体包括以下步骤：
8.根据权利要求1或7所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述对建筑光子进行后处理包括以下步骤：
9.根据权利要求1或7所述的一种顾及不同地形场景下的城市建筑高度提取方法，其特征在于，所述估算准确的建筑高度参数包括以下步骤：
