CN111047695B - 一种城市群高度空间信息及等高线的提取方法 - Google Patents

一种城市群高度空间信息及等高线的提取方法 Download PDF

Info

Publication number
CN111047695B
CN111047695B CN201911221885.7A CN201911221885A CN111047695B CN 111047695 B CN111047695 B CN 111047695B CN 201911221885 A CN201911221885 A CN 201911221885A CN 111047695 B CN111047695 B CN 111047695B
Authority
CN
China
Prior art keywords
urban
image
building
height
point
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.)
Active
Application number
CN201911221885.7A
Other languages
English (en)
Other versions
CN111047695A (zh
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.)
Institute of Geographic Sciences and Natural Resources of CAS
Original Assignee
Institute of Geographic Sciences and Natural Resources of CAS
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 Institute of Geographic Sciences and Natural Resources of CAS filed Critical Institute of Geographic Sciences and Natural Resources of CAS
Priority to CN201911221885.7A priority Critical patent/CN111047695B/zh
Publication of CN111047695A publication Critical patent/CN111047695A/zh
Application granted granted Critical
Publication of CN111047695B publication Critical patent/CN111047695B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Graphics (AREA)
  • Multimedia (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Remote Sensing (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Image Processing (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明涉及一种城市群高度空间信息及等高线的提取方法，针对城市建筑物的提取融合了遥感数据和社会感知数据，建立了地块尺度的高/低密度建筑区，结合多视角差值特征图，实现了城市建筑物的提取；针对城市建筑物高度的提取，利用资源三号前视、后视和正视影像两两匹配生成密集点云，采用TIN线性插值法求取数字表面模型，引入点云分类的思想，获取精确的地面点，建立数字高程模型。本发明充分利用资源三号前视、后视和正视影像，创新性的构建了提取城市建筑物信息及建筑物高度信息的系统框架，并以此为基础获得城市等高线、城市起伏度和城市生长高度的模型，揭示了城市横向面积扩张和纵向高度增长的演变格局。

Description

一种城市群高度空间信息及等高线的提取方法
技术领域
本发明涉及数据信息处理领域，特别是涉及一种信息提出方法，尤其涉及一种城市群高度空间信息及等高线的提取方法。
背景技术
随着中国城市人口的不断增加以及城市化进程的不断加快，城市用地在急剧扩张。城市扩张包括城市建筑物横向的水平扩张以及纵向的高度增长，城市建筑物是城市规划、管理及信息化的重要组成部分，建筑物高度是进行城市监测、城市三维重建、城市污染扩散模拟等过程的基础信息。因此，综合考虑城市横向面积扩张和纵向高度增长对于城市空间形态特征的研究和城市景观的模拟至关重要。
随着高分辨率遥感影像数据的广泛应用，实时的、大面积的提取城市建筑物信息和监测城市横向面积扩张的研究也越来越多地受到重视。然而，大面积的提取城市建筑物高度和监测城市高度增长依然面临着一些挑战。现有的提取建筑物高度的技术主要是利用建筑物阴影进行建筑物高度的提取，但是阴影测高在实际应用中面临着诸多不足。首先，受周围地物的干扰，从图像中提取的建筑物阴影区域可能会存在缺失或变形。其次，在建筑物密集的建成区，建筑物的阴影存在很多遮挡和重叠的问题。最后，阴影区域的形状多不规则，边缘参差不齐，也给建筑物阴影长度的准确量算带来困扰。
资源三号测绘卫星(简称ZY-3)是我国第一颗民用高分辨率光学传输型测绘卫星，它搭载了4台光学相机，包括一台地面分辨率2.1m的正视全色TDI CCD相机、两台地面分辨率3.6m的前视和后视全色TDI CCD相机、一台地面分辨率5.8m的正视多光谱相机。它获取的数据可用于地形图制图、高程建模以及资源调查等。资源三号卫星是民用三线阵立体测图卫星，可以满足1:50000测图精度需要，达到甚至超过国际同类分辨率卫星(SPOT5、ALOS)的几何质量。资源三号卫星上的正、前、后视相机能够便捷地获取到相同地区的立体观测像对，提供了大量的三维信息，具有重要的现实意义。资源三号卫星数据的主要特点有：(1)卫星重访周期为5天，具备立体测绘与资源调查两种模式；(2)定位精度高，影像的控制定位精度优于1个像素；立体像对幅宽为52公里，基高比为0.85-0.95，能够满足1：50000比例尺立体测图需求；正视影像空间分辨率为2.1米,可满足1:25000比例尺地形图更新需求。(3)影像信息量丰富，资源三号卫星的数据量化值为10位，有利于影像的目视解译、分类识别和特征匹配精度的改善。
霍少峰,顾行发et al.(2017)以北京市城区为例，采用资源三号高分辨率卫星前视影像，通过对建筑物侧面及阴影的提取，结合太阳的高度角和方位角，利用建筑物、太阳及卫星的空间几何关系，计算出建筑物侧面与阴影在太阳方位角方向上的长度比例关系，依此比例关系估算出其他建筑物的高度信息。赵志明,周小成et al.(2015)利用资源三号遥感影像构建了形态学建筑物指数MBI和阴影指数MSI，结合多尺度优化度优化分割和面向对象分类技术提取了影像建筑物及其阴影轮廓，通过建立建筑物高度与其阴影成像的几何关系模型，结合太阳和卫星的方位角、高度角反演出了建筑物高度信息。
从目前的研究来看，在数据上，资源三号卫星包含丰富的前视、后视、正视信息，是提取建筑物高度信息的优势数据，但应用尚少。在方法上，目前的研究大多构建各种指数，基于几何关系模型、太阳与卫星的角度来反演建筑物的高度信息。但是，对于提取城市建筑物信息及建筑物高度信息的系统框架几乎没有涉及，不利于研究城市横向面积扩张和纵向高度增长的演变格局。
发明内容
本发明的目的是提出一种城市群高度空间信息及等高线的提取方法，并以此为基础获得城市等高线、城市起伏度和城市生长高度的模型，以揭示城市横向面积扩张和纵向高度增长的演变格局。
为实现上述目的，本发明提供了一种城市群高度空间信息的提取方法，所述提取方法包括如下步骤：
步骤S101，将资源三号卫星多光谱影像用SVM算法提取出植被分布图层，再利用OpenStreetMap路网数据经过预处理得到地块分割图层；
步骤S102，将所述植被分布图层和所述地块分割图层进行叠加分析，得到每个地块的植被占比，并通过建立阈值判别体系得到高密集建筑区和低密度建筑区；
步骤S103，将资源三号卫星三线阵卫星影像用多视角差值法得到多视角差值特征图，并对所述多视角差值特征图进行去除阴影处理；
步骤S104，将高密集建筑区和低密度建筑区与去除阴影后的多视角差值特征图进行叠加分析，并通过建立阈值判别体系提取出城市建筑物信息；
步骤S105，将资源三号卫星三线阵卫星影像用计算连接点和区域网平差方法生成点云数据；
步骤S106，采用TIN线性插值法进行点云插值得到数字表面模型DSM，再对点云数据进行分类得到数字高程模型DEM，之后求取数字表面模型DSM与数字高程模型DEM的差值，得到nDSM；
步骤S107，通过步骤S104中提取出的建筑物信息获得城市建筑物图层，叠加城市建筑物图层和nDSM，采用掩膜处理提取城市建筑物高度信息。
步骤S108，根据城市建筑物信息及建筑物高度信息计算城市生长高度与城市起伏度，获得城市群高度空间信息。
优选地，采用SVM算法提取植被分布图层的步骤包括：在真彩色空间和假彩色空间分别建立一个最优分离面将植被与背景进行分离，分离面方程如下：
假彩色反射率:ω1×ρgreen+ω2ρred+ω3ρnir+b1＝0；
真彩色反射率:ω4×ρblue+ω5ρgreen+ω6ρred+b2＝0；
上式中ρblue、ρgreen、ρred、ρnir分别为blue、green、red、nir波段的反射率，ω1～ω6、b1、b2为系数。
优选地，OpenStreetMap路网数据的预处理步骤包括：
步骤S201，删除不必要的细节，修剪短于500m的道路，延长离要连接的线少于100m的道路；
步骤S202，生成缓冲区，根据标准和道路宽度的调查对路网的原始分类进行重分类，并对重分类的道路建立相应宽度的缓冲区。
优选地，在步骤S105中，点云数据的生成步骤包括：
步骤S301，用自相关二阶矩阵进行特征检测或者描述局部的影像结构；
步骤S302，用局部尺度为微分尺度的高斯核来计算图像特征点的方向梯度，再用尺度为积分尺度的高斯窗口来平滑特征点的邻域范围；
步骤S303，在求得的自相关二阶矩阵的特征值足够大时，将局部极大值点所在的像素检测为角点，之后根据角点响应函数来计算图像每个像素点的角点响应值，得到特征点集；
步骤S304，对提取的特征点集采用SIFT算法中的特征向量进行描述，形成Harris特征描述子，并采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置，之后用最小二乘影像匹配算法完成精确匹配，提取同名像点；
步骤S305，将同名像点与影像自带的RPC参数结合建立有理函数模型，进行三维坐标的解算生成3D的点云数据。
优选地，所述自相关二阶矩阵为：
Figure GDA0002593375210000041
其中，σ1为积分尺度，σD为微分尺度，σD＝SσD，Ix为图像在水平方向梯度的偏导数，Iy为图像在垂直方向梯度的偏导数，x,y为自相关矩阵对应的值；
所述角点响应函数为:
R＝Det(μ(x,y,σ1,σD))/αTrace2(μ(x,y,σ1,σD))+ε，
Det(M)＝λ1λ2＝AB-C2，Trace＝λ1+λ2＝A+B，
其中，Det为矩阵M的行列式的值，Trace为矩阵M的迹，λ1,λ2为自相关矩阵M的特征值，A为前视影像的某个关键点，B为A与正视或后视影像中欧式距离最近的关键点，C为A与正视或后视影像中欧式距离次近的关键点。
优选地，在步骤S304中，采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置的步骤包括：
取前视影像的某个关键点A，用最邻近法搜索关键点A与正视或后视影像中欧式距离最近和次近的两个关键点B、关键点C；在关键点B、关键点C中，如果最临近距离与次临近距离少于设定的比例阈值时，则接受关键点B、关键点C这一对匹配点。
优选地，在步骤S305中，建立有理函数模型包括利用影像附带的辅助参数，将像点坐标(r,c)表示为与相应地面点坐标(X,Y,Z)为自变量的多项式的比值，形式为：
Figure GDA0002593375210000051
其中(rn，cn)、(Xn,Yn,Zn)分别表示影像的像点坐标(r,c)与对应的地面点坐标(X,Y,Z)经过平移和缩放后的正则化地面坐标，Pi(Xn,Yn,Zn),i＝1,2,3,4为多项式，阶数不超过3。
优选地，所述城市生长高度的计算公式为：
ΔH＝H(t+1)-Ht，
其中Ht是指t时期内的城市建筑物高度，H(t+1)是指t+1时期内的城市建筑物高度。
优选地，所述城市起伏度的计算公式为：
R＝Hmax-Hmin，
其中Hmax是指在特定的区域内最高的城市建筑物，Hmin是指在特定的区域内最低的城市建筑物。
本发明的另一目的在于提供一种城市等高线的提取方法，包括采用所述提取方法获得城市建筑物信息及建筑物高度信息，根据城市建筑物信息及建筑物高度信息，将城市中建筑物高度相等的相邻各点连接成闭合曲线。
基于上述技术方案，本发明的优点是：
(1)本发明充分利用资源三号前视、后视和正视影像，创新性的构建了提取城市建筑物信息及建筑物高度信息的系统框架，并以此为基础获得城市等高线、城市起伏度和城市生长高度的模型，揭示了城市横向面积扩张和纵向高度增长的演变格局；
(2)针对城市建筑物的提取，本发明融合了遥感数据和社会感知数据，建立了地块尺度的高/低密度建筑区，结合多视角差值特征图，实现了城市建筑物的提取。
(3)针对城市建筑物高度的提取，本发明利用资源三号前视、后视和正视影像两两匹配生成密集点云；采用TIN线性插值法求取数字表面模型(DSM)，引入点云分类的思想，获取精确的地面点，建立数字高程模型(DEM)。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为城市群高度空间信息的提取方法步骤图；
图2为OpenStreetMap路网数据的预处理步骤图；
图3为点云数据的生成步骤图；
图4为城市群高度空间信息的提取方法流程图；
图5为真彩色反射率空间示意图；
图6为假彩色反射率空间示意图；
图7为实施例1中选取区域的植被图层；
图8为实施例1中选取区域的地块图层；
图9为实施例1中选取区域的多视角差值特征图；
图10为实施例1中选取区域的城市建筑物分布图；
图11为实施例1中选取区域的城市建筑物高度信息图；
图12为实施例1中选取区域的城市等高线分布图；
图13为实施例1中选取区域的城市起伏度分布图；
图14为实施例1中选取区域的城市生长高度图。
具体实施方式
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
本发明提供了一种城市群高度空间信息的提取方法，如图1～图6所示，其中示出了本发明的一种优选实施方式。
如图1所示，所述城市建筑物信息及建筑物高度信息的提取方法包括如下步骤：
步骤S101，将资源三号卫星多光谱影像用SVM算法提取出植被分布图层，再利用OpenStreetMap路网数据经过预处理得到地块分割图层；
步骤S102，将所述植被分布图层和所述地块分割图层进行叠加分析，得到每个地块的植被占比，并通过建立阈值判别体系得到高密集建筑区和低密度建筑区；
步骤S103，将资源三号卫星三线阵卫星影像用多视角差值法得到多视角差值特征图，并对所述多视角差值特征图进行去除阴影处理；
步骤S104，将高密集建筑区和低密度建筑区与去除阴影后的多视角差值特征图进行叠加分析，并通过建立阈值判别体系提取出城市建筑物信息；
步骤S105，将资源三号卫星三线阵卫星影像用计算连接点和区域网平差方法生成点云数据；
步骤S106，采用TIN线性插值法进行点云插值得到数字表面模型DSM，再对点云数据进行分类得到数字高程模型DEM，之后求取数字表面模型DSM与数字高程模型DEM的差值，得到nDSM；
步骤S107，通过步骤S104中提取出的建筑物信息获得城市建筑物图层，叠加城市建筑物图层和nDSM，采用掩膜处理提取城市建筑物高度信息；
步骤S108，根据城市建筑物信息及建筑物高度信息计算城市生长高度与城市起伏度，获得城市群高度空间信息。
资源三号卫星各项参数如下表：
Figure GDA0002593375210000081
本发明的提取方法充分利用资源三号前视、后视和正视影像，创新性的构建了提取城市建筑物信息及建筑物高度信息的系统框架，揭示了城市横向面积扩张和纵向高度增长的演变格局。
如图4所示，本发明的提取方法流程如下：
本发明的提取目标主要包括提取城市建筑物信息和城市建筑物高度信息两个主要的部分。针对城市建筑物的提取，首先利用ZY-3号多光谱影像，采用SVM算法提取植被分布图层；其次利用OSM路网数据得到地块分割图层，将上述两种图层叠加，通过判别地块植被占比得到高密集建筑区和低密度建筑区；然后利用ZY-3三线阵卫星影像，采用多视角差值法得到多视角差值特征图；最后分别将高密集建筑区和低密度建筑区与多视角差值特征图进行叠加分析，建立阈值判别体系，提取城市建筑物信息。
针对城市建筑物高度的提取，首先利用ZY-3三线阵卫星影像，采用计算连接点和区域网平差方法得到点云数据；其次对点云数据进行点云插值得到DSM，对点云分类得到DEM；最后结合城市建筑物数据，提取城市建筑物高度信息，进而分析城市等高线、城市起伏度和城市生长高度。
在机器学习领域，神经网络、支持向量机、决策树都可以用来很好地区分2类地物。本发明采用基于支持向量机算法SVM，提取区分植被和非植被与其它地物的最优分离面。SVM是在N维空间(Rn)空间中，寻找一个最大间隔超平面，因此，SVM也被称作最大间隔分类器。选用SVM用来构建植被和非植被与其它地物的最优分离面具有以下优势：(1)SVM产生的最优分离平面，有助于提升训练出的参数的稳定性；(2)SVM在较少样本条件下较其他分类器分类精度更高。SVM的基本思想是用一个最大间距的最优超平面，将一个输入向量分成2类，即求解下述方程：
Figure GDA0002593375210000091
Figure GDA0002593375210000092
式中：xi∈Rd代表训练的样本向量；yi∈{-1，+1}代表类别标签；K(u,v)代表核函数。
本发明在提取城市植被分布图层中，选用的是线性核，即K(xi,xj)＝xiTxj。此时，超平面的方程可以w×x-b＝0表示，SVM就是训练出超平面的法向量w和位移b的值，对于每个预测数据x带入该平面方程的正负来判断类别，x代表不同波段反射率的3维向量，而w对应每个向量前面的系数。
如图5、图6所示，由于植被跟低值背景(NDVI≤0.25)在假彩色空间具有很好的线性可分性(即可用一个平面进行分离)，植被跟高值背景(NDVI>0.25)在真彩色空间也具有很好的线性。因此，本发明在真彩色空间和假彩色空间分别建立一个最优分离面来实现植被与背景的分离，以达到提高植被提取精度的目的。
优选地，采用SVM算法提取植被分布图层的步骤包括：在真彩色空间和假彩色空间分别建立一个最优分离面将植被与背景进行分离，分离面方程如下：
假彩色反射率:ω1×ρgreen+ω2ρred+ω3ρnir+b1＝0；
真彩色反射率:ω4×ρblue+ω5ρgreen+ω6ρred+b2＝0；
上式中ρblue、ρgreen、ρred、ρnir分别为blue、green、red、nir波段的反射率，ω1～ω6、b1、b2为系数。
OpenStreetMap(简称OSM)是一个存储海量XML数据的数据库，来源于OSM的路网数据是线矢量文件，并包含道路类别属性。利用OSM的路网数据可以得到地块。为能够满足本发明提取方法的要求，需要对OpenStreetMap路网数据进行预处理。如图2所示，OpenStreetMap路网数据的预处理步骤包括：
步骤S201，删除不必要的细节，修剪短于500m的道路，延长离要连接的线少于100m的道路；
步骤S202，生成缓冲区，根据标准和道路宽度的调查对路网的原始分类进行重分类，并对重分类的道路建立相应宽度的缓冲区。
利用OSM路网数据得到地块分割图层，将植被分布图层和地块分割图层进行叠加分析得到每个地块的植被占比，建立阈值判别体系，得到高密集建筑区和低密度建筑区。
由于建筑物在多角度图像中(ZY-3号前视、后视、下视)呈现出明显的角差特征，而道路、土壤等地势较低的物体在多角度图像中没有明显的角差特征。因此，根据这一特征建立了角差特征公式如下：
ADF＝max(Xb-Xn,Xn-Xf,Xb-Xf)，其中，Xf，Xn，Xb分别代表ZY-3号前视，下视和后视影像。采用多视角差值法，本发明从资源三号卫星三线阵卫星影像中提取了多视角差值特征图。
在提取城市建筑物信息时，其步骤包括：首先对多视角差值特征图进行去除阴影的处理，再分别将高密集建筑区和低密度建筑区与多视角差值特征图进行叠加分析，建立阈值判别体系，最终提取城市建筑物信息。
由于Harris算子计算量小，有较好的稳定性，对图像灰度的平移、旋转，以及光照变化等具有不变性，不受摄像机姿态和光照的影响，提取角点分布均匀合理，本系统采用Harris算法提取特征点。如图3所示，点云数据的生成步骤包括：
步骤S301，用自相关二阶矩阵进行特征检测或者描述局部的影像结构；
步骤S302，用局部尺度为微分尺度的高斯核来计算图像特征点的方向梯度，再用尺度为积分尺度的高斯窗口来平滑特征点的邻域范围；
步骤S303，在求得的自相关二阶矩阵的特征值足够大时，将局部极大值点所在的像素检测为角点，之后根据角点响应函数来计算图像每个像素点的角点响应值，得到特征点集；
步骤S304，对提取的特征点集采用SIFT算法中的特征向量进行描述，形成Harris特征描述子，并采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置，之后用最小二乘影像匹配算法完成精确匹配，提取同名像点；
步骤S305，将同名像点与影像自带的RPC参数结合建立有理函数模型，进行三维坐标的解算生成3D的点云数据。
Harris算子是基于图像的自相关二阶矩阵，通常用自相关二阶矩阵进行特征检测或者描述局部的影像结构，找到局部极大值点作为特征点。所以，自相关二阶矩阵必须能够适应尺度的变化。优选地，所述自相关二阶矩阵为：
Figure GDA0002593375210000121
其中，σ1为积分尺度，σD为微分尺度，σD＝SσD，Ix为图像在水平方向梯度的偏导数，Iy为图像在垂直方向梯度的偏导数；上式描述了像素点局部邻域的梯度变化。用局部尺度为σD的高斯核来计算图像特征点的方向梯度，再用尺度为σ1的高斯窗口来平滑特征点的邻域范围。如果求得的自相关矩阵M的特征值λ1、λ2足够大时，就把该像素检测为角点。然后根据角点响应函数来计算图像每个像素点的角点响应值，所述角点响应函数为:
R＝Det(μ(x,y,σ1,σD))/αTrace2(μ(x,y,σ1,σD))+ε，
Det(M)＝λ1λ2＝AB-C2，Trace＝λ1+λ2＝A+B，其中，Det为矩阵M的行列式的值，Trace为矩阵M的迹，λ1,λ2为自相关矩阵M的特征值，A为前视影像的某个关键点，B为A与正视或后视影像中欧式距离最近的关键点，C为A与正视或后视影像中欧式距离次近的关键点。
对提取的特征点集使用SIFT算法中的特征向量进行描述，形成Harris特征描述子，采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置。具体地，采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置的步骤包括：取前视影像的某个关键点A，用最邻近法搜索关键点A与正视或后视影像中欧式距离最近和次近的两个关键点B、关键点C；在关键点B、关键点C中，如果最临近距离与次临近距离少于设定的比例阈值时，则接受关键点B、关键点C这一对匹配点。
之后用最小二乘影像匹配算法完成精确匹配，提取同名点。随后将同名像点与影像自带的RPC参数结合建立有理函数模型，进行三维坐标的解算生成3D点云。优选地，建立有理函数模型包括利用影像附带的辅助参数，将像点坐标(r,c)表示为与相应地面点坐标(X,Y,Z)为自变量的多项式的比值，形式为：
Figure GDA0002593375210000131
其中(rn，cn)、(Xn,Yn,Zn)分别表示影像的像点坐标(r,c)与对应的地面点坐标(X,Y,Z)经过平移和缩放后的正则化地面坐标，Pi(Xn,Yn,Zn),i＝1,2,3,4为多项式，阶数不超过3。
对于提取城市建筑物高度信息，其包括：根据上述步骤生成的点云数据，首先利用TIN线性插值法进行点云插值，得到数字表面模型(DSM)；采用随机森林的分类方法，根据目标对象的颜色和几何特征将点云分成植被、建筑、地面和道路四种类型；利用已分类的地面点和道路，采用TIN线性插值法进行点云插值，得到数字高程模型(DEM)；求取数字表面模型(DSM)与数字高程模型(DEM)的差值，得到nDSM；叠加城市建筑物图层和nDSM，利用掩膜处理的方法得到城市建筑物高度信息。
城市群高度空间信息包含城市起伏度和城市生长高度，因此可根据城市建筑物信息及建筑物高度信息计算城市生长高度与城市起伏度，获得城市群高度空间信息。
在本发明中，城市生长高度的是指同一区域不同时期内城市建筑物高度的变化。本发明还提供了一种城市生长高度的计算方法，包括采用所述提取方法获得城市建筑物信息及建筑物高度信息，根据城市建筑物信息及建筑物高度信息计算城市生长高度，其中城市生长高度的计算公式为：
ΔH＝H(t+1)-Ht；
其中，Ht是指t时期内的城市建筑物高度，H(t+1)是指t+1时期内的城市建筑物高度。
在本发明中，城市起伏度是指一个特定的区域内，最高的城市建筑物与最低的城市建筑物高度的差值。它是描述一个区域地形特征的一个宏观性的指标。本发明还提供了一种城市起伏度的计算方法，包括采用所述提取方法获得城市建筑物信息及建筑物高度信息，根据城市建筑物信息及建筑物高度信息计算城市起伏度，其中城市起伏度的计算公式为：
R＝Hmax-Hmin，
其中Hmax是指在特定的区域内最高的城市建筑物，Hmin是指在特定的区域内最低的城市建筑物。
城市等高线是指建筑物高度相等的相邻各点所连接成的闭合曲线。城市等高线可以直观的反映一个城市的纵向增长的空间格局。本发明提供了一种城市等高线的提取方法，包括采用所述提取方法获得城市建筑物信息及建筑物高度信息，根据城市建筑物信息及建筑物高度信息，将城市中建筑物高度相等的相邻各点连接成闭合曲线。
本发明的提取方法充分利用了资源三号前视、后视和正视影像，创新性的构建了提取城市建筑物信息及建筑物高度信息的系统框架，并以此为基础获得城市等高线、城市起伏度和城市生长高度的模型，揭示了城市横向面积扩张和纵向高度增长的演变格局；针对城市建筑物的提取，本发明融合了遥感数据和社会感知数据，建立了地块尺度的高/低密度建筑区，结合多视角差值特征图，实现了城市建筑物的提取。针对城市建筑物高度的提取，本发明利用资源三号前视、后视和正视影像两两匹配生成密集点云；采用TIN线性插值法求取数字表面模型(DSM)，引入点云分类的思想，获取精确的地面点，建立数字高程模型(DEM)。
实施例1
为进一步说明本发明的城市建筑物信息及建筑物高度信息的提取方法，下面以在北京市选取部分区域作为研究区进行说明。
首先利用SVM算法得到如图7所示的植被图层，再利用道路数据得到如图8所示的地块图层；
之后利用多视角差值法得到多视角差值特征图，从图9可以看出建筑物与其他地类在多视角差值特征图中表现出了明显的差异性，建筑物整体比较明亮，越高的区域颜色越明亮；图10是多视角差值特征图经过去除阴影的处理之后，叠加高/低密度建筑区以及阈值分析得到的最终的城市建筑物分布。
针对建筑物高度的提取，利用Harris算子、欧氏距离等算法提取同名点，构建有理函数模型生成点云，然后分别对点云进行插值和分类后插值得到DSM和DEM数据；根据两者的差值和城市建筑物图层，得到城市建筑物高度，如图11所示。根据本系统城市等高线的定义得到城市等高线分布图，图12中大部分等高线都比较稀疏，表明该区域不同建筑物之间高度差异较小，城市天际线与地面接近水平状态；局部地区出现密集的等高线，表明该区域不同建筑物之间高度差异较大，密集中心区出现了高层建筑。
根据城市起伏度的定义得到城市起伏度分布图，如图13所示。如图14所示，根据2018年和2013年两期的建筑物高度数据得到该区域近5年的城市生长高度图。
最后应当说明的是:以上实施例仅用以说明本发明的技术方案而非对其限制；尽管参照较佳实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者对部分技术特征进行等同替换；而不脱离本发明技术方案的精神，其均应涵盖在本发明请求保护的技术方案范围当中。

Claims (10)

1.一种城市群高度空间信息的提取方法，其特征在于：所述提取方法包括如下步骤：
步骤S101，将资源三号卫星多光谱影像用SVM算法提取出植被分布图层，再利用OpenStreetMap路网数据经过预处理得到地块分割图层；
步骤S102，将所述植被分布图层和所述地块分割图层进行叠加分析，得到每个地块的植被占比，并通过建立阈值判别体系得到高密集建筑区和低密度建筑区；
步骤S103，将资源三号卫星三线阵卫星影像用多视角差值法得到多视角差值特征图，并对所述多视角差值特征图进行去除阴影处理；
步骤S104，将高密集建筑区和低密度建筑区与去除阴影后的多视角差值特征图进行叠加分析，并通过建立阈值判别体系提取出城市建筑物信息；
步骤S105，将资源三号卫星三线阵卫星影像用计算连接点和区域网平差方法生成点云数据；
步骤S106，采用TIN线性插值法进行点云插值得到数字表面模型DSM，再对点云数据进行分类得到数字高程模型DEM，之后求取数字表面模型DSM与数字高程模型DEM的差值，得到nDSM；
步骤S107，通过步骤S104中提取出的建筑物信息获得城市建筑物图层，叠加城市建筑物图层和nDSM，采用掩膜处理提取城市建筑物高度信息；
步骤S108，根据城市建筑物信息及建筑物高度信息计算城市生长高度与城市起伏度，获得城市群高度空间信息。
2.根据权利要求1所述的提取方法，其特征在于：采用SVM算法提取植被分布图层的步骤包括：在真彩色空间和假彩色空间分别建立一个最优分离面将植被与背景进行分离，分离面方程如下：
假彩色反射率:ω1×ρgreen+ω2ρred+ω3ρnir+b1＝0；
真彩色反射率:ω4×ρblue+ω5ρgreen+ω6ρred+b2＝0；
上式中ρblue、ρgreen、ρred、ρnir分别为blue、green、red、nir波段的反射率，ω1～ω6、b1、b2为系数。
3.根据权利要求1所述的提取方法，其特征在于：OpenStreetMap路网数据的预处理步骤包括：
步骤S201，删除不必要的细节，修剪短于500m的道路，延长离要连接的线少于100m的道路；
步骤S202，生成缓冲区，根据标准和道路宽度的调查对路网的原始分类进行重分类，并对重分类的道路建立相应宽度的缓冲区。
4.根据权利要求1所述的提取方法，其特征在于：在步骤S105中，点云数据的生成步骤包括：
步骤S301，用自相关二阶矩阵进行特征检测或者描述局部的影像结构；
步骤S302，用局部尺度为微分尺度的高斯核来计算图像特征点的方向梯度，再用尺度为积分尺度的高斯窗口来平滑特征点的邻域范围；
步骤S303，在求得的自相关二阶矩阵的特征值足够大时，将局部极大值点所在的像素检测为角点，之后根据角点响应函数来计算图像每个像素点的角点响应值，得到特征点集；
步骤S304，对提取的特征点集采用SIFT算法中的特征向量进行描述，形成Harris特征描述子，并采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置，之后用最小二乘影像匹配算法完成精确匹配，提取同名像点；
步骤S305，将同名像点与影像自带的RPC参数结合建立有理函数模型，进行三维坐标的解算生成3D的点云数据。
5.根据权利要求4所述的提取方法，其特征在于：所述自相关二阶矩阵为：
Figure FDA0002593375200000031
其中，σ1为积分尺度，σD为微分尺度，σD＝SσD，Ix为图像在水平方向梯度的偏导数，Iy为图像在垂直方向梯度的偏导数，x,y为自相关矩阵对应的值；
所述角点响应函数为:
R＝Det(μ(x,y,σ1,σD))/αTrace2(μ(x,y,σ1,σD))+ε，
Det(M)＝λ1λ2＝AB-C2，Trace＝λ1+λ2＝A+B，
其中，Det为矩阵M的行列式的值，Trace为矩阵M的迹，λ1,λ2为自相关矩阵M的特征值，A为前视影像的某个关键点，B为A与正视或后视影像中欧式距离最近的关键点，C为A与正视或后视影像中欧式距离次近的关键点。
6.根据权利要求4所述的提取方法，其特征在于：在步骤S304中，采用欧氏距离测度找到一幅影像的特征点在另一幅影像中对应的位置的步骤包括：
取前视影像的某个关键点A，用最邻近法搜索关键点A与正视或后视影像中欧式距离最近和次近的两个关键点B、关键点C；在关键点B、关键点C中，如果最临近距离与次临近距离少于设定的比例阈值时，则接受关键点B、关键点C这一对匹配点。
7.根据权利要求4所述的提取方法，其特征在于：在步骤S305中，建立有理函数模型包括利用影像附带的辅助参数，将像点坐标(r,c)表示为与相应地面点坐标(X,Y,Z)为自变量的多项式的比值，形式为：
Figure FDA0002593375200000032
其中(rn，cn)、(Xn,Yn,Zn)分别表示影像的像点坐标(r,c)与对应的地面点坐标(X,Y,Z)经过平移和缩放后的正则化地面坐标，Pi(Xn,Yn,Zn),i＝1,2,3,4为多项式，阶数不超过3。
8.根据权利要求1所述的提取方法，其特征在于：所述城市生长高度的计算公式为：
ΔH＝H(t+1)-Ht，
其中Ht是指t时期内的城市建筑物高度，H(t+1)是指t+1时期内的城市建筑物高度。
9.根据权利要求1所述的提取方法，其特征在于：所述城市起伏度的计算公式为：
R＝Hmax-Hmin，
其中Hmax是指在特定的区域内最高的城市建筑物，Hmin是指在特定的区域内最低的城市建筑物。
10.一种城市等高线的提取方法，其特征在于：包括采用上述权利要求1～9任一项中所述提取方法获得城市建筑物信息及建筑物高度信息，根据城市建筑物信息及建筑物高度信息，将城市中建筑物高度相等的相邻各点连接成闭合曲线。
CN201911221885.7A 2019-12-03 2019-12-03 一种城市群高度空间信息及等高线的提取方法 Active CN111047695B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911221885.7A CN111047695B (zh) 2019-12-03 2019-12-03 一种城市群高度空间信息及等高线的提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911221885.7A CN111047695B (zh) 2019-12-03 2019-12-03 一种城市群高度空间信息及等高线的提取方法

Publications (2)

Publication Number Publication Date
CN111047695A CN111047695A (zh) 2020-04-21
CN111047695B true CN111047695B (zh) 2020-11-10

Family

ID=70233415

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911221885.7A Active CN111047695B (zh) 2019-12-03 2019-12-03 一种城市群高度空间信息及等高线的提取方法

Country Status (1)

Country Link
CN (1) CN111047695B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111598048B (zh) * 2020-05-31 2021-06-15 中国科学院地理科学与资源研究所 一种融合高分遥感影像和街景影像的城中村识别方法
CN112017285B (zh) * 2020-08-28 2021-09-17 北京国遥新天地信息技术股份有限公司 一种三维gis中实时使地形精确贴合条状模型的方法
CN112598562A (zh) * 2020-12-30 2021-04-02 东南大学建筑设计研究院有限公司 一种基于地块罩面的城市高度精细化管控方法
CN113033484B (zh) * 2021-04-21 2022-11-22 河北工程大学 一种面向无人机应急网络部署的城市分类方法
CN113358091B (zh) * 2021-06-02 2021-11-16 自然资源部国土卫星遥感应用中心 一种利用三线阵立体卫星影像生产数字高程模型dem的方法
FR3123750A1 (fr) * 2021-06-04 2022-12-09 Murmuration Système et procédé pour prédire l’impact environnemental d’au moins une structure à construire sur au moins une zone géographique prédéterminée de la terre
CN113920266B (zh) * 2021-11-03 2022-10-21 泰瑞数创科技(北京)股份有限公司 城市信息模型语义信息人工智能生成方法和系统
CN115100536B (zh) * 2022-06-01 2023-03-28 中科星睿科技(北京)有限公司 建筑物识别方法、装置、电子设备和计算机可读介质
CN115471634B (zh) * 2022-10-28 2023-03-24 吉奥时空信息技术股份有限公司 一种城市绿植孪生的建模方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104463868A (zh) * 2014-12-05 2015-03-25 北京师范大学 一种基于无参数高分辨率影像的建筑高度快速获取方法
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN108871286A (zh) * 2018-04-25 2018-11-23 中国科学院遥感与数字地球研究所 空间大数据协同的城市建成区人口密度估算方法和系统
CN109579784A (zh) * 2018-11-26 2019-04-05 青岛国测海遥信息技术有限公司 基于数字表面模型的城市地区建筑物高度的自动获取方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104463868A (zh) * 2014-12-05 2015-03-25 北京师范大学 一种基于无参数高分辨率影像的建筑高度快速获取方法
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN108871286A (zh) * 2018-04-25 2018-11-23 中国科学院遥感与数字地球研究所 空间大数据协同的城市建成区人口密度估算方法和系统
CN109579784A (zh) * 2018-11-26 2019-04-05 青岛国测海遥信息技术有限公司 基于数字表面模型的城市地区建筑物高度的自动获取方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"On the validation of DEM and FEM/DEM models in 2D and 3D";Jiansheng Xiang等;《Engineering Computations》;20090831;全文 *
"基于规则格网DEM的等高线提取算法的优化与实现";边淑莉等;《测绘科学》;20081130;第33卷(第6期);全文 *
"无参数高分辨率遥感影像的建筑高度";乔伟峰等;《地球信息科学》;20150831;第17卷(第8期);全文 *

Also Published As

Publication number Publication date
CN111047695A (zh) 2020-04-21

Similar Documents

Publication Publication Date Title
CN111047695B (zh) 一种城市群高度空间信息及等高线的提取方法
CN112489212B (zh) 一种基于多源遥感数据的建筑物智能化三维测图方法
CN108573276B (zh) 一种基于高分辨率遥感影像的变化检测方法
CN110263717B (zh) 一种融入街景影像的土地利用类别确定方法
CN108197583B (zh) 基于图割优化和影像结构特征的建筑物变化检测方法
CN109598794B (zh) 三维gis动态模型的构建方法
CN111709981A (zh) 特征线融合的激光点云与模拟图像的配准方法
CN106780524A (zh) 一种三维点云道路边界自动提取方法
WO2018061010A1 (en) Point cloud transforming in large-scale urban modelling
CN110084205A (zh) 一种基于改进的面向对象高分辨率遥感影像分类方法
CN111652241B (zh) 融合影像特征与密集匹配点云特征的建筑物轮廓提取方法
CN110033484A (zh) 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法
CN111323788B (zh) 建筑物变化的监测方法、装置及计算机设备
CN111487643B (zh) 一种基于激光雷达点云和近红外影像的建筑物检测方法
CN112241661A (zh) 一种结合机载LiDAR点云数据和航空影像的城市地物精细化分类方法
CN109492606A (zh) 多光谱矢量图获取方法及系统、三维单体化方法及系统
Wang et al. 3D road information extraction from LiDAR data fused with aerial-images
CN115222884A (zh) 一种基于人工智能的空间对象分析及建模优化方法
CN111458691B (zh) 建筑物信息的提取方法、装置及计算机设备
CN112946679A (zh) 一种基于人工智能的无人机测绘果冻效应检测方法及系统
CN116863357A (zh) 一种无人机遥感堤坝影像定标与智能分割变化检测方法
CN109671109A (zh) 密集点云生成方法及系统
CN115984721A (zh) 基于倾斜摄影和图像识别技术实现乡村风貌管理的方法
Sun et al. Building outline extraction from aerial imagery and digital surface model with a frame field learning framework
CN118015494A (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
GR01 Patent grant
GR01 Patent grant