Nothing Special   »   [go: up one dir, main page]

CN117932333A - 一种顾及不同地形场景下的城市建筑高度提取方法 - Google Patents

一种顾及不同地形场景下的城市建筑高度提取方法 Download PDF

Info

Publication number
CN117932333A
CN117932333A CN202311775357.2A CN202311775357A CN117932333A CN 117932333 A CN117932333 A CN 117932333A CN 202311775357 A CN202311775357 A CN 202311775357A CN 117932333 A CN117932333 A CN 117932333A
Authority
CN
China
Prior art keywords
photons
ground
building
photon
urban
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.)
Pending
Application number
CN202311775357.2A
Other languages
English (en)
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.)
Yunnan Normal University
Original Assignee
Yunnan Normal University
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 Yunnan Normal University filed Critical Yunnan Normal University
Priority to CN202311775357.2A priority Critical patent/CN117932333A/zh
Publication of CN117932333A publication Critical patent/CN117932333A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00Measuring arrangements characterised by the use of optical techniques
    • G01B11/02Measuring arrangements characterised by the use of optical techniques for measuring length, width or thickness
    • G01B11/06Measuring arrangements characterised by the use of optical techniques for measuring length, width or thickness for measuring thickness ; e.g. of sheet material
    • G01B11/0608Height gauges
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/10Pre-processing; Data cleansing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种顾及不同地形场景下的城市建筑高度提取方法,包括以下步骤:从原始光子中识别城市信号光子(即建筑光子和地面光子)。基于自适应局部参数统计方法提取光子局部特征ACDI来识别城市信号光子;从建筑光子和地面光子中提取地面光子。从城市地面光子在水平和垂直两个方向上的分布特征,来准确识别地面光子;识别建筑光子及建筑高度估算。根据地面光子拟合准确的地面曲线即可提取建筑光子,接着对建筑光子进行后处理;最后,结合地面曲线和建筑光子数据可估算准确的建筑高度参数。可在不同地形场景下获取更精准的城市建筑高度数据。

Description

一种顾及不同地形场景下的城市建筑高度提取方法
技术领域
本发明涉及建筑高度估算领域,具体涉及一种顾及不同地形场景下的城市建筑高度提取方法。
背景技术
城市用地仅占地球陆地表面积的3%左右,却造成了全球能源消耗的60-80%和碳排放的75%。若要厘清城市化进程对城市环境的影响,城市三维结构信息至关重要。由于数据获取等限制,现有城市化研究多探讨的是城市水平扩张问题。相比之下,揭示城市垂直扩张分布的研究较少,可见大尺度高分辨率的城市建筑高度数据是实现城市可持续发展等重大需求较紧缺的关键参数。另外,已有研究表明诸多城市物理量与建筑高度存在显著函数关系,例如城市人口分布、物资存量分配、城市热岛效应等。在此背景下,可将城市建筑高度估算作为当下研究的重要课题。
迅速发展的遥感技术为城市建筑高度制图提供了丰富的数据。从获取城市建筑三维结构信息的精细程度,可将其划分为两方面讨论:区域尺度的城市建筑高度反演、精细尺度的城市建筑高度反演。在区域和全球尺度上,SAR估计不同空间分辨率建筑高度的可行性和有效性。然而,对于异质性较高的城市区域,不同建筑物的SAR后向散射信号难以分离、互相干扰,导致建筑高度估计存在明显的不确定性。精细尺度下的城市建筑高度反演,目前可划分为以下四类:
1)基于高分辨率光学影像捕捉建筑物阴影长度来估算。Liasis and Stavrou的研究表明,基于高分辨率光学图像,可根据建筑物的阴影长度及其与太阳和传感器之间的几何关系来推算建筑高度。其算法如下:a.通过空间、光谱分析来自定义滤波器,增强建筑物阴影区域;b.基于活动轮廓模型(Active contour models)以图像分割方式提取阴影,辅以形态学操作剔除阴影杂质;c.结合建筑物阴影长度、预定义或估计的太阳高度角定义几何关系推算建筑高度。但该研究是假设研究区内建筑物方位角一致条件下,仅限于小范围区域。Zhao et al.认为当建筑物方位角小于180°,部分建筑物阴影将存在叠加的情况。于是在搭建建筑物阴影长度、建筑高度之间的定量关系时,纳入了根据建筑物方位角计算的校定参数。综上来看,基于阴影的方法很依赖于阴影分割的准确性,但复杂城市景观中阴影不可避免受不同地物遮挡,进而出现建筑物阴影分布不连续的情况,导致城市建筑高度被低估。
2)基于街景图像(SVIs)来估算:Google maps可高效率、低成本的获取城市街道场景图像(Street View Images,SVIs)。Yan and Huang基于SVIs通过消失点/灭点(Vanishing point)检测、图像线段检测与语义分割等过程分割出图像中所有的垂直线段,最后筛选出最佳的垂直线段来估算建筑高度。但在复杂街道场景下,建筑物之间、建筑物与其他地物之间存在遮挡,导致识别的垂直线段很难匹配到实际对应的建筑物上。与上不同,Zhao et al.提出了新的建筑高度估算策略,同时采用了从沿街道方向的街道场景图像和正视建筑物的图像,并且改识别建筑物垂直线段为识别建筑物角点、屋顶线,联合相机位置来估算城市建筑高度。此方法在复杂街道场景下估算建筑高度稳定性更高,但问题是同一建筑物的街景图像的可用性有限。
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地形校正方法,但在地形复杂的城市场景下仍有待进一步探索。总体上立体图像在城市建筑高度估算上,精度非常可靠,但其数据获取困难、且受云量干扰严重等问题严重阻碍了其大尺度应用。
4)地基、车载和机载等激光雷达数据虽然具备了获取城市建筑物高精准三维结构信息的能力,但成本过高,不便于推广至大尺度应用。
综上可见,区域尺度的城市建筑高度反演研究,其结果数据无法映射每个建筑高度信息且缺乏高精度的建筑高度样本;然而精细尺度的城市建筑高度反演研究,则是不同数据源存在诸多不足,很难推广至区域或全球尺度计算。在未来实现“即见城市又见建筑”尺度的城市建筑高度反演,目前已有研究从GEDI自带高度数据出发,进行初步探索。目标是从密度更高、精度更高的ICESat-2/ATLAS光子点云中获取精准的城市建筑高度样本,以此支撑未来“即见城市又见建筑”尺度的城市建筑高度反演。目前,Lao et al.已基于ICESat-2/ATLAS探究了相应的城市建筑高度提取算法,但在城市光子点云去噪方面仍存在不足、以及并未考虑地形起伏对建筑高度估算的影响。
发明内容
本发明的目的在于:针对目前存在的问题,提供了一种顾及不同地形场景下的城市建筑高度提取方法,基于ICESat-2/atl03光子点云数据提出一种参数自适应的城市光子点云分类算法,可在不同地形场景下获取更精准的城市建筑高度数据。
本发明的技术方案如下:
一种顾及不同地形场景下的城市建筑高度提取方法,包括以下步骤:
从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子;
从建筑光子和地面光子中识别地面光子,从城市地面光子在水平和垂直两个方向上的分布特征,来准确识别地面光子;
识别建筑光子及建筑高度估算,在获取地面曲线后,即获取地面曲线3m以上的信号光子为建筑光子,并对建筑光子进行后处理;
结合地面曲线和建筑光子数据可估算准确的建筑高度参数。
进一步的,所述从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子具体包括以下步骤:
ICESat-2/ALT03光子点云数据读取;
基于地形自适应的高程缓冲带实现粗去噪;
光子点云等距分段,且各分段作边缘缓冲和镜像处理;
基于局部参数统计提取建筑光子及地面光子;
精确识别建筑光子和地面光子。
进一步的,所述基于局部参数统计提取建筑光子及地面光子包括以下步骤:
步骤一、设置数据处理尺度及规则化:以1km为数据处理的基本尺度;通过ACDI识别建筑光子和地面光子时,需凭借移动窗口法在相对沿轨方向作等距分段,并对每个分段进行缓冲和镜像处理;
步骤二、自适应局部参数统计及ACDI计算:在每个分段构建高程频率直方图并获取最高直方柱,并以RANSAC对相应的光子数据作线性拟合;
获取各点至拟合线的平均距离,定义其值为和/>的长度,定义τminor的长度为短轴的十倍;
将光子点云数据置于X-Y-Z三维空间内,椭球体τminor依次绕Z轴、X轴旋转变换,两次变换的坐标计算如公式如下所示:
其中,(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轴旋转β°后坐标;
获取单光子P在三维空间内的最大邻域及垂直邻域,获取垂直邻域的光子数NCvertical和最大邻域的光子数NCHorizonta并结合以下公式计算单光子P的VHDRp指标:
其中,NC(θ),...,NC(θ90°),...,NC(θ180°)是不同角度θ下的单光子的邻域点数;NCvertical是垂直邻域的光子数;NCHorizontal是最大邻域的光子数,NR是NCvertical和NCHorizontal中最大值与最小值的比值;
获取最大邻域内所有的邻域光子,并通过RANSAC对其线性拟合,并计算每个邻域光子至拟合线的距离接着即可结合以下公式计算单光子P的PHALp指标:
其中,是最大邻域中邻域光子qi和邻域中心光子p的欧式距离;distmean则是所有邻域光子至邻域中心光子的平均距离;
综合VHDRp和PHALp两个光子邻域特征值,根据以下公式计算出单光子P的ACDIp
ACDIp=PHALp+VHDRp
基于该指标从总体光子中识别出建筑光子和地面光子;
步骤三、自适应确定去噪阈值:对ACDIp阈值的训练样本范围定在[-1,2]之间,以0.05为跨度依次进行点云分类,在不同阈值i下将获取信号集和噪声集,经PCA计算后获取前三主成分的方差;接着通过以下公式计算信号集和噪声集各主成分方差的标准差和/>当/>和/>的绝对差值最大时对应的i即为ACDI最佳分割阈值,
其中,Var(PCAi)是第i主成分的方差;Varmean是主成分分析的前三主成分方差的平均值;Std是前三主成分方差的标准差;和/>是当ACDI分类阈值为i时,被分类为信号光子和噪声光子的数据,通过主成分分析后前三主成分方差的标准差;另外,当和/>差值达到最大时,以此时的i为ACDI的最佳分割阈值。
进一步的,所述粗去噪包括以下步骤:
以1Km为ICESat-2剖面点云数据处理尺度,根据50m跨度的移动窗口在沿轨方向建立等距分段,在各分段内,在高程范围内划分10个等距分箱并作高程频率统计;获取频数最高的分箱,并以在其上下方各取一分箱的位置为缓冲区,除此之外的分箱将视作噪声被剔除。
进一步的,所述从建筑光子和地面光子中识别地面光子具体包括以下步骤:
移动窗口下通过百分位统计法识别初始地面光子;
基于地面光子的水平分布特征对地面光子初过滤;
区域生长识别所有地面光子,同时以MDNN突变时作为停止区域生长标志;
基于地面光子的垂直分布特征对地面光子终过滤,得到精确的地面光子。
进一步的,所述从建筑光子和地面光子中识别地面光子具体包括以下步骤:
初始地面光子提取和去噪:沿初始地面光子水平分布方向进行迭代RANSAC线性拟合,并以拟合点至拟合线平均距离、拟合线斜率阈值作为分割条件,将初始地面光子聚类为若干线状点簇;过滤线状点簇,包括全局过滤,即根据线状点簇的拟合斜率进行过滤;以及局部过滤;将每个线性点簇当作一个泵,泵可将水流从地势高的地方输送至地势低的地方,也可在一定增压能力下将水流从地势低的地方输送至地势高的地方,泵的增压能力通过地势高程差来衡量;那么可在一定高程差阈值下去除伪地面光子;以连续分布的三个线性点簇为步长,获取中间点簇相对两侧点簇的高度差,相对两侧起伏较大的点簇被滤除,该过程迭代进行直至相邻线性点簇的相对位置趋于平稳;
区域生长完善地面光子识别:对去噪后的初始地面光子构建KD-Tree索引,依次获取单光子的最大邻域方向,结合椭圆结构对单光子进行区域生长,以完善地面光子识别;通过以下公式计算MDNN,并设置MDNN阈值,在区域生长中若MDNN出现突变或零变化则停止区域生长,
其中,N表示椭圆生长域内光子数据;表示第i组最邻近点对pi和qi间的欧式距离;
地面光子二次去噪:需要对已识别的地面光子进行分段,在适宜窗口下进行该去噪处理;在每个窗口内,根据识别的地面光子构建高程频率直方图,将所有的直方柱聚类分组,实现地面光子的点簇分类;定义每个地面点簇为待判断的目标点簇,并根据目标点簇的横坐标范围获取原始数据中的对应光子,构建该部分光子数据的高程频率直方图并对直方柱进行聚类分组,将定位目标点簇在分组中的位置;根据目标点簇在原始数据高程频率直方图中的相对位置来判断目标点簇是否为地面光子;
地面曲线获取:由每个窗口中获取最终的地面光子数据,并联合傅里叶变换和Lowess拟合地面曲线,结合B样条插值算法将地面曲线插值至原始光子分辨率。
进一步的,所述结合地面曲线和建筑光子数据可估算准确的建筑高度参数具体包括以下步骤:
联合Lowess拟合和B样条插值获取信号光子分辨率的地面曲线,得到精准的地面曲线;
获取地面曲线3m以上的信号光子为建筑光子,对建筑光子后处理:以区域生长完善建筑光子,以KNN剔除离群伪建筑光子,得到精准的建筑光子;
结合精准的地面曲线和建筑光子数据可估算准确的建筑高度参数。
进一步的,所述对建筑光子进行后处理包括以下步骤:
获取建筑光子和地面光子集合、地面曲线,定义地面曲线3m以上的光子为初始建筑光子;
以MDNN作为终止区域生长的区域生长法来完善建筑光子;
通过KNN算法来检测离群光子,并采用OTSU算法确定阈值进行剔除,KNN算法的最邻近点数为5个。
进一步的,所述估算准确的建筑高度参数包括以下步骤:
获取建筑光子和地面曲线后,利用B样条插值法将地面曲线插值至建筑光子分辨率;将每个建筑光子将匹配对应位置上的地面光子,计算两者的高程差为建筑高度;根据实景三维模型量算的高精度建筑高度数据为参照,然后分别计算均方根误差、平均绝对误差和决定系数三个指标来评估不同场景下的建筑高度估算精度。
与现有的技术相比本发明的有益效果是:
1、一种顾及不同地形场景下的城市建筑高度提取方法,所提出方法框架可准确识别不同城市地形场景下的建筑光子和地面光子。通过量化不同场景下建筑高度的估算精度,与参考建筑高度具备较高一致性;平坦地形下,RMSE在0.230m-0.315m之间,MAE在0.161m-0.243m之间;起伏地形下RMSE在0.912m-1.424m之间,MAE在0.647m-0.871m之间。与同行研究相比,该方法在不同场景下表现出更良好的适应性;
2、一种顾及不同地形场景下的城市建筑高度提取方法,提出一种基于单光子局部特征ACDI的自适应信号光子提取方法,并且对不同的城市场景作了验证(包括日间-起伏地形、日间-平坦地形、夜间-起伏地形、夜间-平坦地形),相比Lao et al.的研究,所提出的算法能够更准确的捕捉ICESat-2/atlas沿轨方向的信号光子(即建筑光子和地面光子)数据。其中,ACDI是在自适应椭球体结构下,获取单光子邻域参数PHAL和VHDR来计算,最后在PCA辅助下自适应获取分割信号光子和其他光子(即噪声光子和其他地物光子)的阈值,与OTSU算法对比发现,该方法的阈值分割算法效果更佳;
3、一种顾及不同地形场景下的城市建筑高度提取方法,针对起伏地形下,城市地形识别不稳定的问题。根据地面光子的水平分布特征和垂直分布特征提出一种简洁、高效的地面识别方法。相比同行Lao et al.的研究,该方法能更稳定的解决城市地面识别问题。其方法划分三个过程:a.通过移动窗口获取局部高程最低点为初始地面光子;b.借鉴Li etal.研究中根据水泵布设准则剔除伪地面光子的思路,对初始地面光子作了粗去噪;c.将未剔除干净的伪地面光子(源自建筑光子)通过区域生长从“离群点”延伸至“带状建筑光子”,接着基于高程频率直方图,从光子垂直分布方向划分出地面光子簇和伪地面光子簇(即“带状建筑光子”)。结果表明,该方法显然是通用于不同地形场景下的地面光子识别;
4、一种顾及不同地形场景下的城市建筑高度提取方法,在获取准确的地面光子后,联合傅里叶变换和Lowess算法拟合了平稳的地形曲线,接着提取地面曲线3m以上的信号光子为建筑光子。同时,为减少建筑光子的缺失,以区域生长和KNN方法分别对建筑光子进行补齐及剔除离群噪声光子。
附图说明
图1为一种顾及不同地形场景下的城市建筑高度提取方法的流程图。
图2为一种顾及不同地形场景下的城市建筑高度提取方法的不同城市场景下所估算建筑高度的精度验证图一。
图3为一种顾及不同地形场景下的城市建筑高度提取方法的不同城市场景下所估算建筑高度的精度验证图二。
具体实施方式
需要说明的是,术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来,而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下,由语句“包括一个……”限定的要素,并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
下面结合实施例对本发明的特征和性能作进一步的详细描述。
请参阅图1-3,一种顾及不同地形场景下的城市建筑高度提取方法,包括以下步骤:
从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子;
从建筑光子和地面光子中识别地面光子,从城市地面光子在水平和垂直两个方向上的分布特征,来准确识别地面光子;
识别建筑光子及建筑高度估算,在获取地面曲线后,即获取地面曲线3m以上的信号光子为建筑光子,并对建筑光子进行后处理;
结合地面曲线和建筑光子数据可估算准确的建筑高度参数。
从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子具体包括以下步骤:
ICESat-2/ALT03光子点云数据读取;
基于地形自适应的高程缓冲带实现粗去噪;
光子点云等距分段,且各分段作边缘缓冲和镜像处理;
基于局部参数统计提取建筑光子及地面光子;
精确识别建筑光子和地面光子。
从建筑光子和地面光子中识别地面光子具体包括以下步骤:
移动窗口下通过百分位统计法识别初始地面光子;
基于地面光子的水平分布特征对地面光子初过滤;
区域生长识别所有地面光子,同时以MDNN突变时作为停止区域生长标志;
基于地面光子的垂直分布特征对地面光子终过滤,得到精确的地面光子。
结合地面曲线和建筑光子数据可估算准确的建筑高度参数具体包括以下步骤:
联合Lowess拟合和B样条插值获取信号光子分辨率的地面曲线,得到精准的地面曲线;
获取地面曲线2.5m以上的信号光子为建筑光子,对建筑光子后处理:以区域生长完善建筑光子,以KNN剔除离群伪建筑光子,得到精准的建筑光子;
结合精准的地面曲线和建筑光子数据可估算准确的建筑高度参数。
粗去噪包括以下步骤:
以1Km为ICESat-2剖面点云数据处理尺度,根据50m跨度的移动窗口在沿轨方向建立等距分段,在各分段内,在高程范围内划分10个等距分箱并作高程频率统计;获取频数最高的分箱,并以在其上下方各取一分箱的位置为缓冲区,除此之外的分箱将视作噪声被剔除。
基于局部参数统计提取建筑光子及地面光子包括以下步骤:
步骤一、设置数据处理尺度及规则化:以1km为数据处理的基本尺度;通过ACDI识别建筑光子和地面光子时,需凭借移动窗口法在相对沿轨方向作等距分段,并对每个分段进行缓冲和镜像处理;
步骤二、自适应局部参数统计及ACDI计算:在每个分段构建高程频率直方图并获取最高直方柱,并以RANSAC对相应的光子数据作线性拟合;
获取各点至拟合线的平均距离,定义其值为和/>的长度,定义τminor的长度为短轴的十倍;
将光子点云数据置于X-Y-Z三维空间内,椭球体τminor依次绕Z轴、X轴旋转变换,两次变换的坐标计算如公式如下所示:
其中,(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轴旋转β°后坐标;
获取单光子P在三维空间内的最大邻域及垂直邻域,获取垂直邻域的光子数NCvertical和最大邻域的光子数NCHorizonta并结合以下公式计算单光子P的VHDRp指标:
其中,NC(θ),...,NC(θ90°),...,NC(θ180°)是不同角度θ下的单光子的邻域点数;NCvertical是垂直邻域的光子数;NCHorizontal是最大邻域的光子数,NR是NCvertical和NCHorizontal中最大值与最小值的比值;
获取最大邻域内所有的邻域光子,并通过RANSAC对其线性拟合,并计算每个邻域光子至拟合线的距离接着即可结合以下公式计算单光子P的PHALp指标:
其中,是最大邻域中邻域光子qi和邻域中心光子p的欧式距离;distmean则是所有邻域光子至邻域中心光子的平均距离;
综合VHDRp和PHALp两个光子邻域特征值,根据以下公式计算出单光子P的ACDIp
ACDIp=PHALp+VHDRp
基于该指标从总体光子中识别出建筑光子和地面光子;
步骤三、自适应确定去噪阈值:对ACDIp阈值的训练样本范围定在[-1,2]之间,以0.05为跨度依次进行点云分类,在不同阈值i下将获取信号集和噪声集,经PCA计算后获取前三主成分的方差;接着通过以下公式计算信号集和噪声集各主成分方差的标准差和/>当/>和/>的绝对差值最大时对应的i即为ACDI最佳分割阈值,
其中,Var(PCAi)是第i主成分的方差;Varmean是主成分分析的前三主成分方差的平均值;Std是前三主成分方差的标准差;和/>是当ACDI分类阈值为i时,被分类为信号光子和噪声光子的数据,通过主成分分析后前三主成分方差的标准差;另外,当和/>差值达到最大时,以此时的i为ACDI的最佳分割阈值。
从建筑光子和地面光子中识别地面光子具体包括以下步骤:
初始地面光子提取和去噪:沿初始地面光子水平分布方向进行迭代RANSAC线性拟合,并以拟合点至拟合线平均距离、拟合线斜率阈值作为分割条件,将初始地面光子聚类为若干线状点簇;过滤线状点簇,包括全局过滤,即根据线状点簇的拟合斜率进行过滤;以及局部过滤;将每个线性点簇当作一个泵,泵可将水流从地势高的地方输送至地势低的地方,也可在一定增压能力下将水流从地势低的地方输送至地势高的地方,泵的增压能力通过地势高程差来衡量;那么可在一定高程差阈值下去除伪地面光子;以连续分布的三个线性点簇为步长,获取中间点簇相对两侧点簇的高度差,相对两侧起伏较大的点簇被滤除,该过程迭代进行直至相邻线性点簇的相对位置趋于平稳;
区域生长完善地面光子识别:对去噪后的初始地面光子构建KD-Tree索引,依次获取单光子的最大邻域方向,结合椭圆结构对单光子进行区域生长,以完善地面光子识别;通过以下公式计算MDNN,并设置MDNN阈值,在区域生长中若MDNN出现突变或零变化则停止区域生长,
其中,N表示椭圆生长域内光子数据;表示第i组最邻近点对pi和qi间的欧式距离;
地面光子二次去噪:需要对已识别的地面光子进行分段,在适宜窗口下进行该去噪处理;在每个窗口内,根据识别的地面光子构建高程频率直方图,将所有的直方柱聚类分组,实现地面光子的点簇分类;定义每个地面点簇为待判断的目标点簇,并根据目标点簇的横坐标范围获取原始数据中的对应光子,构建该部分光子数据的高程频率直方图并对直方柱进行聚类分组,将定位目标点簇在分组中的位置;根据目标点簇在原始数据高程频率直方图中的相对位置来判断目标点簇是否为地面光子;该步骤提高识别伪地面光子的可能性。
地面曲线获取:由每个窗口中获取最终的地面光子数据,并联合傅里叶变换和Lowess拟合地面曲线,结合B样条插值算法将地面曲线插值至原始数据分辨率。能够针对不同的地形场景实现地面光子提取。
对建筑光子进行后处理包括以下步骤:
获取建筑光子和地面光子集合、地面曲线,定义地面曲线2.5m以上的光子为初始建筑光子;
以MDNN作为终止区域生长的区域生长法来完善建筑光子;
通过KNN算法来检测离群光子,并采用OTSU算法确定阈值进行剔除,KNN算法的最邻近点数为5个。
估算准确的建筑高度参数包括以下步骤:
获取建筑光子和地面曲线后,利用B样条插值法将地面曲线插值至建筑光子分辨率;将每个建筑光子将匹配对应位置上的地面光子,计算两者的高程差为建筑高度;以PR3DM量算的高精度建筑高度数据为参照,然后分别计算均方根误差、平均绝对误差和决定系数三个指标来评估不同场景下的建筑高度估算精度。
实验验证
ICESat-2/ATLAS(Advanced Topographic Laser Altimeter System)采用微脉冲光子计数激光雷达技术,其沿地面轨道每隔0.7m发射波长532nm、频率10kHz的激光脉冲,激光脉冲经分光处理后为3对6束。光束对在横轨方向上的间距为3.3km,每对光束包括一个强光束和一个弱光束(能量比4:1),两束光束的距离为90m。激光脚印的直径约为17m。相比ICESat和GEDI,可沿轨获取高密度、高频率和高精度的光子点云数据,目前已证实在反演植被冠层高度和水下地形等方面具备强大优势,则试图基于ICESat-2光子点云数据开发自适应较高的城市建筑高度提取方法。ICESat-2产品分为三个层次,使用的是二级产品(ATL03)全球定位光子数据。另外,ATL03数据可从国家冰雪数据中心(https://nsidc.org/data/atl03/)以HDF5格式免费下载,其光子数据的坐标参考为1984年世界大地测量系统(WGS84)。所选用的研究数据信息可见表1所示。
表1ICESat-2/ATLAS试验数据信息
为提升昆明的城市管理精细化、动态化和信息化水平,为政府决策提供科学依据,以实现智慧化城市规划建设管理目标。昆明市政府决定启动昆明市城市规划建设管理数字化应用项目。项目核心任务之一便是采用倾斜航空摄影等前沿测绘技术,对昆明市主城建成区430平方公里范围的建筑物、城市用地等进行大范围、可视化动态监测,并生产高精度实景三维模型(Photorealistic 3D model,PR3DM)。本专利的参考城市建筑高度数据来源于该PR3DM。关于倾斜摄影数据采集过程,采用全球最先进的、配置最高的Leica RCD30Oblique Penta(8000万Pixels)的量测型高精度倾斜航空相机并搭载在运-12轻型运输机上展开。考虑到昆明市建筑物较为密集等情况,设置地面分辨率时要求下视镜头不低于0.07m、以及下视镜头航向重叠率不低于80%、旁向重叠率不低于70%,确保城市三维模型质量。基于采集的倾斜摄影数据,依次进行几何校正、多视匹配、三角网(TIN)构建、自动赋予纹理和人工修正模型局部等,输出了昆明主城区高精度的PR3DM。对于PR3DM的核心信息见表1所示。另外,需特别说明的是,ICESat-2/ATL03光子点云的坐标参考是WGS-84,而PR3DM的坐标参考是CGCS2000,虽在不同IRTF框架、历元下,但两者的坐标偏差仅在分米量级,并不影响本专利的城市建筑高度精度验证。
表2城市实景三维模型产品生产规格信息
选自中国云南省昆明市主城区,总面积21473km2,市中心海拔约1891m,总体地势北部高南部低,由北向南呈阶梯状逐渐降低。区域内涵盖了城市中的住宅、平房、商场和工棚等,具有建筑类型多样、建筑分布复杂等特点。而本专利拟验证“日间-起伏地形、日间-平坦地形、夜间-起伏地形、夜间-平坦地形”四个场景下的城市建筑高度估算精度,无论是ICESat-2/atl03数据覆盖情况,或是地形起伏情况,昆明市主城区均能满足本专利方法的测试需求。
根据上述算法所获取的准确地面曲线和准确建筑光子,在插值加密后的地面曲线上提取每个建筑光子对应位置的地面光子,通过计算高程差来获取沿轨方向所分布建筑的高度。图2中描绘了文中若干验证点上由ICESat-2/atlas估算的建筑高度与由PR3DM获取建筑高度之间的回归验证图,其中从“日间-起伏地形、日间-平坦地形、夜间-起伏地形、夜间-平坦地形”四个场景分别验证。由图2线性回归和相关参数可知,ICESat-2/atlas估算的建筑高度和由PR3DM获取建筑高度呈显著线性关系,其中平坦地形下相关性较高,R2逼近于1,RMSE分别为0.315m和0.230m,MAE分别为0.243m和0.161m。相比之下起伏地形下则稍差,R2同样也逼近于1,但RMSE为0.912m和1.424m,MAE为0.647m和0.871m。Lao et al,方法主要针对平坦地形下的城市建筑高度提取,其结果显示RMSE在0.35-0.45m之间,MAE在0.28-0.35m之间,对比结果来看本专利所估算的建筑高度精度有进一步提高。结合上述信号光子识别和地面光子提取结果来看,在起伏地形下其效果很差;而准确的建筑光子和地面曲线是准确估算建筑高度的先决条件,由此可见Lao et al.方法在应用于起伏地形下的建筑高度提取时仍需较大改进。相比之下,该方法在起伏地形下的日间、夜间数据中仍能获取较高精度,R2高达99.196%、99.039%,RMSE也仅为0.912m和1.424m,MAE为0.647m和0.871m。综上可见,该方法在不同城市场景下均可获取较高精度的城市建筑高度数据。
另外,还通过误差棒图来分析城市建筑高度的估算精度(图3),图内柱形表示建筑高度估算值与真实值之间的误差均值,误差棒长度则是由置信度区间来表示误差均值的不确定性程度。根据图内柱形的长度可知四个场景下所估测建筑高度的误差均值:“夜间-起伏地形”>“日间-平坦地形”=“日间-起伏地形”>“日间-起伏地形”,根据误差棒长度来看,起伏地形下建筑高度的误差不确定性明显高于平坦地形。根据柱形的正负向分布、上下误差棒的相对长度来看,日间数据下的建筑高度数据存在被低估、夜间数据下的建筑高度数据存在被高估。无论夜间数据还是日间数据,建筑光子远离地表且其局部密度大,因此易与周边噪声分离;但地面光子不同,地表附近的低矮地物信息对其干扰较大,因此城市建筑高度误差多源自地面识别误差。那么,相对于日间数据建筑高度被低估情况,实际是地面曲线被高估;而夜间数据建筑高度被高估情况,实际是地面曲线被低估。可见基于ICESat-2/atlas光子点云估算城市建筑高度时,如何更精准的识别地面曲线仍是一个需要持续探究的过程,尤其是起伏地形下的城市场景。
以上所述实施例仅表达了本申请的具体实施方式,其描述较为具体和详细,但并不能因此而理解为对本申请保护范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本申请技术方案构思的前提下,还可以做出若干变形和改进,这些都属于本申请的保护范围。

Claims (9)

1.一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,包括以下步骤:
从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子;
从建筑光子和地面光子中识别地面光子,从城市地面光子在水平和垂直两个方向上的分布特征,来准确识别地面光子;
识别建筑光子及建筑高度估算,在获取地面曲线后,即获取地面曲线3m以上的信号光子为建筑光子,并对建筑光子进行后处理;
结合地面曲线和建筑光子数据可估算准确的建筑高度参数。
2.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述从原始光子中识别建筑光子和地面光子,基于局部参数统计原理提取光子局部特征ACDI来自适应的识别建筑光子及地面光子具体包括以下步骤:
ICESat-2/ALT03光子点云数据读取;
基于地形自适应的高程缓冲带实现粗去噪;
光子点云等距分段,且各分段作边缘缓冲和镜像处理;
基于局部参数统计提取建筑光子及地面光子;
精确识别建筑光子和地面光子。
3.根据权利要求2所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述粗去噪包括以下步骤:
以1Km为ICESat-2剖面点云数据处理尺度,根据50m跨度的移动窗口在沿轨方向建立等距分段,在各分段内,在高程范围内划分10个等距分箱并作高程频率统计;获取频数最高的分箱,并以在其上下方各取一分箱的位置为缓冲区,除此之外的分箱将视作噪声被剔除。
4.根据权利要求1或2所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述基于局部参数统计提取建筑光子及地面光子包括以下步骤:
步骤一、设置数据处理尺度及规则化:以1km为数据处理的基本尺度;通过ACDI识别建筑光子和地面光子时,需凭借移动窗口法在相对沿轨方向作等距分段,并对每个分段进行缓冲和镜像处理;
步骤二、自适应局部参数统计及ACDI计算:在每个分段构建高程频率直方图并获取最高直方柱,并以RANSAC对相应的光子数据作线性拟合;
获取各点至拟合线的平均距离,定义其值为和/>的长度,定义τminor的长度为短轴的十倍;
将光子点云数据置于X-Y-Z三维空间内,椭球体τminor依次绕Z轴、X轴旋转变换,两次变换的坐标计算如公式如下所示:
其中,(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轴旋转β°后坐标;
获取单光子P在三维空间内的最大邻域及垂直邻域,获取垂直邻域的光子数NCvertical和最大邻域的光子数NCHorizonta并结合以下公式计算单光子P的VHDRp指标:
其中,NC(θ),…,NC(θ90°),…,NC(θ180°)是不同角度θ下的单光子的邻域点数;NCvertical是垂直邻域的光子数;NCHorizontal是最大邻域的光子数,NR是NCvertical和NCHorizontal中最大值与最小值的比值;
获取最大邻域内所有的邻域光子,并通过RANSAC对其线性拟合,并计算每个邻域光子至拟合线的距离接着即可结合以下公式计算单光子P的PHALp指标:
其中,是最大邻域中邻域光子qi和邻域中心光子p的欧式距离;distmean则是所有邻域光子至邻域中心光子的平均距离;
综合VHDRp和PHALp两个光子邻域特征值,根据以下公式计算出单光子P的ACDIp指标:
ACDIp=PHALp+VHDRp
基于该指标从总体光子中识别出建筑光子和地面光子;
步骤三、自适应确定去噪阈值:对ACDIp阈值的训练样本范围定在[-1,2]之间,以0.05为跨度依次进行点云分类,在不同阈值i下将获取信号集和噪声集,经PCA计算后获取前三主成分的方差;接着通过以下公式计算信号集和噪声集各主成分方差的标准差当/>和/>的绝对差值最大时对应的i即为ACDI最佳分割阈值,
其中,Var(PCAi)是第i主成分的方差;Varmean是主成分分析的前三主成分方差的平均值;Std是前三主成分方差的标准差;和/>是当ACDI分类阈值为i时,被分类为信号光子和噪声光子的数据,通过主成分分析后前三主成分方差的标准差;另外,当和/>差值达到最大时,以此时的i为ACDI的最佳分割阈值。
5.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述从建筑光子和地面光子中识别地面光子具体包括以下步骤:
移动窗口下通过百分位统计法识别初始地面光子;
基于地面光子的水平分布特征对地面光子初过滤;
区域生长识别所有地面光子,同时以MDNN突变时作为停止区域生长标志;
基于地面光子的垂直分布特征对地面光子终过滤,得到精确的地面光子。
6.根据权利要求1或5所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述从建筑光子和地面光子中识别地面光子具体包括以下步骤:
初始地面光子提取和去噪:沿初始地面光子水平分布方向进行迭代RANSAC线性拟合,并以拟合点至拟合线平均距离、拟合线斜率阈值作为分割条件,将初始地面光子聚类为若干线状点簇;过滤线状点簇,包括全局过滤,即根据线状点簇的拟合斜率进行过滤;以及局部过滤;将每个线性点簇当作一个泵;在一定高程差阈值下去除伪地面光子;以连续分布的三个线性点簇为步长,获取中间点簇相对两侧点簇的高度差,相对两侧起伏较大的点簇被滤除,该过程迭代进行直至相邻线性点簇的相对位置趋于平稳;
区域生长完善地面光子识别:对去噪后的初始地面光子构建KD-Tree索引,依次获取单光子的最大邻域方向,结合椭圆结构对单光子进行区域生长,以完善地面光子识别;通过以下公式计算MDNN,并设置MDNN阈值,在区域生长中若MDNN出现突变或零变化则停止区域生长,
其中,N表示椭圆生长域内光子数据;表示第i组最邻近点对pi和qi间的欧式距离;
地面光子二次去噪:需要对已识别的地面光子进行分段,在适宜窗口下进行该去噪处理;在每个窗口内,根据识别的地面光子构建高程频率直方图,将所有的直方柱聚类分组,实现地面光子的点簇分类;定义每个地面点簇为待判断的目标点簇,并根据目标点簇的横坐标范围获取原始数据中的对应光子,构建该部分光子数据的高程频率直方图并对直方柱进行聚类分组,将定位目标点簇在分组中的位置;根据目标点簇在原始数据高程频率直方图中的相对位置来判断目标点簇是否为地面光子;
地面曲线获取:由每个窗口中获取最终的地面光子数据,并联合傅里叶变换和Lowess拟合地面曲线,并结合B样条插值算法将地面曲线插值至原始光子分辨率。
7.根据权利要求1所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述结合地面曲线和建筑光子数据可估算准确的建筑高度参数具体包括以下步骤:
联合Lowess拟合和B样条插值获取信号光子分辨率的地面曲线,得到精准的地面曲线;
获取地面曲线3m以上的信号光子为建筑光子,对建筑光子后处理:以区域生长完善建筑光子,以KNN剔除离群伪建筑光子,得到精准的建筑光子;
结合精准的地面曲线和建筑光子数据可估算准确的建筑高度参数。
8.根据权利要求1或7所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述对建筑光子进行后处理包括以下步骤:
获取建筑光子和地面光子集合、地面曲线,定义地面曲线3m以上的光子为初始建筑光子;
以MDNN作为终止区域生长的区域生长法来完善建筑光子;
通过KNN算法来检测离群光子,并采用OTSU算法确定阈值进行剔除,KNN算法的最邻近点数为5个。
9.根据权利要求1或7所述的一种顾及不同地形场景下的城市建筑高度提取方法,其特征在于,所述估算准确的建筑高度参数包括以下步骤:
获取建筑光子和地面曲线后,利用B样条插值法将地面曲线插值至建筑光子分辨率;将每个建筑光子匹配对应位置上的地面光子,计算两者的高程差为建筑高度;根据实景三维模型量算的高精度建筑高度数据为参照,然后分别计算均方根误差、平均绝对误差和决定系数三个指标来评估不同场景下的建筑高度估算精度。
CN202311775357.2A 2023-12-21 2023-12-21 一种顾及不同地形场景下的城市建筑高度提取方法 Pending CN117932333A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311775357.2A CN117932333A (zh) 2023-12-21 2023-12-21 一种顾及不同地形场景下的城市建筑高度提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311775357.2A CN117932333A (zh) 2023-12-21 2023-12-21 一种顾及不同地形场景下的城市建筑高度提取方法

Publications (1)

Publication Number Publication Date
CN117932333A true CN117932333A (zh) 2024-04-26

Family

ID=90769223

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311775357.2A Pending CN117932333A (zh) 2023-12-21 2023-12-21 一种顾及不同地形场景下的城市建筑高度提取方法

Country Status (1)

Country Link
CN (1) CN117932333A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118519131A (zh) * 2024-07-18 2024-08-20 武汉大学 顾及地形影响双重性的gedi激光雷达冠层高度数据校正方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118519131A (zh) * 2024-07-18 2024-08-20 武汉大学 顾及地形影响双重性的gedi激光雷达冠层高度数据校正方法

Similar Documents

Publication Publication Date Title
CN108510467B (zh) 基于深度可变形卷积神经网络的sar图像目标识别方法
CN111079611B (zh) 一种道路面及其标志线的自动提取方法
CN108197583B (zh) 基于图割优化和影像结构特征的建筑物变化检测方法
CN106204547B (zh) 从车载激光扫描点云中自动提取杆状地物空间位置的方法
CN108710875A (zh) 一种基于深度学习的航拍公路车辆计数方法及装置
CN104361590B (zh) 一种控制点自适应分布的高分辨率遥感影像配准方法
CN106204705A (zh) 一种基于多线激光雷达的3d点云分割方法
CN106970375A (zh) 一种机载激光雷达点云中自动提取建筑物信息的方法
WO2023060632A1 (zh) 基于点云数据的街景地物多维度提取方法和系统
CN110927765B (zh) 激光雷达与卫星导航融合的目标在线定位方法
CN112669333B (zh) 一种单木信息提取方法
CN111007531A (zh) 一种基于激光点云数据的道路边沿检测方法
CN109711256B (zh) 一种低空复杂背景无人机目标检测方法
CN108562885B (zh) 一种高压输电线路机载LiDAR点云提取方法
CN112085778B (zh) 基于超像素和形态学的倾斜摄影违法建筑检测方法及系统
CN118314157A (zh) 基于人工智能的遥感影像智能分割方法
CN117932333A (zh) 一种顾及不同地形场景下的城市建筑高度提取方法
Yao et al. Automatic extraction of road markings from mobile laser-point cloud using intensity data
CN102013014B (zh) 一种高分辨率遥感图像多类目标特征模型的建立方法
CN116579949A (zh) 适用于城市多噪声环境下机载点云地面点滤波方法
CN106709432B (zh) 基于双目立体视觉的人头检测计数方法
CN118918482B (zh) 基于遥感图像的自然资源测量方法及系统
Dong et al. Pixel-level intelligent segmentation and measurement method for pavement multiple damages based on mobile deep learning
Huang et al. Urban building height extraction accommodating various terrain scenes using ICESat-2/ATLAS data
CN115236643B (zh) 一种传感器标定方法、系统、装置、电子设备及介质

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