地理格网下不易获取地理要素空间渐变信息的数字制图方法

摘要

地理格网下不易获取地理要素空间渐变信息的数字制图方法，利用不易获取地理要素与易获取地理要素间协同变化关系，在地理格网下由易获取地理要素数据集生成反映不易获取地理要素空间渐变信息的数字地图。技术方案为：准备与不易获取地理要素存在协同变化的易获取地理要素数据集；确定不易获取地理要素类型划分；针对各类型构建反映全局协同变化和反映局域协同变化的关系曲线；利用关系曲线和已有数据集生成不易获取地理要素各类型模糊隶属度图；综合各模糊隶属度图生成不易获取地理要素的数字属性图和类型图。本发明克服了目前制图方法所需样本量大、应用条件苛刻等问题，可用于自然资源制图、自然灾害危险性制图等领域。

说明书

技术领域

本发明涉及一种地理格网下对不易获取的地理要素的空间渐变信息进行数字制图的方法，技术应用领域为自然资源制图、自然灾害危险性制图等。

背景技术

反映地理要素空间渐变信息的数字地图是自然资源管理、土地利用规划、灾害应急管理所需要的重要基础数据。从数字制图难易程度上来说，一些易获取地理要素可以通过遥感等现代观测手段得到(如地形，许多植被要素)，或是由于其空间连续性强而可以通过少量观测点的数据进行内插而得到(如温度场)，因而对其空间渐变信息进行制图较容易；与此相对的，一些不易获取地理要素不容易被直接观测到(如滑坡危险度)，或是由于其空间变化复杂、野外采样的成本过高(如土壤)，在实际应用中无法利用遥感技术和空间内插方法有效精细地定量其空间变化，针对这些地理要素的数字制图也较为困难，目前对这些地理要素所产生的数字地图无论在空间分辨率上还是在属性的详细程度上均难以达到易获取地理要素数字地图的精细程度，难以支持许多实际应用需求。

目前，已有的对不易获取地理要素空间渐变信息进行制图的技术方法主要有两类：一是基于大量野外样本的地统计及空间插值方法；一是基于自然过程机理进行模拟的方法。前者在原理上要求所应用的对象需满足一定的前提假设(例如平稳假定)，而这种前提假设对于许多不易获取地理要素常常无法满足，此外，这一类方法在实际应用中需要大量样本，而通常对不易获取地理要素的采样代价很大，当对空间及属性精度有很高的要求时，相应的高昂采样代价显著降低了这类方法对不易获取地理要素进行制图的实用性。利用自然过程机理进行模拟的方法需要对于待刻画的不易获取地理要素的形成机理有足够的理解，通过建立机理过程方程模拟不易获取地理要素的空间变化从而成图。而实际情况多是人们对机理过程难以达到这类方法所要求的理解程度，因而这种方法的在实际制图中的实用性受到很大限制。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种地理格网下不易获取地理要素空间渐变信息的数字制图方法，该方法利用易获取地理要素(或称环境条件)对应的格网数据，通过建立反映易获取地理要素和不易获取地理要素间协同变化的关系曲线，可以充分利用现有数据资源，快速、准确的产生数字形式的、反映不易获取地理要素的空间渐变信息的属性图和类型图。

本发明的技术解决方案为：地理格网下不易获取地理要素空间渐变信息的数字制图方法，步骤如下：

(1)首先选择与不易获取地理要素存在协同变化的易获取地理要素，并根据最后制图的需求确定易获取地理要素的格网数据的分辨率，即像素大小，将已有易获取地理要素数据转换为该分辨率下的格网格式，形成易获取地理要素数据集；所述的易获取地理要素指可以通过遥感、地理信息分析技术等直接得到的地理要素，如地形、植被类型，所述的不易获取地理要素指难以通过遥感、地理信息技术等直接获取的地理要素，如土壤、滑坡危险度等；

(2)对不易获取地理要素进行类型划分，即确定产生的数字地图要反映哪些类型的空间渐变；

(3)针对步骤(2)划分的每一类型，逐一建立反映易获取地理要素与该类型之间全局协同变化和局域协同变化的关系曲线，从而构建关系曲线库；所述的关系曲线有三种类型：

基于连续函数，其关系曲线公式如下：

$f\left(x\right)=\frac{1}{1+{\left(\frac{x-b}{d}\right)}^2}$

其中，f是模糊隶属度值；x是待制图位置处的易获取地理要素值；b是模糊隶属度为1，即出现典型不易获取地理要素类型时的易获取地理要素值；d是当模糊隶属度为0.5，即出现不易获取地理要素类型过渡情况时的易获取地理要素值与b之差，d≠0；

基于三次样条函数，将模糊隶属度已经确定的关键点连接起来形成关系曲线；

基于二值阶跃函数的关系曲线，仅取0和1两个值；

所述的关系曲线库的组织方式是：一个不易获取地理要素类型对应于一个或多个环境条件组合，每一个环境条件组合相应包括了一条或多条关系曲线，环境条件组合分为全局环境条件组合和局域环境条件组合两类，全局环境条件组合普遍作用于整个制图区域，局域环境条件组合作用于以某点为中心的一定地理区域范围，该范围用中心点空间坐标、作用半径和距离衰减率确定；

(4)将步骤(3)确定的关系曲线库，结合步骤(1)确定的易获取地理要素数据集进行计算，产生不易获取地理要素每一类型的模糊隶属度图；所述的计算过程是：对于一个类型下各个全局环境条件组合，首先计算制图区域某一像素位置上相对于该全局环境条件组合下每条关系曲线的模糊隶属度，然后取模糊隶属度最小值作为该像素位置上相对于该全局环境条件组合的模糊隶属度，最后取各全局环境条件组合下得出的模糊隶属度的最大值作为该像素位置相对于所有全局环境条件组合所形成集合的模糊隶属度，如果该类型下没有局域环境条件组合，则该模糊隶属度即为该像素位置相对于该类型的最终模糊隶属度；对于一个类型下的各个局域环境条件组合，如果该像素落入其作用半径内，首先按照与全局环境条件组合相同的办法求得该像素位置上相对于该环境条件组合的模糊隶属度，然后根据该像素与局域环境条件组合中心点的距离，利用距离衰减率对所得模糊隶属度进行纠正作为该像素位置上相对于该局部环境条件组合的模糊隶属度，若像素落入多个局域环境条件组合作用半径内，则取该像素位置相对于各局域环境条件组合的模糊隶属度的最大值作为该像素位置相对于该类型的最终模糊隶属度；循环对制图区域各像素位置进行计算，最后每个像素位置都会获得对于各个类型的模糊隶属度，从而形成每一类型的模糊隶属度图；

(5)将步骤(4)中各个类型的模糊隶属度图中的模糊隶属度作为权值，对每个类型对应的典型属性进行加权平均，以此得到每个像素所在位置的属性值，从而得到地理格网数据结构下不易获取地理要素的数字属性图；或对步骤(4)产生的各个类型的模糊隶属度图进行硬化处理，即将某位置的最终类型赋为该位置上拥有最大模糊隶属度的类型，从而得到地理格网数据结构下不易获取地理要素的数字类型图。

本发明的原理是：(1)虽然不易获取地理要素的空间渐变信息难以直接定量得到，但是它们与许多易获取地理要素存在着协同变化，可以利用这种协同变化关系进行制图；(2)易获取地理要素的数据源方面，迅速发展的遥感技术和地理信息技术可以提供具有较高精细程度、地域覆盖面广的数据，为本发明步骤1中的易获取地理要素数据集收集与准备提供了保证；(3)反映不易获取地理要素与易获取地理要素之间协同变化的关系已经以多种形式存在，这些关系可以通过转化为本发明步骤3所述的关系曲线的形式后用于制图生产，且利用关系曲线来制图考虑了地理现象的过渡性，这与大多数地理要素在空间呈现连续渐变特点是相吻合的；(4)步骤3中局域环境条件组合的地学含义是当易获取地理要素的某一种类型除了在整个区域上的分布规律外，还在某一具体空间位置上出现分布特例时，该特例出现位置周围的局部空间范围上的环境组合与该特例所对应的环境组合越相似，并且空间距离越近，则模糊隶属度越高；(5)地理要素受多个因子影响时，限制因子往往起决定作用，因此，本发明步骤4中求得某一像素位置上相对于某环境条件组合下各条关系曲线的模糊隶属度后，取模糊隶属度最小值作为该像素位置相对于该环境条件组合的模糊隶属度。

本发明与现有技术相比的优点在于：

(1)现有大多数不易获取地理要素空间渐变信息的制图方法对已有的地理要素数据没有充分利用。随着遥感技术和地理信息技术的迅速发展，易获取地理要素的数据精细程度、丰富程度上都有了很大提高，本发明充分利用了现有数据资源，可以在地理格网下准确、快速地对不易获取地理要素空间渐变信息进行数字制图。

(2)现有技术中，基于地统计与空间插值的方法对样点数目、应用前提等都有较高要求，而基于自然过程机理模拟的方式又需要对机理过程有深刻把握，因此这些技术的可应用性受到很大限制。本发明充分利用了不同地理要素间的协同变化关系，并以关系曲线的形式将其表达，能够和已有数据结合完成制图，应用的要求较低，实用性强。

(3)本发明在实施过程中所积累的关系曲线具有较好的适用性，可以在不同地区进行重复利用或经过调整后再次利用，从而避免了基于统计的方法中每换一个区域都要进行大量采样来构建统计关系的问题，大大加快了制图流程。

(4)本发明由一套数据产出多种数字地图产品，发明实施的最终结果是地理格网数据结构下的不易获取地理要素的数字属性图和数字类型图，实施过程中还能附带产生不易获取地理要素各类型模糊隶属度图，多种形式的产品使得结果的应用范围更加广泛。

附图说明

图1为本发明方法的流程图；

图2为本发明中的钟形、S形、Z形关系曲线示意图；

图3为本发明的关系曲线库的组织结构图。

具体实施方式

如图1所示的流程图，以数字土壤制图为例，本发明的具体实施方法如下：

1.土壤为不易获取地理要素，选择与制图区域土壤存在协同变化的易获取地理要素(如高程、坡度、坡向等)，根据最后制图的需求确定格网数据分辨率(如10m)，将已有易获取地理要素数据转换成该分辨率下的格网格式，形成易获取地理要素数据集；

2.对制图区域的土壤进行类型划分，即确定产生的数字地图要反映哪些土壤类型的空间渐变；

3.针对步骤2划分的每一种土壤类型，逐一建立反映易获取地理要素与该类型之间全局协同变化和局域协同变化的关系曲线，从而构建关系曲线库；具体方法为：对于每一种土壤类型，指定该种类型土壤分布于哪些环境条件组合，并确定环境条件组合是全局的还是局域的，对于局域作用的环境条件组合，指定其中心点空间坐标、作用半径以及距离衰减率；对于每个环境条件组合，建立一条或多条关系曲线，每条关系曲线反映的是该土壤类型与一个易获取地理要素(如高程)之间的协同变化关系，多条关系曲线反映一种土壤类型与一种环境条件组合情况(如高程、坡度、坡向)之间的协同变化关系。

关系曲线主要有三种类型，即基于连续函数的、基于三次样条函数的、和基于二值阶跃函数的：

(1)基于连续函数的关系曲线公式如下：

$f\left(x\right)=\frac{1}{1+{\left(\frac{x-b}{d}\right)}^2}$

其中，f是模糊隶属度值；x是待制图位置处的易获取地理要素(如坡度)的值；b是模糊隶属度为1，即出现典型的某一土壤类型时的易获取地理要素值(如坡度为6°)；d是当模糊隶属度为0.5，即出现土壤类型过渡时的易获取地理要素值与b之差，d≠0。当d值不同时，连续函数表现出的曲线形式不同，如图2所示，当x>b和x<b时，d取相同有限值的情况下，连续函数表现为对称的钟形曲线，而当x>b和x<b时，d取不同有限值的情况下，连续函数表现为不对称的钟形曲线，钟形曲线表示当易获取地理要素的数值偏离隶属度为1的值时，土壤相对于典型土壤的模糊隶属度随之降低；当x>b时，d取为无穷大，并且x<b时，d取值为有限值的情况下，连续函数表现为S形曲线，S形曲线表示仅当易获取地理要素的数值小于对应隶属度为1的最小值时，土壤相对于典型土壤的模糊隶属度小于1，此时偏离距离越远，模糊隶属度越低；当x<b时，d取值为无穷大，并且x>b时，d取值为有限值的情况下，连续函数表现为Z形曲线，Z形曲线表示仅当易获取地理要素的数值大于对应隶属度为1的最大值时，土壤相对于典型土壤的模糊隶属度小于1，此时偏离距离越远，模糊隶属度越低。

如果可以得到上述公式中的b，d两个参数，则可以直接构建基于连续函数的关系曲线。例如已知某种土壤类型典型出现于坡度为6°处，在坡度为9°和2°处呈现明显过渡特性，则公式中b取为6，d在x>b和x<b时分别取3和4(9-6＝3，6-2＝4)，从而形成不对称钟形关系曲线。

(2)如果无法得到出现典型的某一土壤类型时易获取地理要素的值，或者该土壤类型与易获取地理要素协同变化关系较为复杂时，用基于连续函数的关系曲线难以表达，此时利用反映该土壤类型与易获取地理要素关系的一些模糊隶属度已经确定的关键点，采用三次样条函数将各关键点连接起来形成关系曲线。

(3)如果选择的易获取地理要素是岩石类型、植被类型等一些类型要素，则建立二值阶跃函数类型的关系曲线，函数仅取0和1两个值。如在岩石类型为1，3，4(1，3，4为岩石类型数字编号)的条件下产生某一土壤类型，则该土壤类型对应岩石类型这一易获取地理要素的关系曲线仅在1，3，4下取值为1，其余情况下取值为0。

各条关系曲线按图3所示的层次结构进行组织。一个土壤类型下有一个或多个环境条件组合，每一个环境条件组合相应包括了一条或多条关系曲线，环境条件组合分为全局环境条件组合和局域环境条件组合两类，全局环境条件组合普遍作用于整个制图区域，局域环境条件组合应用于以某点为中心的一定地理区域范围，该范围用中心点空间坐标、作用半径和距离衰减率确定

4.将步骤3建立的关系曲线库，结合步骤1中易获取地理要素数据集进行计算，产生每一土壤类型的模糊隶属度图，计算过程如下：

对于某一土壤类型下各个全局环境条件组合，首先计算制图区域某一像素位置上相对于该全局环境条件组合下每条关系曲线的模糊隶属度，然后取模糊隶属度最小值作为该像素位置相对于该全局环境条件组合的模糊隶属度(如某一全局环境条件组合有高程、坡度、坡向三条关系曲线，将某像素位置的高程值、坡度值、坡向值分别代入这三条关系曲线中计算得到三个模糊隶属度值，取其中最小值作为该像素位置相对与该环境条件组合的模糊隶属度)，最后取各全局环境条件组合下得出的模糊隶属度的最大值作为该像素位置相对于所有全局环境条件组合所形成集合的模糊隶属度(如某一土壤类型下有两个全局环境条件组合，每一个像素将会获得相对与这两个全局条件组合的两个模糊隶属度值，假设两个值分别是0.8，0.3，则取较大的0.8作为像素相对与该土壤类型下所有全局环境条件所形成集合的模糊隶属度)，如果该土壤类型下没有局域环境条件组合，则该模糊隶属度(如上述例子中的0.8)即为该像素位置相对于该土壤类型的最终模糊隶属度；

对于一个类型下的各个局域环境条件组合，如果该像素落入其作用半径内，首先按照与全局环境条件组合相同的办法求得该像素位置上相对于该环境条件组合的模糊隶属度，然后根据该像素与局域环境条件组合中心点的距离，利用距离衰减率对所得模糊隶属度进行纠正作为该像素位置上相对于该局部环境条件组合的模糊隶属度，计算公式如下：

$S_{\mathrm{ij},t}\prime =\frac{S_{\mathrm{ij},t}\left(\frac{1}{{(d_{\mathrm{ij},t}\prime )}^r}\right)+S_k\left(\frac{1}{{\left(d_{\mathrm{ij},k}^{\prime }\right)}^r}\right)}{\left(\frac{1}{{(d_{\mathrm{ij},t}\prime )}^r}\right)+\left(\frac{1}{{\left(d_{\mathrm{ij},k}^{\prime }\right)}^r}\right)}$

其中，S_ij，t’表示该土壤类型在像素(i，j)位置上相对于局域环境条件组合t的模糊隶属度值；S_ij，t表示按照与全局环境条件组合相同方法计算得到的像素(i，j)位置上相对于局域环境条件组合t纠正前的模糊隶属度；d_ij，t’表示像素(i，j)与局域环境条件组合t中心点的地表距离；根据所设置的局域环境条件组合作用半径w_t，可以计算出其的作用边界，局域环境条件组合t中心点和像素(i，j)可以形成一条直线，该直线与作用边界有两个交点，k是与像素(i，j)在局域环境条件组合t中心点同一侧的交点，d_ij，k’表示像素(i，j)位置与k的地表距离；S_k是k所在位置相对于所有全局环境条件所形成集合的模糊隶属度值，该值在针对全局环境条件进行计算时已经得到；r表示距离衰减率；若像素落入多个局域环境条件组合作用半径内，则取该像素位置相对于各局域环境条件组合的模糊隶属度的最大值作为该像素位置相对于该类型的最终模糊隶属度；

循环对制图区域各像素位置进行计算，最后每个像素位置都会获得对于各个类型的模糊隶属度，从而形成每一类型的模糊隶属度图；

5.对步骤4得到的各个土壤类型的模糊隶属度图进行不同处理，可以分别得到土壤属性图或土壤类型图。例如要对有机碳含量这一土壤属性进行制图，通过查阅文献或检索知识表得到每个土壤类型对应的有机碳含量典型值，将步骤4中产生的土壤类型模糊隶属度图中的模糊隶属度值作为权值对各类型对应的有机碳含量典型值进行线性加权平均，这样就可以估算出每个像素所在位置上的有机碳含量，从而得到地理格网数据结构下的土壤有机碳含量属性图，线性加权平均采用如下公式：

$V_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\Sigma }}{S^k}_{\mathrm{ij}}·V^K}{\underset{k=1}{\overset{n}{\Sigma }}{S_{\mathrm{ij}}}^k}$

其中，V_ij表示像素(i，j)位置的有机碳含量，V^k代表土壤类型k有机碳含量的典型值，S_ij^k代表土壤类型k在像素(i，j)位置的模糊隶属度值，n是制图区域内土壤类型数目；如果要产生土壤类型图，则对步骤4产生的各个土壤类型的模糊隶属度图进行硬化处理，即将某像素位置的最终类型赋为该像素位置上拥有最大模糊隶属度的土壤类型，从而得到地理格网数据结构下的土壤类型图。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。