CN111651545A - 一种基于多源数据融合的城市边缘区提取方法 - Google Patents

一种基于多源数据融合的城市边缘区提取方法 Download PDF

Info

Publication number
CN111651545A
CN111651545A CN202010493552.6A CN202010493552A CN111651545A CN 111651545 A CN111651545 A CN 111651545A CN 202010493552 A CN202010493552 A CN 202010493552A CN 111651545 A CN111651545 A CN 111651545A
Authority
CN
China
Prior art keywords
mutation
urban
cluster
point
mutation 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.)
Granted
Application number
CN202010493552.6A
Other languages
English (en)
Other versions
CN111651545B (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.)
Nanjing Tushentu Information Technology Co ltd
Original Assignee
Nanjing Tushentu Information Technology Co ltd
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 Nanjing Tushentu Information Technology Co ltd filed Critical Nanjing Tushentu Information Technology Co ltd
Priority to CN202010493552.6A priority Critical patent/CN111651545B/zh
Publication of CN111651545A publication Critical patent/CN111651545A/zh
Application granted granted Critical
Publication of CN111651545B publication Critical patent/CN111651545B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/148—Wavelet transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A30/00—Adapting or protecting infrastructure or their operation
    • Y02A30/60—Planning or developing urban green infrastructure

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Human Resources & Organizations (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Strategic Management (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Remote Sensing (AREA)
  • Geometry (AREA)
  • Development Economics (AREA)
  • Educational Administration (AREA)
  • Economics (AREA)
  • Computer Graphics (AREA)
  • Game Theory and Decision Science (AREA)
  • Algebra (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Traffic Control Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于多源数据融合的城市边缘区提取方法，该方法选取多项城市边缘特征数据作为数据源，并分别确定各项城市边缘特征数据的城市边缘区识别指标因子；对各项城市边缘区识别指标因子分别进行数据标准化处理后，计算各项城市边缘区识别指标因子的指标值和指标权重，并通过指标值和指标权重计算得到城市化特征值；然后根据城市化特征值进行城市边缘区突变点识别和提取，得到城市边缘区突变点集群；通过对城市边缘区突变点集群进行空间位置和属性特征聚类，并对每个聚类簇进行内外边界提取，得到城市边缘区边界。采用本发明的基于多源数据融合的城市边缘区提取方法，实现了城市边缘区边界自动提取，且提高了边缘区边界识别精度。

Description

一种基于多源数据融合的城市边缘区提取方法
技术领域
本发明涉及城市规划技术领域，具体而言，涉及一种基于多源数据融合的城市边缘区提取方法。
背景技术
城市边缘区位于城市与乡村之间的过渡交接地带，是城市扩张的前沿，也是城乡建设和用地置换中最具活力的地区。准确地识别城市边缘区，有助于从城乡对比的角度来衡量城市化程度及其对生态环境的影响，对于推进存量规划下的城市更新、空间优化、城市增长边界精细划定、多情景下的用地模拟等有着重要的现实意义。
城市边缘区的空间范围在早期是通过定性的方式来确定，较为常见的是以距中心城区的距离进行划分，随着遥感与GIS的快速发展，通过定量的方式识别城市边缘区空间范围的研究也日益受到重视。现有技术中，对于城市边缘区识别提取方法中的空间数据类型，常见的主要有遥感影像数据、城市空间轨迹数据和社会经济数据等。现有的城市边缘区识别提取方法主要分为综合分析法、阈值法、突变检测法等。其中，综合分析法是通过对区域内的土地利用、人口、经济等特征表现进行界定，如有学者人通过综合分析法分析各项指标梯度在城市空间的情况从而对城市边缘区进行划分，但此种方法较为依赖人的主观性；阈值法是通过距离、密度等测度的值域范围进行判断，如结合信息熵理论计算出城市边缘区内外边界的阈值区间范围，但算法中依旧存在主观判断的问题；而突变检测法是通过计算指标突变值进行边界划分，有的基于小波模极大值检测原理检测土地利用程度综合指数的突变点群，进而识别城乡边缘带边界；有的采用滑动t检验突变检测方法检测夜间灯光数据的突变点，提取了城市边缘区的内外边界；有的寻找建设用地比率指标的突变点，划分了城市城乡过渡带；还有通过确定土地利用动态度的高低值转折点来划定城乡交错带边界等。
现有技术方法仍存在以下的问题：
1、缺少多源数据融合的综合指标定量化判别方法
现有对城市边缘区范围识别主要是利用土地利用数据对空间范围进行识别，辅之以社会经济数据对得到的空间范围进行修正，而这些单一指标不能准确反映城市边缘的内部非均质性，导致识别精度不高。基于遥感影像数据、基于社会经济统计数据和基于城市空间轨迹数据等方法，虽然某一指标精度较高，但不能全面反映于城市边缘区多方面特征，造成识别结果误差较大。因此需要融合多源数据并建立边缘区综合指标，通过不同来源的数据设定不同指标因子权重，从而对城市边缘区范围量化精准识别。
2、现有的边界识别提取方法在城市边缘区识别中存在局限性
现有的边界判别法主要分为综合分析法、阈值法、突变检测法等，其中，综合分析法是通过对区域内的土地利用、人口、经济等特征表现进行界定；阈值法是通过距离、密度等测度的值域范围进行判断；突变检测法是计算指标突变值进行边界划分。以上方法，综合分析法和阈值法的人为主观性较强，而突变检测法的不足主要包括两方面：一是计算突变点时，通常会将突变性最大的两个位置作为城市边缘区的内外边界，符合理想化的“城市建成区——城市边缘区——乡村腹地”呈“同心圆”的发展形态，但现实是，由于建设开发强度、人口集聚等影响，城市边缘区并不是闭合的环状形态，带来了突变点个数的不确定性；二是识别出突变点后，通常是按就近原则直接手动连接相邻突变点获得边界，导致边界结果不准确，同时顾及空间和属性的自动提取边界的方法较少。
发明内容
为解决上述问题，本发明通过设计城市边缘区识别指标因子，采用合适的方法对各项指标进行量化，得到一个综合的城市化特征值，提供了一种有效的城市边缘区突变点识别方法以及城市边缘区边界自动化提取方法。本发明融合多类型的多源数据城市化特征指标因子，并在统一标准化确定各指标因子权重，得到用于判别边缘区特征值综合指标，本发明充分考虑了城市边缘区的过渡性、动态性和非均质性的特点。在此基础上，对于常见的边缘区边界识别提取方法的不足，本发明提出了双约束条件下的Delaunay三角网对城市边缘区突变点群进行空间和属性聚类，并利用约束聚类簇边界识别提取算法对每个聚类簇进行边界提取，最终得到城市边缘区边界。
为实现上述目的，本发明提供了一种基于多源数据融合的城市边缘区提取方法，其特征在于，该方法包括以下步骤：
S1，选取多项城市边缘特征数据作为数据源，并分别确定各项城市边缘特征数据的城市边缘区识别指标因子；
S2，对步骤S1确定的各城市边缘区识别指标因子分别进行数据标准化处理，计算各城市边缘区识别指标因子的指标值和指标权重，并通过上述指标值和指标权重计算得到城市化特征值；
S3，根据步骤S3计算得到的城市化特征值进行城市边缘区突变点识别和提取，得到城市边缘区突变点集群；
S4，通过对步骤S3得到的城市边缘区突变点集群进行聚类处理，并对聚类处理得到的每个聚类簇进行边界提取，得到城市边缘区边界。
作为本发明进一步的改进，步骤S1中，所述多项城市边缘特征数据包括：土地利用数据、道路数据、夜间灯光数据和POI数据，所述多项城市边缘特征数据对应的城市边缘区识别指标因子分别为确定建设用地密度、道路密度、夜间灯光强度和服务型POI密度。
作为本发明进一步的改进，步骤S1中，确定建设用地密度、道路密度、夜间灯光强度和服务型POI密度的具体方式如下，其中，
确定建设用地密度，包括：将所述土地利用数据分为建设用地、非建设用地和水域三大类，并与一格网做相交运算，得到格网中每个网格的土地利用数据信息，并对其进行格网统计，每个网格中的建设用地密度采用如下公式计算，
Figure BDA0002521976620000031
式中，buildDensityi为第i个网格中的建设用地密度，∑Si为第i个网格中的建设用地面积之和，S0为网格面积；
确定道路密度计算，包括：对所述道路数据进行了预处理，剔除国道、省道和高速公路，保留快速路、主干路、次干路和支路，将道路数据与一格网进行相交运算，并分类统计每个网格对应的道路数据，得到道路数据的格网统计结果，每个网格中的道路密度采用如下公式计算，
Figure BDA0002521976620000032
式中，roadDensityi为第i个网格中的道路密度，∑Li为第i个网格中的道路长度之和，S0为网格面积；
确定夜间灯光强度计算，包括：对所述夜间灯光数据进行重采样预处理，将预处理后的夜间灯光数据与一用格网进行相交运算，用格网对夜间灯光数据进行分区统计，计算每个网格中夜间灯光数据的平均值，得到格网统计下的每个网格的夜间灯光强度；
确定服务型POI密度计算，包括：对所述POI数据进行坐标转换，并与一格网做相交运算，得到每个网格中的POI数据统计结果，每个网格中POI密度采用如下公式计算，
Figure BDA0002521976620000041
式中，POIDensityi为第i个网格中的服务型POI密度，∑Ni为第i个网格中的POI数量，S0为网格面积。
确定建设用地密度、道路密度、夜间灯光强度和服务型POI密度均采用同一尺度下的格网进行相交运算。
作为本发明进一步的改进，步骤S2中，数据标准化处理采用极差法对所述各城市边缘区识别指标因子做进行无量纲化处理，并通过采用熵权法确定所述各城市边缘区识别指标因子的各指标权重。
作为本发明进一步的改进，步骤S3中，基于小波分析的突变点检测方法，对得到的城市化特征值进行城市边缘区突变点识别与提取，具体步骤包括：
S301，提取城市化特征值空间序列曲线：
基于各项数据源，获取城市化特征图层；
以城市中心为起点，城市外围的网格中心为终点作多条采样线，得到采样线图层，并在同方向上的采样线上等距离取点，以使采样线经过的格网中的每个网格中至少有一个采样点；
将采样线图层与城市化特征值图层做相交运算，得到城市化特征值的点图层；
在所述点图层中，确定对每条采样线上从所述城市中心至所述网格中心的各采样点对应的城市化特征值，得到城市化特征值空间序列曲线；
S302，
确定所述城市化特征值空间序列曲线的小波系数：
通过db1小波基函数对所述城市化特征值空间序列曲线进行多尺度小波变换；
确定所述多尺度小波变换后各采样线所对应的各小波系数模极大值；
S303，提取所述城市化特征值空间序列曲线上的突变点：
对每条采样线，将对应的小波系数模极大值位置作为候选突变点，通过两倍标准差法剔除所述候选突变点中的伪突变点；
将所述候选突变点中剔除所述伪突变点后剩余的突变点映射到二维空间中该采样线所经过的网格中心；
依次类推，提取所述各采样线对应的突变点；S304，生成城市边缘区突变点图层：
将城市化特征值空间序列曲线上所述各采样线的突变点映射为网格中心；
对映射的网格中心进行去重处理，得到城市边缘区突变点集群。
作为本发明进一步的改进，步骤S4包括：
S401，构建Delaunay三角网：对城市边缘区突变点群构建Delaunay三角网；
S402，基于空间位置的聚类：对所述Delaunay三角网进行双阈值规则约束，得到所述城市边缘区突变点群中各突变点的空间邻近关系；
S403，基于属性特征的聚类：在空间邻近关系的基础上，基于属性权重法进行属性聚类，对所述城市边缘区突变点群中各突变点的连接边间的连接边进一步约束，得到多个聚类簇；
S404，聚类簇边界提取：基于约束聚类簇边界识别提取算法识别聚类簇的内外边界，并对每个聚类簇进行内外边界提取，得到城市边缘区边界。
作为本发明进一步的改进，步骤402中，规则约束采用双阈值方法对所述Delaunay三角网的冗余边进行剔除，
得到城市边缘区突变点群间中各突变点间的空间邻近关系，其中包括：
对所述Delaunay三角网中的全局冗余边和局部冗余边进行剔除，得到城市边缘区突变点群间的空间邻近关系；，
其中，
全局冗余边的约束规则为：
Figure BDA0002521976620000061
式中，Global_LongEdges(p)为p点处需要删除的全局冗余边的集合，GlobalMean为三角网的所有边的平均长度，PartialMean(p)为与p点直接相连边的平均长度，GlobalSD为三角网的所有边长的标准差，ei为任意的点p与其直接相连的边；
局部冗余边的约束规则为：
Figure BDA0002521976620000062
式中，Partial_LongEdges(p)为p点处需要删除的局部冗余边的集合，PartialMean 2(p)为p点2阶邻域内所有边长的平均值，PartialSD(Pj)为与p点直接相连边的标准差，ei为与p相连的2阶邻域内的边。
作为本发明进一步的改进，步骤S403中，基于属性权重法进行属性聚类，对每一个城市边缘区突变点进行如下操作步骤：
a1，根据城市边缘区突变点的属性对聚类结果进行分类，并计算每一个城市边缘区突变点的邻域属性权值，计算公式如下所示，
Figure BDA0002521976620000063
其中，DAEnei(O)为突变点O的邻域属性权值，Eoc表示突变点O与聚类簇C的属性相似性，聚类簇C由n个点{C1,C2,C3,...,Cn}组成，聚类簇C是中心突变点一阶邻域范围内的突变点集；
a2，选取邻域属性权值最大的城市边缘区突变点作为起始突变点，将起始突变点作为中心突变点，将中心突变点周围的每个突变点视为一个单独的聚类簇C，计算中心突变点与周围每个突变点之间的属性相似性Eoc，由属性权值最大的突变点和周围属性相似度Eoc最大的突变点组成初始聚类簇，将初始聚类簇一阶邻域内的突变点设为候选突变点，利用如下公式继续度量每个候选突变点与初始聚类簇间的属性相似性；
Figure BDA0002521976620000064
Figure BDA0002521976620000065
式中，Eoc为中心点与周围每个点之间的属性相似性，vi为聚类簇中每个突变点的属性值，vn+1为中心突变点的属性值；
a3，对属性相似性进行量化，若突变点与聚类簇之间的相似性量化值大与预设阈值，则将该突变点加入到聚类簇C中，否则，跳过该突变点，重复步骤a2，继续扩大聚类簇，直至聚类簇的一阶邻域内没有类似的突变点，得到第一个聚类簇；
a4，在其余的突变点中将属性权值最大的突变点作为第二簇的起始点，重复步骤a2和a3，使所有城市边缘区突变点群中的各突变点分成不同的聚类簇。
作为本发明进一步的改进，步骤a3中，属性相似性量化方式为：
Figure BDA0002521976620000071
Figure BDA0002521976620000072
式中，θ为用于比较相似度Eoc的标准化参数，
Figure BDA0002521976620000073
为突变点O与聚类簇C间的最大属性权值，
其中，标准化参数的取值通过PBM指数确定，所述PBM指数采用如下公式确定：
Figure BDA0002521976620000074
Figure BDA0002521976620000075
式中，Nc为聚类簇的个数，Ni为聚类簇Ci中的突变点个数，zj为聚类簇Ci的质心坐标，xj为第j个突变点的坐标；
当PBM值最大时，得到最优θ值。
作为本发明进一步的改进，步骤S404中，约束聚类簇边界识别提取方法包括轮廓边界提取和孔洞边界提取，并对每个聚类簇进行内外边界提取，得到城市边缘区边界，包括：
基于轮廓边界提取得到城市边缘区突变点集群的外边界；其中，所述提取所述城市边缘区突变点集群的外边界包括：
判断所述Delaunay三角网最外层的边界边对应的角度是否大于阈值R，若所述边界边对应的角度大于阈值R，则删除所述边界边，依次类推，得到所述城市边缘区突变点集群的外边界；
基于孔洞边界提取得到城市边缘区突变点集群的内边界；其中，基于所述城市边缘区突变点集群的外边界提取所述城市边缘区突变点集群的内边界，包括：
计算每个聚类簇的松散度，将每个聚类簇中的突变点根据松散度进行排序，基于松散度阈值将每个聚类簇划分为两个子簇，提取两个子簇中松散度大于所述松散度阈值的突变点，其中，所述松散度阈值通过PBM指数确定；计算每个聚类簇的各突变点的冗余边准则，若突变点p的第i条边大于预设阈值时，则将该边删除，得到所述城市边缘区突变点集群的内边界；
其中，松散度采用如下公式计算，
F(p)＝PartialSD(p)/PartialMean(p)
式中，p为突变点，F(p)为突变点p的松散度，PartialMean(p)为与p相连的边界边长的平均值，PartialSD(p)为与p相连的边界边长的标准差；
冗余边准则采用如下公式计算
T(p)＝{ei|ei＞PartialMean(p)+GlobalSD}
式中，p为突变点，T(p)为突变点p的冗余边准则，PartialMean(p)为与p相连的边界边长的平均值，GlobalSD为所有边的边长标准差，ei为与p相连的第i条边的长度。
本发明的有益效果为：
1、采用多数据源建立城市边缘区识别综合指标因子，融合多源异构数据对城市边缘区进行定量化判别，利用多指标相结合的方式，更准确反映城市边缘区异质性特征，提高了城市边缘区识别精度。
2、本发明综合考虑了突变点群的空间位置邻近关系和属性特征相似性，利用约束型Delaunay三角网聚类方法进行点群聚类，采用约束聚类簇边界识别算法提取出聚类簇边界，从而实现对边缘区突变点群边界的自动化提取，本发明也适用于非闭合环装的城市边缘区突变点的识别，有效避免了人为主观因素的影响，提高了突变点识别精度，进一步提高了城市边缘区边界的识别精确度。
附图说明
图1为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的流程图；
图2为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的图城市边缘区识别流程示意图；
图3为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的城市化特征值空间序列曲线构建过程示意图；
图4为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的采样线的城市化特征值空间序列曲线示意图；
图5为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的采样线对应城市化特征值空间序列曲线小波变换分析示意图；
图6为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的采样线对应特征值曲线突变点检测结果示意图；
图7为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的城市边缘区突变点集群示意图；
图8为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的双约束条件下的突变点群边界提取流程示意图；
图9为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的删除全局冗余边后聚类簇分布情况示意图；
图10为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的删除局部冗余边后聚类簇分布情况示意图；
图11为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的基于约束型Delaunay三角网聚类的聚类结果示意图；
图12为本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法的实验数据城市边缘区提取结果示意图。
具体实施方式
下面通过具体的实施例并结合附图对本发明做进一步的详细描述。
本发明实施例所述的一种基于多源数据融合的城市边缘区提取方法，其特征在于，该方法包括以下步骤：
S1，选取多项城市边缘特征数据作为数据源，并分别确定各项城市边缘特征数据的城市边缘区识别指标因子；
S2，对步骤S1确定的各项城市边缘区识别指标因子分别进行数据标准化处理后，计算各项城市边缘区识别指标因子的指标值和指标权重，并通过上述指标值和指标权重计算得到城市化特征值；
S3，根据步骤S3计算得到的城市化特征值进行城市边缘区突变点识别和提取，得到城市边缘区突变点集群；
S4，通过对步骤S3得到的城市边缘区突变点集群进行空间位置和属性特征聚类，并对每个聚类簇进行内外边界提取，得到城市边缘区边界。
进一步的，由于城市边缘区地区土地价格低廉，土地开发利用粗放低效，建设用地密度低于城区但高于乡村，因此利用建设用地密度指标可直接反映区域土地利用结构和开发强度。道路密度能够反映区域内道路的密集程度，密度越大说明道路对区域发展的承载能力越高，城市中心性越强，城市边缘区的道路密度低于城区但高于乡村，利用道路密度指标可直接反映区域交通通达度。夜间灯光数据可以反映和衡量区域经济活动的空间分布，并且突破了行政区划的限制。城市边缘区的夜间灯光强度低于城区但高于乡村，利用夜间灯光强度可以间接反映区域人口和经济要素分布状况。由遥感影像解译得到的土地利用数据，由于其精度的限制，很难对建设用地类型进一步划分，不同类型的建设用地在城市中的分布密度是不同的，而夜间灯光亮度来自夜间活动，但夜间的活跃程度与产业结构、风俗文化等社会因素有巨大关系，并难以揭示丰富的经济结构内涵。上述数据的缺陷，可以通过POI数据进行弥补，在现代社会中，居民生活不可避免地产生对不同类型服务的需求，推动不同服务实体(例如加油站，便利店，学校和医院)的出现。POI数据中存在大量服务设施点数据，服务设施点的数量与居民数量密切相关。城市边缘区的服务型POI密度低于城区但高于乡村，利用服务型POI密度可以间接反映人群活动强度。因此，本实施例选择待提取城市边缘区域的土地利用数据、道路数据、夜间灯光数据和POI数据作为数据源，确定建设用地密度、道路密度、夜间灯光强度和服务型POI密度作为城市边缘区识别指标因子。其中，
建设用地密度计算：将土地利用数据分为建设用地、非建设用地和水域三大类，并与某一尺度下的格网做相交运算，得到格网中每个网格的土地利用数据信息，并对其进行格网统计，每个网格中的建设用地密度采用式(1)计算，
Figure BDA0002521976620000111
式中，buildDensityi为第i个网格中的建设用地密度(m2/m2)，∑Si为第i个网格中的建设用地面积之和，S0为网格面积；
道路密度计算：在计算道路密度之前，对道路数据进行了预处理，剔除了国道、省道和高速公路等市外交通，保留了快速路、主干路、次干路和支路等市内交通，将道路数据与某一尺度下的格网进行相交运算，并分类统计每个网格对应的道路数据，得到道路数据的格网统计结果，每个网格中的道路密度采用式(2)计算，
Figure BDA0002521976620000112
式中，roadDensityi为第i个格网中的道路密度，∑Li为第i个格网中的道路长度之和，S0为格网面积；
夜间灯光强度计算：对夜间灯光数据进行重采样预处理工作，用某一尺度下的格网对夜间灯光栅格数据进行分区统计，计算每个网格中夜间灯光数据的平均值，得到格网统计下的夜间灯光强度；
服务型POI密度计算：对POI数据进行坐标转换，并与某一尺度下的格网做相交运算，得到每个网格中的POI数据统计结果，每个网格中POI密度采用式(3)计算，
Figure BDA0002521976620000113
式中，POIDensityi为第i个网格中的服务型POI密度，∑Ni为第i个网格中的POI数量，S0为网格面积。
在确定了城市边缘区识别指标因子的基础上，获取同一年份相应的数据源，采用上述计算方法依次进行数据格网化处理、指标权重计算、格网尺度等工作。将城市边缘区识别指标因子量化处理后，便可以通过熵权法确定各项指标权重，计算出每个网格中城市化的特征值，具体包括数据标准化处理和指标权重的确定两部分。数据标准化处理主要是对以上不同来源的数据做无量纲化处理，本发明采用极差法对指标无量纲化处理。对于指标因子权重的确定，采用常见的熵权法对上述获得的指标因子赋予权重。其中，熵权法计算步骤如式(4)、(5)、(6)所示：
Figure BDA0002521976620000121
Figure BDA0002521976620000122
Figure BDA0002521976620000123
式中，Ej为第j项指标因子的熵，pij为第j项指标下第i个网格的值在整个j项指标序列集合中的比重，wj为第j项指标因子的权重；
在数据格网化处理、指标权重计算后，最终得到多项指标融合后的综合值，称为城市化特征值，其函数形式如式(7)所示：
V＝w1*buildDensity+w2*roadDensity+w3*light+w4*POIDensity (7)
式中，V为城市化特征值，buildDensity为建设用地密度，roadDensity为道路密度，light为夜间灯光强度，POIDensity为服务型POI密度，w1、w2、w3和w4分别为各项指标权重。
进一步的，步骤S3中，基于小波分析的突变点检测方法，对得到的城市化特征值进行城市边缘区突变点识别与提取，具体包括：
空间连续小波变换SCWT(Spatial Continuous Wavelet Transform)是异常点检测的有效工具，它通过增大或减小伸缩尺度a来得到信号的高频或低频信息，然后分析信号的近似信息或细节信息，进而分析不同时间(或空间)尺度上的信号局部特征，连续小波变换公式如式(8)所示：
Figure BDA0002521976620000124
式中，S(x)为位置序列信号，
Figure BDA0002521976620000125
为小波，两者内积可以得到小波系数，a为尺度因子，反映小波的周期长度，b为平移因子，反映时间上的平移。
针对城市边缘区的空间突变特征，本实施例选择db1作为小波基函数。
检测主要分为以下几个步骤：
(1)提取不同方向的特征值空间序列曲线。构建过程如图3所示，对样区的四类原数据进行融合，得到样区的城市化特征值图层；然后以城市中心(即城市计算轮廓的特征矩的中心)为起点，以研究区格网外围的网格中心为终点作采样线，得到采样线图层；将点图层与城市化特征值图层相交，可以得到带城市化特征值的点图层，统一按采样线的ID进行分类，相应得到每条线上的数值序列，进而得到城市化特征值空间序列曲线。为了更加清楚地了解城市化特征值空间序列曲线突变情况，选择其中一条曲线进行说明，如图4所示，将其进行一维投影映射，可以得到城市化特征值序列曲线。由图4可知，该曲线的有效采样点个数为46个，该曲线在样点12～29之间发生了显著波动，在样点12之前和样点29之后有较小的波动，其城市化特征值符合理想状态下城乡地域结构中要素城市性的分布情况。
(2)计算小波系数。使用db1小波作为小波基函数，并利用小波方差确定该曲线的小波变换尺度；然后使用确定好的小波基函数和小波变换尺度对该曲线进行突变检测。对图4所示城市化特征值空间序列曲线进行突变检测的结果如图5所示，其中实心圆符号表示极大值，+字形符号表示极小值，都是代表小波系数模极大值的位置，对比图4，可以发现只要是图4中出现波动的地方，都在图5中以模极大值的形式体现出来。对比图4可以发现只要是5中出现波动的地方，都在图6中以模极大值的形式体现出来，但如果简单的将模极大值出现的位置直接作为突变点映射到图层中，就会出现一些“伪突变点”，这些点并不是真正反映城市边缘区边界的点群，因此需要对这些“伪突变点”进一步筛选剔除。本实施例利用2倍标准差的方法对这些突变性不明显或娇弱的突变点进行剔除。如图6中图b所示，两条标准差区间内的突变点即为“伪突变点”，剔除掉后将剩余的突变点映射到二维空间中该采样线所经过的格网中心点，如图7所示，从图中可以看到突变点出现的位置基本是城市化特征值波动较大的区域。
(3)按照上述过程对每一条采样线进行突变检测，提取城市化特征值空间序列曲线上的突变点，小波系数模极大值位置即为突变点的原则获取突变点。
(4)生成城市边缘区突变点图层，将城市化特征值空间序列曲线上的突变点映射为格网中心，并进行去重处理，得到城市边缘区突变点，结果如图6，从图中可以发现每条线上的突变点组成了突变点群，主要围绕着城市中心分布，城市化特征值越紊乱的地方，突变点出现得越多，主要出现在城市边缘区；反之，特征值越稳定，突变点出现得越少，高值稳定区处于城市建成区，低值稳定区处于乡村腹地。
进一步的，本实施例采用一种顾及空间位置和属性特征共同约束作用下的约束型Delaunay三角网聚类方法，Delaunay三角网作为一种空间聚类的方法，是通过一定的约束准则将一组具有相关性的空间实体划分成若干个具有一定意义的空间簇，进而获得空间群体的结构特征。因此首先对突变点群构建Delaunay三角网，对Delaunay三角网进行双阈值的规则约束，得到点群间的空间邻近关系，考虑到点群含有城市化特征值属性信息，故在计算出的空间邻近关系的基础上进一步顾及城市化特征值的属性约束，得到突变点聚类簇。具体方案设计流程如图8所示。具体步骤包括：
S401，构建Delaunay三角网：对城市边缘区突变点群构建Delaunay三角网；
S402，基于空间位置的聚类：对步骤S401构建的三角网进行双阈值规则约束，得到城市边缘区突变点群间的空间邻近关系；
S403，基于属性特征的聚类：在空间邻近关系的基础上，基于属性权重法进行属性聚类，对城市边缘区突变点群间的连接边进一步约束，得到多个聚类簇；
S404，聚类簇边界提取：基于约束聚类簇边界识别提取算法识别聚类簇的内外边界，对每个聚类簇进行内外边界提取，得到城市边缘区边界。
其中，顾忌空间位置约束的聚类中，Delaunay三角网聚类过程中不需要人为设置阈值参数，利用三角网自身的角度和边长提取参数信息，相比需要设置密度阈值、核密度窗口大小的聚类方法而言，Delaunay三角网聚类有利于减少人为误差的产生。本发明采用一种双阈值约束方法对Delaunay三角网的冗余边进行剔除，得到突变点群在空间上的邻近关系。
1、全局冗余边去除
建立突变点群的三角网后，首先从全局的角度删除三角网的冗余边，三角网的边长平均值和标准差能够较好地反映点群分布的整体特征，所以删除这些冗余边后能粗略地反映空间簇的情况。对于任意的点p，与其直接相连的边为ei，全局冗余边的约束规则如式(9)所示：
Figure BDA0002521976620000151
式中，Global_LongEdges(p)为p点处需要删除的全局冗余边的集合，GlobalMean为三角网的所有边的平均长度，PartialMean(p)为与p点直接相连边的平均长度，GlobalSD为三角网的所有边长的标准差，ei为任意的点p与其直接相连的边；
利用该约束规则对模拟数据的全局冗余边进行删除，如图9所示。
2、局部冗余边去除
当全局冗余边去除并进行拓扑检查后，某些局部冗余边仍然存在，如图9中圈出的区域所示。对于该情况，需要设计相应的局部冗余边约束规则，对于任意点，在其2阶邻域内考虑局部因素的影响，即考虑到该点距离等于1或2的其余点。对于任意点p，与其2阶邻域内相连的边为ei，则局部冗余边的约束规则可以定义为如式(10)所示：
Figure BDA0002521976620000152
式中，Partial_LongEdges(p)为p点处需要删除的局部冗余边的集合，PartialMean 2(p)为p点2阶邻域内所有边长的平均值，PartialSD(Pj)为与p点直接相连边的标准差，ei为与p相连的2阶邻域内的边。
删除局部冗余边后，再进行一次拓扑检查，最终结果如图10所示，从中可以看出突变点群的空间邻近域及其聚类簇的分布情况。
进一步的，构建突变点群空间邻近域之后，可以通过属性特征(城市化特征值)进行聚类。目前顾及空间邻近和属性相似的聚类方法主要有三种：1)双重距离聚类方法，该方法将属性距离和空间距离均作为聚类变量，利用二元关系来判断实体间的相似性，即单个实体与单个实体间的关系，而二元相似性应用于复杂数据获取有意义的数据结构时会比较困难，这是因为局部区域内会出现差异性很小的相邻实体，二元相似关系在处理这种情况时，误差会在聚类过程中不断累积变大，造成最后的聚类结果无法反映实体在空间分布中的过渡性。2)融合法，该方法将空间位置与属性信息通过加权融合构建距离函数，再通过特定的空间聚类方法进行聚类。由于空间位置与属性信息无可比性，所以该方法在确定两者权值时存在一定的人为主观性。3)空间位置上的属性相似性聚类，该方法先计算出实体间的空间位置邻近关系，再用属性特征进行约束，但是这种方法设定了固定阈值，此方法虽然揭示了实体的空间分布规律，但是不适用于属性分布不均的情况。本实施例将属性权重引入到空间聚类算法中。属性权值可用于定量化描述数据的紊乱程度，属性权值越小，说明数据越有序，反之，越无序。将属性权重运用到聚类中时，可以发现当同一类簇中属性越相似，数据就会越有序。因此属性权重聚类应当能让聚类后的簇内熵尽量小，簇间的熵尽量大。本发明提出基于属性权重的属性聚类方法，进一步根据突变点的属性(城市化特征值)对聚类结果进行分类，属性权重是表示一阶邻域内中心突变点与周围其它突变点的属性相似性指标，
计算公式如式(11)所示，
Figure BDA0002521976620000161
其中，DAEnei(O)为突变点O的邻域属性权值，Eoc表示突变点O与聚类簇C的属性相似性，聚类簇C由n个点{C1,C2,C3,...,Cn}组成，聚类簇C是中心突变点一阶邻域范围内的突变点集；
该值越大，表示中心点与其它连接点差异越小。
计算出每一个突变点的属性权值后，将属性权值最大的突变点作为起始点，然后利用式(12)和(13)，将起始点作为中心点，将周围每个点视为一个单独的聚类簇C，计算中心点与周围每个点之间的属性相似性Eoc，初始聚类簇由属性权值最大的突变点和周围属性相似度Eoc最大的点组合而成。将初始聚类簇一阶邻域内的突变点设为候选点，继续利用式(12)和(13)来度量每个候选点与初始聚类簇间的属性相似性。式(14)和(15)是使相似度Eoc在不同情况下具有可比性的量化方法：
Figure BDA0002521976620000162
Figure BDA0002521976620000163
式中，Eoc为中心点与周围每个点之间的属性相似性；vi为聚类簇中每个突变点的属性值，vn+1为中心突变点的属性值，
Figure BDA0002521976620000171
Figure BDA0002521976620000172
式中，θ为用于比较相似度Eoc的标准化参数，
Figure BDA0002521976620000173
为突变点O与聚类簇C间的最大属性权值，这是假设突变点的属性值和聚类簇点群的属性值完全相等时计算出的最大属性权值情况。
显然θ∈[0,1]，θ值越大，突变点O与聚类簇之间的相似性越大。我们需要设定一个阈值，如果θ大于阈值，则允许突变点O加入到聚类簇C中，如果不是，则跳过该点；重复上述候选点筛选步骤，继续扩大聚类簇，直至簇的一阶邻域内没有类似的突变点(即θ大于阈值的)，从而得到第一个簇；然后在其余的突变点中将属性权值最大的突变点作为第二簇的起始点。重复上述处理过程，直至将所有突变点分成不同的簇。
θ阈值的设定是本部分的关键技术，它能够影响聚类的结果。本发明引入PBM指数来确定θ的取值，如式(16)和(17)所示，PBM指数是一个计算聚类簇的紧致性和分离性的一个相对评价指标，该值越大，聚类结果越可靠。随着θ的不断变化，我们用迭代法来计算每个聚类簇的PBM指数，当PBM值最大时，得到最优θ值。
Figure BDA0002521976620000174
Figure BDA0002521976620000175
式中，Nc为聚类簇的个数，Ni为聚类簇Ci中的突变点个数，zj为聚类簇Ci的质心坐标，xj为第j个突变点的坐标。
经过上述一系列步骤，对城市化特征值突变点群进行聚类，得到模拟数据的空间聚类结果，如图11所示。
进一步的，本实施例采用一种约束聚类簇识别提取算法对每个聚类簇进行内外边界提取，此算法的优点在于可以构建有空洞对象和无空洞对象两种情况下的点群边界。城市边缘区由内外边界组合而成，从几何角度来看，城市边缘区是一个含空洞的几何空间。该方法由两部分组成，分为轮廓边界提取和孔洞边界提取。
轮廓边界提取是通过迭代删除Delaunay三角网最外层的边界边，判断条件是：边界边对应的角度是否大于阈值R，将大于阈值R的对应边删除，通过轮廓边界提取得到城市边缘区突变点集群的外边界；
轮廓边界处理完后就能得到点集的外边界，然而这种删除处理不能解决空洞情况。当点集存在空洞时，则会存在内边界，此时需要通过孔洞边界提取算法进一步处理。孔洞边界提取是基于轮廓边界提取结果，计算每个聚类簇的松散度和冗余边准则，将每个聚类簇中的突变点根据松散度进行排序，按松散度阈值将每个聚类簇划分为两个子簇，阈值通过PBM指数确定，提取松散度大于阈值的突变点；对每个突变点应用冗余边准则，当突变点p的第i条边大于一定值时，则将该边删除，通过孔洞边界提取得到城市边缘区突变点集群的内边界；
其中，松散度采用式(18)计算，
F(p)＝PartialSD(p)/PartialMean(p) (18)
式中，p为突变点，F(p)为突变点p的松散度，PartialMean(p)为与p相连的边界边长的平均值，PartialSD(p)为与p相连的边界边长的标准差；
冗余边准则采用式(19)计算
T(p)＝{ei|ei＞PartialMean(p)+GlobalSD} (19)
式中，p为突变点，T(p)为突变点p的冗余边准则，PartialMean(p)为与p相连的边界边长的平均值，GlobalSD为所有边的边长标准差，ei为与p相连的第i条边的长度。
当突变点同时满足轮廓准则和孔洞准则时，处理完后的结果即为城市边缘区的内外边界结果。
图12本发明实施例最终实验数据得到的城市边缘区提取结果示意图。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于多源数据融合的城市边缘区提取方法，其特征在于，该方法包括以下步骤：
S1，选取多项城市边缘特征数据作为数据源，并分别确定各项城市边缘特征数据的城市边缘区识别指标因子；
S2，对步骤S1确定的各城市边缘区识别指标因子分别进行数据标准化处理，计算各城市边缘区识别指标因子的指标值和指标权重，并通过上述指标值和指标权重计算得到城市化特征值；
S3，根据步骤S3计算得到的城市化特征值进行城市边缘区突变点识别和提取，得到城市边缘区突变点集群；
S4，通过对步骤S3得到的城市边缘区突变点集群进行聚类处理，并对聚类处理得到的每个聚类簇进行边界提取，得到城市边缘区边界。
2.根据权利要求1所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S1中，所述多项城市边缘特征数据包括：土地利用数据、道路数据、夜间灯光数据和POI数据，所述多项城市边缘特征数据对应的城市边缘区识别指标因子分别为建设用地密度、道路密度、夜间灯光强度和服务型POI密度。
3.根据权利要求2所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S1中，确定建设用地密度、道路密度、夜间灯光强度和服务型POI密度的指标值的具体方式如下，其中，
确定建设用地密度，包括：将所述土地利用数据分为建设用地、非建设用地和水域三大类，并与一格网做相交运算，得到格网中每个网格的土地利用数据信息，并对其进行格网统计，每个网格中的建设用地密度采用如下公式计算，
Figure FDA0002521976610000011
式中，buildDensityi为第i个网格中的建设用地密度，∑Si为第i个网格中的建设用地面积之和，S0为网格面积；
确定道路密度，包括：对所述道路数据进行了预处理，剔除国道、省道和高速公路，保留快速路、主干路、次干路和支路，将预处理后的道路数据与一格网进行相交运算，并分类统计每个网格对应的道路数据，得到道路数据的格网统计结果，每个网格中的道路密度采用如下公式计算，
Figure FDA0002521976610000021
式中，roadDensityi为第i个网格中的道路密度，∑Li为第i个网格中的道路长度之和，S0为网格面积；
确定夜间灯光强度，包括：对所述夜间灯光数据进行重采样预处理，将预处理后的夜间灯光数据与一格网进行相交运算，并对夜间灯光数据进行分区统计，计算每个网格中夜间灯光数据的平均值，得到格网统计下的每个网格的夜间灯光强度；
确定服务型POI密度，包括：对所述POI数据进行坐标转换，并与一格网做相交运算，得到每个网格中的POI数据统计结果，每个网格中POI密度采用如下公式计算，
Figure FDA0002521976610000022
式中，POIDensityi为第i个网格中的服务型POI密度，∑Ni为第i个网格中的POI数量，S0为网格面积。
4.根据权利要求1所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S2中，数据标准化处理采用极差法对所述各城市边缘区识别指标因子进行无量纲化处理，并通过熵权法确定所述各城市边缘区识别指标因子的各指标权重。
5.根据权利要求1所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S3中，基于小波分析的突变点检测方法，对得到的城市化特征值进行城市边缘区突变点识别与提取，具体步骤包括：
S301，提取城市化特征值空间序列曲线：
基于各项数据源，获取城市化特征图层；
以城市中心为起点，城市外围的网格中心为终点作多条采样线，得到采样线图层，并在同方向上的采样线上等距离取点，以使采样线经过的格网中的每个网格中至少有一个采样点；
将采样线图层与城市化特征值图层做相交运算，得到城市化特征值的点图层；
在所述点图层中，确定每条采样线上从所述城市中心至所述网格中心的各采样点对应的城市化特征值，得到城市化特征值空间序列曲线；
S302，确定所述城市化特征值空间序列曲线的小波系数：
通过db1小波基函数对所述城市化特征值空间序列曲线进行多尺度小波变换；
确定所述多尺度小波变换后各采样线所对应的各小波系数模极大值；
S303，提取所述城市化特征值空间序列曲线上的突变点：
对每条采样线，将对应的小波系数模极大值位置作为候选突变点，通过两倍标准差法剔除所述候选突变点中的伪突变点；
将所述候选突变点中剔除所述伪突变点后剩余的突变点映射到二维空间中该采样线所经过的网格中心；
依次类推，提取所述各采样线对应的突变点；
S304，生成城市边缘区突变点图层：
将城市化特征值空间序列曲线上所述各采样线的突变点映射为网格中心；
对映射的网格中心进行去重处理，得到城市边缘区突变点集群。
6.根据权利要求1所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S4包括：
S401，构建Delaunay三角网：对城市边缘区突变点群构建Delaunay三角网；
S402，基于空间位置的聚类：对所述Delaunay三角网进行双阈值法规则约束，得到所述城市边缘区突变点群中各突变点的空间邻近关系；
S403，基于属性特征的聚类：在空间邻近关系的基础上，基于属性权重法进行属性聚类，对所述城市边缘区突变点群中各突变点的连接边进一步约束，得到多个聚类簇；
S404，聚类簇边界提取：基于约束聚类簇边界提取方法识别聚类簇的内外边界，并对每个聚类簇进行内外边界提取，得到城市边缘区边界。
7.根据权利要求6所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤402中，规则约束采用双阈值约束方法对所述Delaunay三角网的冗余边进行剔除，得到城市边缘区突变点群中各突变点间的空间邻近关系，包括：
对所述Delaunay三角网中的全局冗余边和局部冗余边进行剔除，得到城市边缘区突变点群间的空间邻近关系；
其中，
全局冗余边约束规则为：
Figure FDA0002521976610000041
式中，Global_LongEdges(p)为p点处需要删除的全局冗余边的集合，GlobalMean为三角网的所有边的平均长度，PartialMean(p)为与p点直接相连边的平均长度，GlobalSD为三角网的所有边长的标准差，ei为任意的点p与其直接相连的边；
局部冗余边的约束规则为：
Figure FDA0002521976610000042
式中，Partial_LongEdges(p)为p点处需要删除的局部冗余边的集合，PartialMean 2(p)为p点2阶邻域内所有边长的平均值，PartialSD(Pj)为与p点直接相连边的标准差，ei为与p相连的2阶邻域内的边。
8.根据权利要求6所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S403中，基于属性权重法进行属性特征聚类，对每一个城市边缘区突变点进行如下操作步骤：
a1，根据城市边缘区突变点的属性对聚类结果进行分类，并计算每一个城市边缘区突变点的邻域属性权值，邻域属性权值采用如下公式计算，
Figure FDA0002521976610000043
式中，DAEnei(O)为突变点O的邻域属性权值，Eoc表示突变点O与聚类簇C的属性相似性，聚类簇C由n个点{C1,C2,C3,...,Cn}组成，聚类簇C是中心突变点一阶邻域范围内的突变点集；
a2，选取邻域属性权值最大的城市边缘区突变点作为起始突变点，将起始突变点作为中心突变点，将中心突变点周围的每个突变点视为一个单独的聚类簇C，计算中心突变点与周围每个突变点之间的属性相似性Eoc，由属性权值最大的突变点和周围属性相似度Eoc最大的突变点组成初始聚类簇，将初始聚类簇一阶邻域内的突变点设为候选突变点，利用如下公式继续度量每个候选突变点与初始聚类簇间的属性相似性；
Figure FDA0002521976610000051
Figure FDA0002521976610000052
式中，Eoc为中心点与周围每个点之间的属性相似性；vi为聚类簇中每个突变点的属性值，vn+1为中心突变点的属性值；
a3，对属性相似性进行量化，若突变点与聚类簇之间的相似性量化值大与预设阈值，则将该突变点加入到聚类簇C中，否则，跳过该突变点，重复步骤a2，继续扩大聚类簇，直至聚类簇的一阶邻域内没有相似性量化值大于预设阈值的突变点，得到第一个聚类簇；
a4，在其余的突变点中将属性权值最大的突变点作为第二簇的起始点，重复步骤a2和a3，使所有城市边缘区突变点群中的各突变点分成不同的聚类簇。
9.根据权利要求8所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤a3中，属性相似性量化方式为：
Figure FDA0002521976610000053
Figure FDA0002521976610000054
式中，θ为用于比较相似度Eoc的标准化参数，
Figure FDA0002521976610000055
为突变点O与聚类簇C间的最大属性权值，
其中，标准化参数的取值通过PBM指数确定，所述PBM指数采用如下公式确定：
Figure FDA0002521976610000056
Figure FDA0002521976610000057
式中，Nc为聚类簇的个数，Ni为聚类簇Ci中的突变点个数，zj为聚类簇Ci的质心坐标，xj为第j个突变点的坐标；
当PBM值最大时，得到最优θ值。
10.根据权利要求6所述的基于多源数据融合的城市边缘区提取方法，其特征在于，步骤S404中，约束聚类簇边界识别提取方法包括轮廓边界提取和孔洞边界提取，并对每个聚类簇进行内外边界提取，得到城市边缘区边界，包括：
基于轮廓边界提取得到城市边缘区突变点集群的外边界；其中，所述提取所述城市边缘区突变点集群的外边界包括：
判断所述Delaunay三角网最外层的边界边对应的角度是否大于阈值R，若所述边界边对应的角度大于阈值R，则删除所述边界边，依次类推，得到所述城市边缘区突变点集群的外边界；
基于孔洞边界提取得到城市边缘区突变点集群的内边界；其中，基于所述城市边缘区突变点集群的外边界提取所述城市边缘区突变点集群的内边界，包括：
计算每个聚类簇的松散度，将每个聚类簇中的突变点根据松散度进行排序，基于松散度阈值将每个聚类簇划分为两个子簇，提取两个子簇中松散度大于所述松散度阈值的突变点，其中，所述松散度阈值通过PBM指数确定；计算每个聚类簇的各突变点的冗余边准则，若突变点p的第i条边大于预设阈值时，则将该边删除，得到所述城市边缘区突变点集群的内边界；
其中，松散度采用如下公式计算，
F(p)＝PartialSD(p)/PartialMean(p)
式中，p为突变点，F(p)为突变点p的松散度，PartialMean(p)为与p相连的边界边长的平均值，PartialSD(p)为与p相连的边界边长的标准差；
冗余边准则采用如下公式计算
T(p)＝{ei|ei＞PartialMean(p)+GlobalSD}
式中，p为突变点，T(p)为突变点p的冗余边准则，PartialMean(p)为与p相连的边界边长的平均值，GlobalSD为所有边的边长标准差，ei为与p相连的第i条边的长度。
CN202010493552.6A 2020-06-03 2020-06-03 一种基于多源数据融合的城市边缘区提取方法 Active CN111651545B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010493552.6A CN111651545B (zh) 2020-06-03 2020-06-03 一种基于多源数据融合的城市边缘区提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010493552.6A CN111651545B (zh) 2020-06-03 2020-06-03 一种基于多源数据融合的城市边缘区提取方法

Publications (2)

Publication Number Publication Date
CN111651545A true CN111651545A (zh) 2020-09-11
CN111651545B CN111651545B (zh) 2024-02-13

Family

ID=72345003

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010493552.6A Active CN111651545B (zh) 2020-06-03 2020-06-03 一种基于多源数据融合的城市边缘区提取方法

Country Status (1)

Country Link
CN (1) CN111651545B (zh)

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112199566A (zh) * 2020-09-27 2021-01-08 成都房联云码科技有限公司 基于地产大数据的城市更新效果评估方法与系统
CN113159406A (zh) * 2021-04-14 2021-07-23 同济大学 基于多源数据的城镇实体范围识别及划定信息处理方法
CN113660623A (zh) * 2021-08-16 2021-11-16 智慧足迹数据科技有限公司 基于手机信令的边界识别方法及相关装置
CN113790759A (zh) * 2021-09-16 2021-12-14 滁州职业技术学院 一种多传感器融合式建筑测量机器人及其控制方法
CN113988198A (zh) * 2021-11-03 2022-01-28 中国地质大学（武汉） 一种基于地标约束的多尺度城市功能分类方法
CN114266138A (zh) * 2021-11-29 2022-04-01 西南大学 一种利用云端数据进行城市边缘区识别与验证方法
CN114580975A (zh) * 2022-05-05 2022-06-03 成都数联云算科技有限公司 街区活力获得方法及系统及装置及介质
CN115564122A (zh) * 2022-10-14 2023-01-03 中国测绘科学研究院 优化建设阈值选择的城市实体地域划定方法、装置及设备
CN115660306A (zh) * 2022-07-28 2023-01-31 深圳大学 一种城市中心识别方法、终端及存储介质
CN115661670A (zh) * 2022-12-14 2023-01-31 四川师范大学 一种城乡自然景观梯度上过渡区识别方法及装置
CN116486031A (zh) * 2023-05-25 2023-07-25 南京大学 时序遥感数据与地理大数据融合的城市更新区域识别方法
CN116645012A (zh) * 2023-07-27 2023-08-25 河北工业大学 一种城市边缘区空间范围的高精度动态识别方法
CN117217556A (zh) * 2023-09-11 2023-12-12 中国农业科学院农业资源与农业区划研究所 多指标间权衡的生产可能性边界生成及阈值识别方法和系统
CN117237780A (zh) * 2023-11-15 2023-12-15 人工智能与数字经济广东省实验室(深圳) 多维数据特征图构建方法、系统、智能终端和介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103745470A (zh) * 2014-01-08 2014-04-23 兰州交通大学 基于小波的多边形轮廓演化医学ct图像交互式分割方法
CN106570468A (zh) * 2016-10-25 2017-04-19 中国人民解放军空军工程大学 一种重建LiDAR原始点云建筑物轮廓线的方法
CN110533038A (zh) * 2019-09-04 2019-12-03 广州市交通规划研究院 一种基于信息数据的城市活力区和中心城区边界识别的方法
CN110634299A (zh) * 2019-10-25 2019-12-31 福州大学 基于多源轨迹数据的城市交通状态精细划分与识别方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103745470A (zh) * 2014-01-08 2014-04-23 兰州交通大学 基于小波的多边形轮廓演化医学ct图像交互式分割方法
CN106570468A (zh) * 2016-10-25 2017-04-19 中国人民解放军空军工程大学 一种重建LiDAR原始点云建筑物轮廓线的方法
CN110533038A (zh) * 2019-09-04 2019-12-03 广州市交通规划研究院 一种基于信息数据的城市活力区和中心城区边界识别的方法
CN110634299A (zh) * 2019-10-25 2019-12-31 福州大学 基于多源轨迹数据的城市交通状态精细划分与识别方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张成虎等: "基于贝叶斯分类的可疑金融交易识别研究" *
王磊等: "基于小波变换和Teager能量算子指纹图像增强" *

Cited By (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112199566B (zh) * 2020-09-27 2021-07-30 成都房联云码科技有限公司 基于地产大数据的城市更新效果评估方法与系统
CN112199566A (zh) * 2020-09-27 2021-01-08 成都房联云码科技有限公司 基于地产大数据的城市更新效果评估方法与系统
CN113159406A (zh) * 2021-04-14 2021-07-23 同济大学 基于多源数据的城镇实体范围识别及划定信息处理方法
CN113660623A (zh) * 2021-08-16 2021-11-16 智慧足迹数据科技有限公司 基于手机信令的边界识别方法及相关装置
CN113790759B (zh) * 2021-09-16 2023-09-29 滁州职业技术学院 一种多传感器融合式建筑测量机器人及其控制方法
CN113790759A (zh) * 2021-09-16 2021-12-14 滁州职业技术学院 一种多传感器融合式建筑测量机器人及其控制方法
CN113988198A (zh) * 2021-11-03 2022-01-28 中国地质大学（武汉） 一种基于地标约束的多尺度城市功能分类方法
CN113988198B (zh) * 2021-11-03 2024-05-14 中国地质大学（武汉） 一种基于地标约束的多尺度城市功能分类方法
CN114266138A (zh) * 2021-11-29 2022-04-01 西南大学 一种利用云端数据进行城市边缘区识别与验证方法
CN114580975A (zh) * 2022-05-05 2022-06-03 成都数联云算科技有限公司 街区活力获得方法及系统及装置及介质
CN115660306A (zh) * 2022-07-28 2023-01-31 深圳大学 一种城市中心识别方法、终端及存储介质
CN115564122A (zh) * 2022-10-14 2023-01-03 中国测绘科学研究院 优化建设阈值选择的城市实体地域划定方法、装置及设备
CN115661670A (zh) * 2022-12-14 2023-01-31 四川师范大学 一种城乡自然景观梯度上过渡区识别方法及装置
CN115661670B (zh) * 2022-12-14 2023-04-07 四川师范大学 一种城乡自然景观梯度上过渡区识别方法及装置
CN116486031A (zh) * 2023-05-25 2023-07-25 南京大学 时序遥感数据与地理大数据融合的城市更新区域识别方法
CN116645012A (zh) * 2023-07-27 2023-08-25 河北工业大学 一种城市边缘区空间范围的高精度动态识别方法
CN116645012B (zh) * 2023-07-27 2023-10-10 河北工业大学 一种城市边缘区空间范围的高精度动态识别方法
CN117217556A (zh) * 2023-09-11 2023-12-12 中国农业科学院农业资源与农业区划研究所 多指标间权衡的生产可能性边界生成及阈值识别方法和系统
CN117217556B (zh) * 2023-09-11 2024-05-14 中国农业科学院农业资源与农业区划研究所 多指标间权衡的生产可能性边界生成及阈值识别方法和系统
CN117237780A (zh) * 2023-11-15 2023-12-15 人工智能与数字经济广东省实验室(深圳) 多维数据特征图构建方法、系统、智能终端和介质
CN117237780B (zh) * 2023-11-15 2024-03-15 人工智能与数字经济广东省实验室(深圳) 多维数据特征图构建方法、系统、智能终端和介质

Also Published As

Publication number Publication date
CN111651545B (zh) 2024-02-13

Similar Documents

Publication Publication Date Title
CN111651545A (zh) 一种基于多源数据融合的城市边缘区提取方法
CN112949413B (zh) 基于街景图片的城市景观要素分类及地方性度量方法
CN111652198B (zh) 一种城市边缘区识别方法及系统
CN109493119B (zh) 一种基于poi数据的城市商业中心识别方法及系统
CN111489008B (zh) 一种交通事故影响因素空间效应分析方法及其应用
CN110929939B (zh) 一种基于聚类-信息量耦合模型下的滑坡灾害易发性空间预测方法
CN111985389B (zh) 一种基于流域属性距离的流域相似判别方法
CN111814596B (zh) 一种遥感影像与出租车轨迹融合的自动城市功能分区方法
CN111651502B (zh) 一种基于多子空间模型的城市功能区识别方法
CN111797188B (zh) 一种基于开源地理空间矢量数据的城市功能区定量识别方法
Haining Spatial statistics and the analysis of health data
CN111984701A (zh) 乡村聚落演化的预测方法、装置、设备及存储介质
CN112365595A (zh) 基于报警数据点要素的鉴别交通事故多发点的分析方法
CN114997534A (zh) 基于视觉特征的相似降雨预报方法和设备
CN114565207A (zh) 一种融合属性数据和流动数据的城市群高质量发展监测评估方法
CN113379269A (zh) 多因素空间聚类的城市商业功能区划方法、装置及介质
CN116110210B (zh) 复杂环境下数据驱动的滑坡灾害辅助决策方法
CN115186074B (zh) 基于Meta分析的模拟土壤pH值空间分布格局的方法
Zhou et al. Mapping urban landscape heterogeneity: agreement between visual interpretation and digital classification approaches
CN116013084A (zh) 交通管控场景确定方法、装置、电子设备及存储介质
CN111710157B (zh) 一种出租车热点区域的提取方法
CN115879594A (zh) 一种基于地理探测器的城市定居人口分布趋势预测方法
Liu et al. Road density analysis based on skeleton partitioning for road generalization
CN112925784B (zh) 一种实有人口数据多尺度空间化方法
CN109190783A (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