一种空间变量相关性的精细表达方法

摘要

本发明属于预测变量与响应变量的空间关系分析技术领域，具体公开了一种空间变量相关性的精细表达方法，适用土壤环境、农业、气象、水文、地理科学等涉及空间变量关系分析和表达的研究。本发明方法以预测变量与响应变量的空间关系为参照，用局部线性模型拟合预测变量与响应变量之间的空间关系，并预测和输出回归系数。本发明的优点在于：用基于预测变量相似度量的局部线性模型对空间关系进行建模和拟合，保证了空间变量相关性的精细表达。

说明书

技术领域

本发明涉及预测变量与响应变量的空间关系分析的技术领域，更具体涉及一 种空间变量相关性的精细表达方法，适用土壤环境、农业、气象、水文、地理科 学等涉及空间变量关系分析和表达的研究。

背景技术

如何准确、精细地分析和表达预测变量和响应变量的空间互作关系，一直是 人们感兴趣的问题。目前，描述预测变量与响应变量的空间关系的方法有两种： 全局线性模型和局部线性模型。全局线性模型假定预测变量和响应变量的空间关 系是平稳不变，这实际上违背了预测变量和响应变量的空间互作的局部变异性。 局部线性模型考虑空间变量之间关系的局部变异，但是该模型中预测变量对响应 变量影响的强弱依赖空间对象之间的空间距离；即空间邻近则预测变量对响应变 量的影响越强，反之，则越弱。这实际也有违“空间邻近，空间对象的非空间属 性未必相似”的客观事实。用以下示例来说明该问题。给定20个4维（地理坐标 X和Y及两个非空间属性）空间样点，分布在A、B、C三个空间异质的单元中，即 三个单元之间的非空间属性差异明显。分别按空间距离相似及非空间属性相似模 式搜索a、b、c三个样点的5个邻近样点，相应结果见图1-a和图1-b。由图1-a可知， 基于空间距离邻近不能保障邻近b或c的5个样点的非空间属性相似；而基于非空 间属性相似则保障了邻近b或c的5个样点的非空间属性相似（图1-b）。这表明基 于空间邻近的局部线性模型不能可靠地描述变量之间的空间关系，而用非空间属 性相似有利于准确描述变量之间的空间关系。

为此，本发明拟提出一种空间变量相关性的精细表达方法。该方法以预测变 量与响应变量的空间互作关系为参照，构建基于预测变量相似的局部线性模型对 变量的空间关系进行建模，实现预测变量与响应变量的空间关系的精细表达。

发明内容

针对现有技术中存在的不足，本发明的目的在于提供了一种空间变量相关性 的精细表达方法，方法易行，操作简便。处理的空间数据对象要求包含地理坐标 和其他非空间属性。预测变量和响应变量都属于非空间属性，同时预测变量等同 于统计学上的自变量而响应变量则对应因变量。

为了实现上述的目的，本发明采用以下技术措施：

本发明的内容包括：提取预测变量、邻近点搜索、局部线性模型构建及回归 系数求解和输出四部分。

总体上，本发明方法的流程如下：

1、分析所有候选预测变量与响应变量的相关性，并进行F检验，与响应变 量相关性显著（显著水平p<0.05）的候选预测变量作为预测变量；

2、网格化研究区域，建立待观测点集，对每个待观测点，搜索它的15个最 邻近的已知观测点；

3、对每个待观测点，用本发明提出的局部线性模型对其最邻近的15个已知 观测点的预测变量与响应变量的空间关系进行拟合，建立预测模型；

4、用预测模型预测和输出所有待观测点的预测变量与响应变量的回归系数， 将输出结果转成栅格数据，以栅格图层显示。

具体的，一种空间变量相关性的精细表达方法，其步骤是：

（1）提取预测变量

在实际研究中，与响应变量相关的预测变量可能有很多。将与响应变量相关 性弱的预测变量引入模型，会降低模型的准确性。因此，需要剔除与响应变量相 关性弱的预测变量。

给定观测样点集S＝{S₁,S₂,...,S_n}，n×(p+3)矩阵M＝{G,X,Y}，任意观测样 点S_i＝{G_m,X_m,Y_m}(1≤m≤n)，其中，n是观测样点总数，p是候选预测变量数目； G是所有观测样点的地理坐标构成的n×2矩阵，G_m是观测样点m的地理坐标组成 的2维行向量；X是所有观测样点的p个候选预测变量值构成的n×p矩阵，X_m是观 测样点m的p个候选预测变量值构成的p维行向量；Y是所有观测样点的响应变量 值构成的n维列向量；Y_m是观测样点m的响应变量值。

设置候选预测变量矩阵X＝{P₁,P₂,...,P_p}，P_j是是第j个候选预测变量值构成的 n维列向量(1≤j≤p)。在显著水平P<0.05对预测变量P_j与Y的相关性进行显著性 检验（F检验），如果P_j与Y相关性显著，则保留P_j，否则剔除P_j。

（2）邻近点搜索

邻近点搜索涉及网格化研究区域，获取待观测样点集U以及搜索最邻近未观 测样点的一定量的观测样点。待观测样点集定义为所有未取样位点的集合，网格 单元大小约定为20米×20米，邻近观测点数为15。首先按20米×20米网格大小对 研究区域进行网格化，每个网格单元的中心位点作为一个待观测样点，提取所有 网格单元的中心位点构成未观测样点集U＝{U₁,U₂,...,U_N}（N>n,N是待观测样点 总数，n是观测样点总数，N和n的大小取决于应用实例；进一步将待观测点集U 与预测变量图层叠加以获取相应的预测变量；最后，对任意待观测点U_i(1≤i≤N), 基于预测变量相似搜索最邻近U_i的15个观测样点集O。定义预测变量相似为点对 在预测变量空间上的欧几里德距离例如，U_i到它一个邻近点O_j(1≤j≤15)在预测变 量空间上距离dist(U_i,O_j)，按公式（1）计算：

$\mathrm{dist}(U_i,O_j)=\sqrt{[U_i\left(X_i\right)-O_j\left(X_j\right)]{\left[\right(U_i\left(X_i\right)-O_j\left(X_j\right)]}^T}---\left(1\right)$

其中，U_i(X_i)表示待观测样点U_i的预测变量值构成的行向量，O_j(X_j)表示观测 样点O_i的预测变量值构成的行向量，T表示矩阵转置操作。

（3）局部线性模型构建

对待观测点U_i，用公式(2)拟合它的预测变量与响应变量的空间关系

$Y\left(U_i\right)={\beta }_0\left(U_i\right)+\underset{l=1}{\overset{k}{\Sigma }}{\beta }_l\left(U_i\right)P_l\left(U_i\right)(1\le l\le k<p,1\le j\le 15)---\left(2\right)$

公式(2)中，k是预测变量数目，p表示候选预测变量数目；Y(U_i)表示U_j的响 应变量，P_l(U_i)是U_i的第l个预测变量，β₀(U_i)是截距，β_l(U_i)是U_i的第l个预测变 量与响应变量的回归系数；

β₀(U_i)和β_l(U_i)是未知的，需要用最邻近U_i的15个观测样点，并通过用公式 (3)进行求解：

$\mathrm{WSRS}\left(\beta \right)=\underset{j=1}{\overset{15}{\Sigma }}W\left(U_j^i\right){[Y\left(O_j\right)-{\beta }_0\left(U_i\right)-\underset{l=1}{\overset{k}{\Sigma }}{\beta }_l\left(U_i\right)P_l\left(O_j\right)]}^2---\left(3\right)$

公式(3)中，Y（O_j）是邻近点O_j的响应变量值，P_l(O_j)是邻近点O_j的第l个预 测变量，是邻近点O_j对U_i的权重系数，按公式（4）计算：

$W\left(U_j^i\right)=\exp \{-0.5*{[\mathrm{dist}(U_i,O_j)/\alpha ]}^2\}---\left(4\right)$

公式(4)中，dist(U_i,O_j)按公式(1)计算，α参数取U_i与最邻近的15个观测点 的距离的平均值；

对公式（3）两端求导，可导出回归系数的解的矩阵表达式如下：

$\hat{\beta }={\left[P^{T}W\left(U_i\right)P\right]}^{-1}{P}^{T}W\left(U_i\right)Y---\left(5\right)$

公式（5）中，W(U_i)是一个15×15矩阵，它的最后一行元素值为其 它行元素值皆为0；Y＝[Y(O₁),Y(O₂),...,Y(O_j)]^T是最邻近U_i的15个观测点的响应变 量构成的列向量；P是15×(k+1)预测变量矩阵，P＝[X(O₁),X(O₂),...,X(O_j)]^T，， 行向量X(O_j)由邻近点O_j的预测变量值和整数1构成，且 X(O_j)＝[1,P₁(O_j),P₂(O_j),...,P_l(O_j)]^T，P_l(O_j)是邻近点O_j的第l个预测变量值；i、 j及l的取值范围分别为1≤i≤N，1≤j≤15，1≤l≤k，k是预测变量数目，N是待 观测样点总数。

（4）回归系数求解及输出

利用步骤（3）建立的回归系数求解公式（5），计算所有未观测点的预测变 量与响应变量的回归系数。将回归系数进行可视化输出是反映预测变量与响应变 量的空间互作关系的最佳模式。回归系数可视化过程不在本发明的保护范围内， 可借助相关软件系统实现，如：ArcGIS、MapGIS等商业软件都能有效地实现回归 系数的可视化表达。

本发明的优点和有益效果在于：以预测变量与响应变量的空间关系为参照， 用基于预测变量相似度量的局部线性模型对空间关系进行建模和拟合，保证了空 间变量关系的精细表达。

附图说明

图1为空间距离邻近和非空间属性邻近分析（a表示空间距离邻近，b表示非 空间属性邻近）；

图2为预测变量（高程、土地利用方式）和响应变量（有机碳含量）分布图；

图3为预测变量的回归系数分布图。

具体实施方式

以下具体的实施例为将上述发明内容中记载的一种空间变量相关性的精细 表达方法的具体方法中的步骤（1）-（4）应用于具体对象，用于详细说明本发 明的技术方案，但不用来限制本发明的应用范围。

实施例1

本实施例来源自湖南省长沙市金井流域，总面积为135km²。该流域属于典 型的亚热带红壤丘陵区，主要土地利用方式为茶园、林地和稻田。2010年在该流 域取土壤样1033份，监测土壤有机碳含量。调查表明该区域内土壤有机碳含量与 土地利用方式有关，三种土地利用方式下有机碳含量的大小依次为：稻田>林地> 茶园，此外有机碳含量也与地形因子有关。为了精细表征该流域土地利用方式、 地形因子对土壤有机碳含量的影响，采用本发明方法分析有机碳含量与土地利用 方式、地形因子的空间关系。

（1）提取预测变量

在金井流域土壤有机碳含量与土地利用方式、高程、坡度及地形湿度指数有 关。分析1033个土壤监测位点的有机碳含量与预测变量（土地利用方式、高程、 坡度及地形湿度指数）的相关性，在小于0.05的显著水平对结果进行显著性方差 检验（F检验），结果显示有机碳含量与土地利用方式和高程显著相关，而与坡 度和地形湿度指数的相关性较弱。故，保留土地利用方式和高程作为预测变量。 预测变量（土地利用方式、高程）及响应变量（土壤有机碳含量）的分布图见图 2所示。

（2）邻近点搜索

邻近点搜索涉及网格化金井流域、提取待观测点集及搜索邻近待观测点集的 观测点。本实施例中，网格单元大小限定为20米×20米，邻近观测点数为15。金 井流域网格化后，生成需要预测的待观测点集U＝{U₁,U₂,...,U_n}(N＝480000)。给 定需要预测的待观测点U_i，基于预测变量相似原则搜索出最邻近U_i的15个观测 点集（O）。

（3）模型构建及回归系数求解和输出

对需要预测的待观测点U_i，将其邻近的观测点集（O）数据代入发明内容中 的公式（5），求解出U_i的预测变量与响应变量的回归系数。要求将此计算过程对 金井流域的所有待观测点执行一遍，以输出所有待观测点的回归系数。为了可视 化回归系数的空间分布，本实施例用ArcGIS软件将回归系数转化成栅格数据后， 以栅格图层进行可视化。

（4）结果比较

为了体现本发明的技术优势，将本实施例的方法与当前应用广泛的地理加 权回归方法（GWR）进行比较。GWR方法对金井流域土壤有机碳与土地利用方 式和高程的空间关系的预测结果见图3-a,b,c，本实施例的预测结果见图3-d,e,f。 由图3可知，GWR不能精细地反映金井流域土地利用方式和高程影响土壤有机碳 含量变异的空间格局，其预测结果总体上趋于平滑、辨识度低。而本实施例则更 精细地反映出了金井流域土壤有机碳含量变异与土地利用方式和高程的空间关 系。例如，本实施例预测的土地利用方式和有机碳含量变异的空间关系分布图（图 3-e）清晰地表征了金井流域土地利用方式影响有机碳含量变异的空间格局，即 回归系数的高值主要分布在稻田区，而其低值则相对分布在林地和茶园地区。同 时，本实施例预测的高程和有机碳含量变异的空间关系分布图（图3-f）也清晰 地反映了金井流域高程和有机碳含量呈负相关的客观事实。比较而言，GWR的 预测结果图则明显不清晰和精细。这表明本发明方法可精细表征预测变量与响应 变量的空间关系，效果优于GWR。

以上为本发明的最佳实施方式，凡依据本发明公开的内容，本领域的普通技 术人员能够显而易见地想到一些雷同、替代方案，均应落入本发明权利要求书的 保护范围。