CN104318066A

CN104318066A - 一种自然地表特征表征方法

Info

Publication number: CN104318066A
Application number: CN201410508901.1A
Authority: CN
Inventors: 刘耀林; 刘艳芳; 焦利民; 何建华
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2014-09-26
Filing date: 2014-09-26
Publication date: 2015-01-28

Abstract

本发明公开了一种自然地表特征表征方法，包括自然地形评价、自然地貌评价、自然地表覆盖评价、自然地理国情综合评价、自然地表格局特征表征；本发明使用综合评估模型对区域地理国情要素中的自然地表要素及其对区域生产、生活和生态环境等活动的影响进行综合定量评估。同时使用景观格局分析方法，设计了包括斑块密度、多样性指数等面向自然地表要素空间格局特征表征的系列指数。克服了现有区域自然地表特征计算与表征技术在表征区域自然地表特征时缺乏综合评估模型、无法表征区域自然地表空间特征等不足。本发明提出的技术方案将服务于即将在全国范围内开展地理国情统计分析工作，同时也是面向地理国情时空数据挖掘的技术基础。

Description

一种自然地表特征表征方法

技术领域

本发明属于地理国情统计分析与时空数据挖掘技术领域，特别是涉及一种自然地表特征表征方法。

背景技术

地理国情是地形地貌、地表覆盖、交通网络、城镇布局等多种自然和人文要素在宏观层面的综合表达，反映了国家自然资源、经济社会和人文在地理空间上的时空特征与基本规律，是基本国情的重要组成部分。

地理国情根据其表达的专题要素内容的不同，可分为自然地表、地理区位和生态覆被等多个专题内容以反映地理国情的不同方面。其中，自然地表主要反映区域地形、地貌、地表覆盖和地表格局等自然地理国情要素的时空格局及变化特征；地理区位则主要从区位的角度主要反映区域的自然资源环境、城镇布局和交通网络等区位特征的时空格局及变化特征；生态覆被则主要从生态环境的角度对区域的植被覆盖、水域湿地、生物丰度和土地退化等特征进行评价和表征。

为全面掌握我国地理国情，满足社会经济和生态文明发展的决策需求，国务院提出了于2013年至2015年在全国范围内开展第一次地理国情普查。地理国情信息的普查、监测和统计分析工作的开展对地理空间信息的采集、存储、可视化与统计分析提出了更高要求。遥感技术、地理信息技术、全球定位技术、并行计算技术和时空数据挖掘等技术被广泛用于地理国情要素信息的采集、存储和统计分析。通过开展地理国情普查和监测工作，能够准确、及时的掌握区域地理国情要素(耕地、园地、公路、河流、湖泊等)的空间位置、数量和形状等信息。地理国情普查与监测是地理国情评价、表征与应用的前提，为地理国情的表征和决策应用提供重要的数据基础。

随着全国第一次地理国情普查试点工作的完成和后续工作的推进，如何基于地理国情普查和监测基础数据，结合必要的社会经济统计数据，准确、客观、合理的表征区域地理国情，服务政府、部门和公众的多层次决策，已成为地理国情普查监测成果应用的迫切需求。

基于我国地理普查监测工作的实践，地理国情评价与表征相关的技术主要可分为以下几种：

(一)基于点、线、面的地理国情矢量图形表征技术

地理国情普查监测数据集中，地理国情相关的要素对象被抽象为点、线、面等几何形态进行表征。如，高速公路出入口、学校、医院等被抽象为带属性的点状要素；公路、铁路、河渠被抽象为带有属性的线状要素；湖泊、耕地、园地、城市建成区被抽象为带属性信息的面要素。以点、线、面表征地理国情要素能够非常详细、具体和直观的表达地理国情要素的空间位置、形状。其缺点主要在于缺乏对区域地理国情信息的宏观整体把握，如无法直接表达区域内各类地理国情要素的面积、长度、个数、密度等信息。

(二)基于汇总统计方法的地理国情评价与表征技术

汇总统计方法是以地理国情普查/监测数据为统计对象，针对不同的统计单元，根据地理国情实体的空间形态特征，将其划分为点状地理国情实体、线状地理国情实体和面状地理国情实体，以数量、密度、长度、面积、距离等作为统计指标，开展统计分析。以行政区划为单元，以地理国情要素为统计对象进行的地理国情汇总统计方法和步骤具体见国务院第一次全国地理国情普查领导小组办公室发布的系列技术规范。

基于汇总统计方法的地理国情表征技术一定程度能够从整体上反映区域地理国情要素的数量特征，如获得某个行政区范围内的耕地面积、园地面积、各等级道路的长度、学校设施的密度等。然而，地理国情是一种与空间位置相关的基本国情，以汇总统计的方法表征地理国情只能从数量上反映区域某一类地理国情要素的情况，无法实现对区域地理国情的综合评价，也不能反映区域地理国情的空间格局特征。

自然地表特征是区域地形、地貌、地表覆盖和地表格局特征的综合体现，反映区域的基本自然地理客观条件和地表空间格局。自然地表特征是自然条件和长期人类活动共同作用的结果，也是未来人类社会经济发展物质基础和背景条件，是基本的地理国情。设计一种全面、客观、综合的刻画区域自然地表特征、反映区域自然地理条件和地表空间格局的评价与表征技术，进而服务于社会经济发展决策，是当前地理国情普查与监测中亟待解决的关键技术问题。

现有技术方案无法对区域自然地表特征和格局进行综合、全面的描述和表征。

发明内容

为了解决上述的技术问题，本发明面向区域自然地表特征综合表达与应用需求，提出一种自然地表特征表征方法，为地理国情普查监测的后续成果统计分析、数据挖掘与应用工作奠定技术基础，推进地理国情普查监测成果服务于区域发展决策水平。

本发明所采用的技术方案是：一种自然地表特征表征方法，其特征在于，包括以下步骤：

步骤1：确定要进行自然地表特征评价与表征的空间范围和统计单元；

步骤2：从地理国情普查或监测数据集中读取计算范围内地理国情普查的10米空间分辨率数字高程模型(DEM)数据；其中，DEM以栅格像元的方式进行存储，每个栅格像元表示10米*10米的一个空间位置，栅格像元的灰度值表示了该位置的高程(海拔)；遍历计算范围内的所有统计单元，依次计算各统计单元的地形指数(TI)；其中，地形指数用于综合刻画和反映统计单元内的自然地形条件和特征；

步骤3：遍历计算范围内的所有统计单元，依次计算各统计单元的地貌指数(FI)；

步骤4：遍历范围内的所有统计单元，依次计算各统计单元的地表覆盖指数(CI)；

步骤5：在步骤2、步骤3、步骤4计算结果的基础上，遍历计算范围内的所有统计单元，计算所有统计单元的自然地表指数(LI)以综合评估区域自然地理国情综合条件对于生产、生活和生态环境保护的优势；

步骤6：遍历区域内的所有统计单元，从斑块格局、类型格局和区域总体格局三个不同层次，依次提取各统计单元的自然地理国情要素的空间格局分布特征。

作为优选，步骤1中所述的空间范围是指县及县级以上的行政区或特定的流域、城市圈、城市群、经济带类在地域空间上连片的、封闭的一个区域；所述的统计单元包括标准网格单元和行政区划单元，是地理国情评价与表征的基本单位。

作为优选，步骤2中所述的地形指数(TI)的计算方法包括以下子步骤：

步骤2.1：利用地理信息技术(GIS)，将各个统计单元和区域DEM空间叠置，获取落在各统计单元多边形内的DEM像元集；在此基础上统计出每个统计单元内栅格像元集的最大值Emax，最小值Emin，平均值Eavg；

步骤2.2：利用地理信息技术(GIS)，将各个统计单元和区域DEM叠置，计算各统计单元的平面投影面积Ap及其在DEM上的表面面积As；

步骤2.3：计算统计单元内的地形起伏度TA，宏观反映统计单元内的自然地形高低起伏程度，其中地形起伏度的计算公式为:TA＝Emax-Emin；

步骤2.4：计算统计单元内的地表粗糙度TR，宏观反映统计单元内的自然地表粗糙程度，其中地表粗糙度的计算公式为：TR＝As/Ap；

步骤2.5：计算统计单元内的地表切割深度TP，直观反映统计单元内地表被侵蚀切割的情况，并对这一地学现象进行了量化，其计算公式为TP＝Eavg-Emin；

步骤2.6：对区域内各个单元的Eavg、TA、TR和TP进行标准化，使各项指标的值的取值区间在1-100之间；各项指标的标准化的方法分别为：

f_{u} = (1 - \frac{t_{u}}{T}) \times 100;

其中：f_u为统计单元u标准化后的指标值，t_u为统计单元u标准化前的原始指标值，T为该项指标在统计范围内各单元的最大值；

步骤2.7：计算各统计单元内的地形指数TI，综合反映统计单元内的自然地形特征；对于任意统计单元u的地形指数TI计算公式为：

TI＝W1×f_ele+W2×f_TA+W3×f_TR+W4×T_TP；

其中W1～W4分别对应平均高程、地形起伏度、地表粗糙度和地表切割深度指标的权重，取值范围为(0-100)；f_ele为经过标准化后的统计单元u的平均高程指标值；f_TA为经过标准化后的统计单元u的地形起伏度指标值；f_TR为经过标准化后的统计单元u的地表粗糙度指标值；f_TP为经过标准化后的统计单元u的地表切割深度指标值。

作为优选，步骤2.7中所述的权重采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。

作为优选，步骤2.7中所述的权重：W1取10，W2取30，W3取40，W4取20。

作为优选，步骤3中所述的地貌指数计算公式如下：

FI＝(1×f_平原+0.8×f_盆地+0.6×f_丘陵-f_高原-f_山地+1)×50；

其中，f_平原、f_盆地、f_丘陵、f_高原、f_山地分别表示各种地貌在统计单元内所占的面积比例，取值为0-100％之间；上述面积比例值直接来自于地理国情普查工作的结果数据，计算方法为：f_i＝a_i/A；其中，A指统计单元的面积，a_i指该统计单元内某种地貌类型的面积，f_i即为该地貌类型的面积比例。

作为优选，步骤4中所述的地表覆盖指数计算公式如下：

CI＝(f_耕+f_园+0.5×(f_林+f_草+f_水)-f_荒+1)×50；

其中，f_耕、f_园、f_林、f_草、f_水、f_荒分别表示各种耕地、园地、林地、草地、水域和荒地在统计单元内所占的面积比例，取值为0-100％之间；上述面积比例值直接来自于地理国情普查工作的结果数据，计算方法为：f_i＝a_i/A；其中，A指统计单元的面积，a_i指该统计单元内某种地表覆盖类型的面积，f_i即为该地表覆盖类型的面积比例。

作为优选，步骤5中所述的自然地表指数计算公式如下：

LI＝E_t×TI+W_f×FI+W_c×CI；

其中W_t、W_f和W_c分别对应地形指数、地貌指数和地表覆盖的权重，取值范围为(0-100)。

作为优选，步骤5中所述的权重采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。

作为优选，步骤5中所述的权重：W_t取30，W_f取30，W_c取40。

作为优选，步骤6中所述的区域总体格局表征包括地表斑块密度、地表边缘密度、地表多样性指数、地表优势度指数、地表均匀度指数和地表破碎度指数共6个定量指标；

所述的地表斑块密度等于统计单元内地表覆盖图层中所有的斑块总数除以统计单元的总面积，反映统计单元内地理国情要素斑块在空间上的分布密度；其计算公式为：

PD = \frac{1}{A} Σ_{i = 1}^{M} N_{i};

式中:PD表示统计单元内的地表斑块密度，M表示统计单元内的地表覆盖类型的个数，N_i表示统计单元内第i种地表覆盖类型的斑块个数(在景观格局分析和空间分析中，斑块可以理解为图斑或地块，即用一个封闭多边形表示空间上连片的、具有相同地表覆盖类型的区域，对应地理国情普查数据集中地表覆盖要素图层中的一个图斑)，A表示统计单元的地表总面积；

所述的地表边缘密度反映统计单元内不同类型地理国情要素斑块间的边界长度，其计算公式为：

ED = \frac{1}{A} Σ_{i = 1}^{M} Σ_{j = 1}^{M} P_{ij};

其中，P_ij表示景观中第i类地理国情要素斑块与相邻第j类地理国情要素斑块间的边界长度；

所述的地表多样性指数是基于信息论基础上，用来度量一定范围内地理国情要素系统的空间结构组成复杂程度的指数，其计算公式为：

H = - Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

式中，H为多样性指数，m为地表覆盖类型总数；p_k为地表覆盖类型面积比例，即某类地表覆盖类型斑块的总面积与统计单元内所有地表覆盖类型斑块总面积的比值；当m＝1，说明一个地区只有一种地表覆盖类型，此时H＝0，多样性指数最小；如果某一个地区地表覆盖类型越多，则m越大；当m→∞时，多样性指数H达到最大值；

所述的地表优势度指数用于测定地表空间格局构成中一种或几种类型支配整个地表构成的程度，计算公式为：

D = H_{\max} + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k}) = {Iog}_{2} (m) + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

式中，D为优势度指数；H是多样性指数，H_max是其最大值，是给定区域内地表覆盖最大可能多样性指数，p_k为地表覆盖类型面积比例，m为地表覆盖类型总数；

所述的地表均匀度指数用以反映地表覆盖空间格局构成中各类型分布的均匀程度，均匀度可表达为：

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)};

所述的地表破碎度指数反映地表覆盖景观的破碎化程度，即地表覆盖被分割的破碎程度，计算公式为：

C＝N/A；

式中，N为区域地表覆盖图斑总数，A为整个区域土地总面积。

作为优选，步骤6中所述的类型格局表征包括类型斑块密度、类型边缘密度、类型多样性指数、类型优势度指数、类型均匀度指数、类型破碎度指数共6个定量指标；

所述的类型斑块密度等于类型斑块总数除以类型总面积，用于反映单元内某种地表覆盖类型的斑块分布状态，计算公式为：

{PD}_{i} = \frac{N_{i}}{A_{i}};

式中：PD_i表示类型斑块密度，N_i表示研究范围内类型斑块总数，A_i表示研究范围某类型总面积；

所述的类型边缘密度为统计单元内单位面积上不同类型景观要素斑块间的边缘长度，计算公式为：

{ED}_{i} = \frac{1}{A_{i}} Σ_{j = 1}^{M} P_{ij};

式中：P_ij表示景观中第i类景观要素斑块与相邻第j类景观要素斑块间的边界长度；

所述的类型多样性指数是基于信息论基础上，用来度量类型结构组成复杂程度的指数，计算公式为：

H = - Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

式中，H为多样性指数，m为某类型图斑总数，p_k为斑块面积比例，即某类地表覆盖类型斑块面积与该类型斑块总面积的比值；

所述的类型优势度指数用于测定某类型图斑构成中一个或几个图斑支配整个类型构成的程度，计算公式为：

D_{t} = H_{t \max} + Σ_{i = 1}^{n} P_{i} {Iog}_{2} (P_{i}) = {Iog}_{2} (n) + Σ_{i = 1}^{n} P_{i} {Iog}_{2} (P_{i});

式中，D_t为类型t的优势度指数；n为类型t的图斑个数；H_max是Shannon多样性指数最大值，P_i为某类型图斑i占该地类总面积的比例；

所述的类型均匀度指数用以反映某类型构成中各图斑大小分布的均匀程度，均匀度可表达为：

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)};

所述的类型破碎度指数反映地表覆盖景观的破碎化程度，即某类地表覆盖被分割的破碎程度，计算公式为：

C＝N_i/A_i；

式中，N_i为区域某类图斑总数，A_i为某类地表总面积。

作为优选，步骤6中所述的斑块格局表征包括斑块平均面积、最大和最小斑块面积、斑块面积标准差、斑块面积变动系数、斑块圆度、斑块矩度、斑块形状复杂度共7个定量指标；

所述的斑块平均面积表示斑块平均大小，计算公式为：

\overset{&OverBar;}{A_{i}} = \frac{1}{N_{i}} Σ_{j = 1}^{N_{i}} A_{ij};

式中，N_i表示第i类斑块总数；A_ij表示第i类第i个斑块的面积；

所述的最大和最小斑块面积表示统计单元内斑块面积的最大最小值；

所述的斑块面积标准差表示统计单元内某类斑块面积的统计标准差；

所述的斑块面积变动系数表示统计单元内斑块规模的变异程度，计算公式为：

C_{i} = \frac{S_{i}}{A_{i}} \times 100 %;

所述的斑块圆度表示斑块与圆的偏离程度，计算公式为：

D = P / 2 \sqrt{π 4};

所述的斑块矩度表示斑块与正方形的偏离程度，计算公式为：

D = 0.25 P / \sqrt{A};

所述的斑块形状复杂度采用斑块分维数表示斑块的形状复杂程度，计算公式为：

P = {kA}^{D / 2}, D = 2 \ln (\frac{P}{k}) / \ln (A);

式中，P是斑块的周长，A是斑块的面积，D是分维数，k是常数。

相对于已有的区域自然地表特征表征技术，本发明使用综合评估模型对区域地理国情要素中的自然地表要素及其对区域生产、生活和生态环境等活动的影响进行综合定量评估。同时使用景观格局分析方法，设计了包括斑块密度、多样性指数等面向自然地表要素空间格局特征表征的系列指数。克服了现有区域自然地表特征计算与表征技术在表征区域自然地表特征时缺乏综合评估模型、无法表征区域自然地表空间特征等不足。本发明提出的技术方案将服务于即将在全国范围内开展地理国情统计分析工作，同时也是面向地理国情时空数据挖掘的技术基础。

附图说明

图1：本发明实施例的流程图。

图2：本发明实施例的统计单元的高程信息提取原理图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明的理论基础为：

(1)多因素综合评价模型

多因素综合评价模型是通过系统、综合地的分析各类因素和因子对区域自然地表特征的影响和作用强度，进而基于定量的方法对评价区域的自然地表特征进行系统、综合的评价，反映区域自然地表对土地利用、经济发展、生态环境保护等人类生产生活活动的难易程度和优劣度。多因素综合评价模型的指导思想是从影响区域自然地表特征或人类活动难易程度的原因出发，采取由原因到结果的评价方法。采用多因素综合评价模型，综合系统的评价区域自然地表特征的关键步骤在于如何科学、合理的确定评价指标体系。

(一)、自然地表特征评价指标体系设计原则

基于多因素综合评价模型的基本原理，构建面向区域自然地表特征评价的自然地表特征评价指标体系的基本原则如下：

1.以地理国情普查和监测数据为基础。区域自然地表特征评价是地理国情普查和监测工作的重要组成部分。地理国情普查和监测是区域自然地表特征评价的主要数据来源。因此，区域自然地表特征评价所设计的评价指标的数据来源也必须以地理国情普查和监测数据为基础，保证设计的指标可获取，具有可行性。

2.具有时空动态变化特征。区域自然地表特征同时具有空间特征和时态特征，表现出时空动态特性。区域自然地表特征在年度或若干年度内具有可监测的变化，通过时序分析可反映区域自然地表特征的时空动态变化状况和变化趋势。

3.综合分析和主导因素分析相结合。自然地表特征具有高度综合性。自然地表特征的评价需要考虑多项指标进行综合计算。但是需要注意构成具体地理国情的指标可能不是同等重要的，在计算时需要考虑因素重要性进行选取或合理确定权重。

(二)、自然地表特征评价指标体系设计方案

基于上述原则，构建区域自然地表特征评价指标体系方案如表1；

表1 自然地表特征综合评价指标体系

(2)景观格局分析方法

景观是在气候、地形、地貌、土壤、植被、水文等自然要素及人文干扰作用下形成的综合有机体。景观格局，一般指景观的空间格局，由一系列大小和形状各异，排列不同的景观斑块在景观空间的排列。景观格局分析是景观生态学的重要组成部分，通过一系列景观格局指数反映区域景观要素在空间上的异质性。由此可见，景观格局分析方法是用于表征自然地表要素在空间上的分布特征和规律。

一、景观格局分析的对象

①斑块。斑块是景观格局的基本组成单元，是指不同于周围背景的、相对均质区域。斑块反映了系统内部和系统间的相似性或相异性。斑块的大小、形状、边界以及斑块之间的距离等空间分布特征构成了不同的生态带，形成了生态系统的差异，对区域生态系统过程产生重要影响。地理国情普查与监测中，斑块对应地表覆盖要素图层中的某一个具体的多边形图形，如某一具体的耕地图斑、园地图斑等。

②类型。景观格局中要素类型相同的所有斑块的集合。例如，在地理国情普查监测中所指的类型包括各类地表覆盖中的各类地理要素的集合，如耕地、园地、林地、草地、水体、裸地等要素。针对某一要素类型的统计即是针对某一种地理国情要素的所有斑块的统计。

③景观格局。区域内由大小、形状、类型各不相同的斑块构成了区域整体的景观格局。对于自然地表综合评价与表征问题，针对景观格局整体的计算就是对评价单元内所有斑块整体所呈现出的特征进行表征。

二、景观格局分析的方法与指标

①斑块密度。反映一定范围内地理国情要素斑块在空间上的分布密度。计算公式为：

PD = \frac{1}{A} Σ_{i = 1}^{M} N_{i};

式中：PD为统计单元内的地表斑块密度，M为统计单元内的地表覆盖类型个数，N_i为统计单元内的地表覆盖类型i的斑块个数(在景观格局分析和空间分析中，斑块可以理解为图斑或地块，即用一个封闭多边形表示空间上连片的、具有相同地表覆盖类型的区域)，A为统计单元内的地表总面积。

②边缘密度。边缘密度反映一定范围内单位面积上不同类型地理国情要素斑块间的边界长度。计算公式为：

式中：P_ij为景观中第i类地理国情要素斑块与相邻第j类地理国情要素斑块间的边界长度。

③多样性指数。基于信息论基础上，用来度量一定范围内地理国情要素系统的空间结构组成复杂程度的指数。计算公式为：

式中，H为多样性指数，m为地表覆盖类型总数。p_k为地表覆盖类型面积比例，即某类地表覆盖类型斑块的总面积与整个区域所有地表覆盖类型斑块总面积的比值。当m＝1，说明一个地区只有一种地表覆盖类型，此时H＝0，多样性指数最小。如果某一个地区地表覆盖类型越多，则m越大。当m→∞时，多样性指数H达到最大值。

④优势度指数。优势度指数用于测定地表空间格局构成中一种或几种类型支配整个地表构成的程度，计算公式为：

D = H_{\max} + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k}) = {Iog}_{2} (m) + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

式中，D为优势度指数；H是多样性指数，H_max是其最大值，是给定区域内地表覆盖最大可能多样性指数，p_k为地表覆盖类型面积比例，m为地表覆盖类型总数。

⑤均匀度指数。用以反映地表覆盖空间格局构成中各类型分布的均匀程度，均匀度可表达为：

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)} .

⑥破碎度指数。地表破碎度指数反映地表覆盖景观的破碎化程度，即地表覆盖被分割的破碎程度。公式为：C＝N/A；其中，N为区域地表覆盖图斑总数，A为整个区域土地总面积。

基于自然地表特征的基本内涵和地理国情普查监测工作的数据特点，设计自然地表的景观格局指标体系如下表2；

表2 自然地表空间特征表征指标体系

如上表所示，分别对统计单元从斑块、类型、景观整体三个不同尺度分别反映统计单元内自然地表的空间特征。

基于上述理论基础，请见图1、图2，本发明所采用的技术方案是：一种自然地表特征表征方法，具体包括以下步骤：

一、确定要进行自然地表特征评价与表征的空间范围和统计单元。其中，所指范围一般是指县及县级以上的行政区或特定的流域、城市圈、城市群、经济带等在地域空间上连片的、封闭的一个区域；统计单元通常包括标准网格单元和行政区划单元等，是地理国情评价与表征的基本单位。计算范围和统计单元的类型、划分方法可参见技术规范《地理国情普查基本统计技术规定(GDPJ 02-2013)》。设统计范围内共有N个统计单元。为了便于理解，本方案描述的范围可理解为县级行政区，统计单元可理解为乡镇级行政区，但具体实施时不限于以上内容，可按照《地理国情普查基本统计技术规定(GDPJ 02-2013)》确定。

二、从地理国情普查或监测数据集中读取计算范围内地理国情普查的10米空间分辨率数字高程模型(DEM)数据。其中，DEM以栅格像元的方式进行存储。每个栅格像元表示10米*10米的一个空间位置，栅格像元的灰度值表示了该位置的高程(海拔)。遍历计算范围内的所有统计单元，依次计算各统计单元的地形指数(TI)。其中，地形指数用于综合刻画和反映统计单元内的自然地形条件和特征。请见图2，地形指数(TI)的计算方法包括以下子步骤：

①利用地理信息技术(GIS)，将各个统计单元和区域DEM空间叠置，获取落在各统计单元多边形内的DEM像元集。在此基础上统计出每个统计单元内栅格像元集的最大值Emax，最小值Emin，平均值Eavg。

②利用地理信息技术(GIS)，将各个统计单元和区域DEM叠置，计算各统计单元的平面投影面积Ap及其在DEM上的表面面积As。多边形的投影面积和表面面积计算方法可利用现有GIS软件实现，非本发明的关注重点，具体计算方法不做阐述。

③计算统计单元内的地形起伏度TA，宏观反映统计单元内的自然地形高低起伏程度。其中地形起伏度的计算公式为:TA＝Emax-Emin。

④计算统计单元内的地表粗糙度TR，宏观反映统计单元内的自然地表粗糙程度。其中地表粗糙度的计算公式为：TR＝As/Ap。

⑤计算统计单元内的地表切割深度TP，直观反映统计单元内地表被侵蚀切割的情况，并对这一地学现象进行了量化，其计算公式为TP＝Eavg-Emin。

⑥对区域内各个单元的Eavg,TA,TR和TP进行标准化，使各项指标的值的取值区间在1-100之间。各项指标的标准化的方法分别为：

f_{u} = (1 - \frac{t_{u}}{T}) \times 100

其中：f_u为统计单元u标准化后的指标值，t_u为统计单元u标准化前的原始指标值，T为该项指标在统计范围内各单元的最大值。

⑦计算各统计单元内的地形指数TI,综合反映统计单元内的自然地形特征。对于任意统计单元u的地形指数TI计算公式为：

TI＝W1×f_ele+W2×f_TA+W3×f_TR+W4×f_TP

其中W1-W4分别对应平均高程、地形起伏度、地表粗糙度和地表切割深度指标的权重，取值范围为(0-100)。f_ele为经过标准化后的统计单元u的平均高程指标值；f_TA为经过标准化后的统计单元u的地形起伏度指标值；f_TR为经过标准化后的统计单元u的地表粗糙度指标值；f_TP为经过标准化后的统计单元u的地表切割深度指标值；

权重一般可采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。也可采用本发明提供的参考权重值：W1取10，W2取30，W3取40，W4取20.

三、遍历计算范围内的所有统计单元，依次计算各统计单元的地貌指数(FI)。地貌指数综合反映了统计单元内的地貌类型条件对生产生活等人类活动的便利度。FI的值越高，则表明该单元内的地貌条件越有利于人类活动的开展，对该区域实施开发利用的成本越低。地貌指数计算公式如下：

FI＝(1×f_平原+0.8×f_盆地+0.6×f_丘陵-f_高原-f_山地+1)×50

其中，f_平原，f_盆地，f_丘陵，f_高原，f_山地分别表示各种地貌在统计单元内所占的面积比例，取值为0-100％之间。上述面积比例值直接来自于地理国情普查工作的结果数据，计算方法为：f_i＝a_i/A.其中，A指统计单元的面积，a_i指该统计单元内某种地貌类型的面积，f_i即为该地貌类型的面积比例。

四、遍历范围内的所有统计单元，依次计算各统计单元的地表覆盖指数(CI)。地表覆盖指数综合反映了统计单元内的各种地表覆盖条件对生产生活和生态环境保护的有利程度。FI的值越高，则表明该单元内的地表覆盖条件越有利于人类活动的开展或越有利于生态环境保护。地表覆盖指数计算公式如下：

CI＝(f_耕+f_园+0.5×(f_林+f_草+f_水)-f_荒+1)×50

其中，f_耕，f_园，f_林，f_草，f_水，f_荒分别表示各种耕地、园地、林地、草地、水域和荒地在统计单元内所占的面积比例，取值为0-100％之间。上述面积比例值直接来自于地理国情普查工作的结果数据，计算方法为：f_i＝a_i/A.其中，A指统计单元的面积，a_i指该统计单元内某种地表覆盖类型的面积，f_i即为该地表覆盖类型的面积比例。

五、在二、三、四计算结果的基础上，遍历计算范围内的所有统计单元，计算所有统计单元的自然地表指数(LI)以综合评估区域自然地理国情综合条件对于生产、生活和生态环境保护的优势。计算公式为：

LI＝W_t×TI+W_f×FI+W_c×CI

其中W_t、W_f和W_c分别对应地形指数、地貌指数和地表覆盖的权重，取值范围为(0-100)。权重一般可采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。也可采用本发明提供的参考权重值：W_t取30，W_f取30，W_c取40。

六、遍历区域内的所有统计单元，从斑块格局、类型格局和区域总体格局三个不同层次，依次提取各统计单元的自然地理国情要素的空间格局分布特征。地表格局特征用于反映统计区域地理国情要素在空间上的分布特征，分别用地表斑块密度、地表边缘密度、地表多样性指数、地表优势度指数、地表均匀度指数和地表破碎度指数、类型斑块密度、类型边缘密度、类型多样性指数、类型优势度指数、类型均匀度指数、类型破碎度指数、斑块平均面积、最大和最小斑块面积、斑块面积标准差、斑块面积变动系数、斑块圆度、斑块矩度、斑块形状复杂度，共19个定量指标进行表征。基本计算方法是，从地理国情普查组中获取统计单元内的地表覆盖图层，并遍历图层中的各个要素矢量图斑，分别读取各矢量图斑的几何形状和要素类型信息，按照以下公式计算上述指标值。

(1)地表斑块密度等于统计单元内地表覆盖图层中所有的斑块总数除以统计单元的总面积，反映统计单元内地理国情要素斑块在空间上的分布密度。计算公式为：

PD = \frac{1}{A} Σ_{i = 1}^{M} N_{i};

式中：PD为统计单元内的地表斑块密度，M为统计单元内的地表覆盖类型的个数，N_i为统计单元内第i种地表覆盖类型的斑块个数(在景观格局分析和空间分析中，斑块可以理解为图斑或地块，即用一个封闭多边形表示空间上连片的、具有相同地表覆盖类型的区域，对应地理国情普查数据集中地表覆盖要素图层中的一个图斑)，A为统计单元的地表总面积。

(2)地表边缘密度反映统计单元内不同类型地理国情要素斑块间的边界长度。计算公式为：

ED = \frac{1}{A} Σ_{i = 1}^{M} Σ_{j = 1}^{M} P_{ij};

(3)地表多样性指数是基于信息论基础上，用来度量一定范围内地理国情要素系统的空间结构组成复杂程度的指数。计算公式为：式中，H为多样性指数，m为地表覆盖类型总数。p_k为地表覆盖类型面积比例，即某类地表覆盖类型斑块的总面积与统计单元内所有地表覆盖类型斑块总面积的比值。当m＝1，说明一个地区只有一种地表覆盖类型，此时H＝0，多样性指数最小。如果某一个地区地表覆盖类型越多，则m越大。当m→∞时，多样性指数H达到最大值。

(4)地表优势度指数用于测定地表空间格局构成中一种或几种类型支配整个地表构成的程度，计算公式为：

D = H_{\max} + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k}) = {Iog}_{2} (m) + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

(5)地表均匀度指数用以反映地表覆盖空间格局构成中各类型分布的均匀程度，均匀度可表达为：

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)} .

(6)地表破碎度指数反映地表覆盖景观的破碎化程度，即地表覆盖被分割的破碎程度。公式为：C＝N/A；N为区域地表覆盖图斑总数，A为整个区域土地总面积。

(7)类型斑块密度等于类型斑块总数除以类型总面积，用于反映单元内某种地表覆盖类型的斑块分布状态。计算公式为：

式中：PD_i为类型斑块密度，N_i为研究范围内类型斑块总数，A_i为研究范围某类型总面积。

(8)类型边缘密度为统计单元内单位面积上不同类型景观要素斑块间的边缘长度。计算公式为：

式中：P_ij为景观中第i类景观要素斑块与相邻第j类景观要素斑块间的边界长度。

(9)类型多样性指数是基于信息论基础上，用来度量类型结构组成复杂程度的指数。计算公式为：

式中，H为多样性指数，m为某类型图斑总数。p_k为斑块面积比例，即某类地表覆盖类型斑块面积与该类型斑块总面积的比值。

(10)类型优势度指数用于测定某类型图斑构成中一个或几个图斑支配整个类型构成的程度。计算公式为：

D_{t} = H_{t \max} + Σ_{i = 1}^{n} P_{i} {Iog}_{2} (P_{i}) = {Iog}_{2} (n) + Σ_{i = 1}^{n} P_{i} {Iog}_{2} (P_{i});

式中，D_t为类型t的优势度指数；n为类型t的图斑个数；H_max是Shannon多样性指数最大值，P_i为某类型图斑i占该地类总面积的比例。

(11)类型均匀度指数用以反映某类型构成中各图斑大小分布的均匀程度，均匀度可表达为：

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)};

(12)类型破碎度指数反映地表覆盖景观的破碎化程度，即某类地表覆盖被分割的破碎程度。公式为：C＝N_i/A_i；

式中，N_i为区域某类图斑总数，A_i为某类地表总面积。

(13)斑块平均面积表示斑块平均大小。公式为：

式中，N_i为第i类斑块总数；A_ij为第i类第i个斑块的面积。

(14)最大和最小斑块面积表示统计单元内斑块面积的最大最小值。

(15)斑块面积标准差表示统计单元内某类斑块面积的统计标准差。

(16)斑块面积变动系数表示统计单元内斑块规模的变异程度。公式为：

C_{i} = \frac{S_{i}}{A_{i}} \times 100 %;

(17)斑块圆度表示斑块与圆的偏离程度。公式为：

(18)斑块矩度表示斑块与正方形的偏离程度。公式为：

(19)斑块形状复杂度采用斑块分维数表示斑块的形状复杂程度。公式为：

P = {kA}^{D / 2}, D = 2 \ln (\frac{P}{k}) / \ln (A);

式中，P是斑块的周长，A是斑块的面积，D是分维数，k是常数。一般地说，欧几里德几何形状的分维为1，具有复杂边界斑块的分维则大于1，但小于2。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims

1.一种自然地表特征表征方法，其特征在于，包括以下步骤：

2.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤1中所述的空间范围是指县及县级以上的行政区或特定的流域、城市圈、城市群、经济带类在地域空间上连片的、封闭的一个区域；所述的统计单元包括标准网格单元和行政区划单元，是地理国情评价与表征的基本单位。

3.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤2中所述的地形指数(TI)的计算方法包括以下子步骤：

f_{u} = (1 - \frac{t_{u}}{T}) \times 100;

TI＝W1×f_ele+W2×f_TA+W3×f_TR+W4×f_TP；

4.根据权利要求3所述的自然地表特征表征方法，其特征在于：步骤2.7中所述的权重采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。

5.根据权利要求3所述的自然地表特征表征方法，其特征在于：步骤2.7中所述的权重：W1取10，W2取30，W3取40，W4取20。

6.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤3中所述的地貌指数计算公式如下：

FI＝(1×f_平原+0.8×f_盆地+0.6×f_丘陵-f_高原-f_山地+1)×50；

7.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤4中所述的地表覆盖指数计算公式如下：

CI＝(f_耕+f_园+0.5×(f_林+f_草+f_水)-f_荒+1)×50；

8.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤5中所述的自然地表指数计算公式如下：

LI＝W_t×TI+W_f×FI+W_c×CI；

其中W_t、W_f和W_e分别对应地形指数、地貌指数和地表覆盖的权重，取值范围为(0-100)。

9.根据权利要求8所述的自然地表特征表征方法，其特征在于：步骤5中所述的权重采用特尔斐法(专家打分法)或层次分析法(AHP法)进行确定。

10.根据权利要求8所述的自然地表特征表征方法，其特征在于：步骤5中所述的权重：W_t取30，W_f取30，W_c取40。

11.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤6中所述的区域总体格局表征包括地表斑块密度、地表边缘密度、地表多样性指数、地表优势度指数、地表均匀度指数和地表破碎度指数共6个定量指标；

PD = \frac{1}{A} Σ_{i = 1}^{M} N_{i};

ED = \frac{1}{A} Σ_{i = 1}^{M} Σ_{j = 1}^{M} P_{ij};

H = - Σ_{k = 1}^{m} P_{k} I {og}_{2} (P_{k});

D = H_{\max} + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k}) = {Iog}_{2} (m) + Σ_{k = 1}^{m} P_{k} {Iog}_{2} (P_{k});

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}}{\lg (m)};

C＝N/A；

12.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤6中所述的类型格局表征包括类型斑块密度、类型边缘密度、类型多样性指数、类型优势度指数、类型均匀度指数、类型破碎度指数共6个定量指标；

{PD}_{i} = \frac{N_{i}}{A_{i}};

E D_{i} = \frac{1}{A_{i}} Σ_{j = 1}^{M} P_{ij};

H = - Σ_{k = 1}^{m} P_{k} I {og}_{2} (P_{k});

D_{t} = H_{t \max} + Σ_{i = 1}^{m} P_{i} {Iog}_{2} (P_{i}) = {Iog}_{2} (n) + Σ_{i = 1}^{n} P_{i} {Iog}_{2} (P_{i});

E = \frac{H}{H_{\max}} = \frac{- \lg [Σ_{k = 1}^{m} {(P_{k})}^{2}]}{\lg (m)};

C＝N_i/A_i；

式中，N_i为区域某类图斑总数，A_i为某类地表总面积。

13.根据权利要求1所述的自然地表特征表征方法，其特征在于：步骤6中所述的斑块格局表征包括斑块平均面积、最大和最小斑块面积、斑块面积标准差、斑块面积变动系数、斑块圆度、斑块矩度、斑块形状复杂度共7个定量指标；

所述的斑块平均面积表示斑块平均大小，计算公式为：

{\overset{&OverBar;}{A}}_{i} = \frac{1}{N_{i}} Σ_{j = 1}^{N_{i}} A_{ij};

C_{i} = \frac{S_{i}}{A_{i}} \times 100 %;

所述的斑块圆度表示斑块与圆的偏离程度，计算公式为：

D = P / 2 \sqrt{πA};

D = 0.25 P / \sqrt{A};

P = {kA}^{D / 2}, D = 2 \ln (\frac{P}{k}) / \ln (A);