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

Abstract

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

Description

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

技术领域

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

背景技术

如何准确、精细地分析和表达预测变量和响应变量的空间互作关系，一直是人们感兴趣的问题。目前，描述预测变量与响应变量的空间关系的方法有两种：全局线性模型和局部线性模型。全局线性模型假定预测变量和响应变量的空间关系是平稳不变，这实际上违背了预测变量和响应变量的空间互作的局部变异性。局部线性模型考虑空间变量之间关系的局部变异，但是该模型中预测变量对响应变量影响的强弱依赖空间对象之间的空间距离；即空间邻近则预测变量对响应变量的影响越强，反之，则越弱。这实际也有违“空间邻近，空间对象的非空间属性未必相似”的客观事实。用以下示例来说明该问题。给定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)={\mathrm{\&beta;}}_0\left(U_i\right)+\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&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(\mathrm{\&beta;}\right)=\underset{j=1}{\overset{15}{\mathrm{\&Sigma;}}}W\left(U_j^i\right){[Y\left(O_j\right)-{\mathrm{\&beta;}}_0\left(U_i\right)-\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&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)/\mathrm{\&alpha;}]}^2\}---\left(4\right)$

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

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

$\widehat{\mathrm{\&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。

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

Claims (1)

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

（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)，约定在显著水平<0.05的情况下对预测变量P_j与Y的相关性进行显著性检验，如果P_j与Y相关性显著，则保留P_j，否则剔除P_j；

（2）邻近点搜索

网格化研究区域，获取待观测样点集U以及搜索最邻近未观测样点的一定量的观测样点；待观测样点集定义为所有未取样位点的集合，网格单元大小约定为20米×20米，邻近观测点数为15；首先按20米×20米网格大小对研究区域进行网格化，每个网格单元的中心位点作为一个待观测样点，提取所有网格单元的中心位点构成未观测样点集U＝{U₁,U₂,...,U_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[(U_i\left(X_i\right)-O_j(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)={\mathrm{\&beta;}}_0\left(U_i\right)+\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_l\left(U_i\right)P_l\left(U_i\right)(1\&le;l\&le;k<p)---\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(\mathrm{\&beta;}\right)=\underset{j=1}{\overset{15}{\mathrm{\&Sigma;}}}W\left(U_j^i\right){[Y\left(O_j\right)-{\mathrm{\&beta;}}_0\left(U_i\right)-\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&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)/\mathrm{\&alpha;}]}^2\}---\left(4\right)$

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

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

$\widehat{\mathrm{\&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)，计算所有未观测点的预测变量与响应变量的回归系数，将回归系数进行可视化输出。