基于雪水当量栅格数据的空间回归建模方法

摘要

本发明提供一种基于雪水当量栅格数据的空间回归建模方法，包括：获取雪水当量栅格数据产品并对所有数据进行预处理；对所有预处理后的数据进行划分；对分块区域构建普通最小二乘线性回归模型：筛选分块残差；步骤5：构建分块残差模型；对分块的普通最小二乘线性回归模型残差构建空间误差模型或空间滞后模型；将分块残差模型结果并入到步骤3构建的普通最小二乘线性回归模型中，得到最终的空间回归模型；若该分块未执行步骤5，则将步骤3构建的普通最小二乘线性回归模型作为该分块的最终模型；对模型进行评价和分析。本发明消除了空间自相关的影响，提高模型的精度，探究不同气象因素和地形因素对雪水当量的影响，以供后续研究和分析。

说明书

技术领域

本发明属于地学统计与空间分析的技术领域，具体涉及一种基于雪水当量栅格数据的空间回归建模方法。

背景技术

作为气候系统中的重要成员之一，冰冻圈不仅会受气候的影响，也能够影响局部地区乃至全球的气候。而其中，冰冻圈的一个重要元素——积雪，也是一种特殊的下垫面，对地表辐射平衡、能量循环和水循环也有着非常重要且显著的影响。除此之外，积雪反复地覆盖和融化不仅直接影响土壤状况，还影响了许多动物的活动节律，直接或间接地对植物群落产生了深刻的影响。然而近年来，随着全球气候不断变暖，冰雪融化速度加快，雪线持续上升，生态环境不断受到影响而恶化。且在积雪强烈消融期间，极易引起与之相关的融雪洪水、泥石流、雪崩等冰冻圈灾害。雪水当量是能够反映地表积雪量变化的重要因子，且对估算地表水资源非常重要，指当积雪完全融化后，所得到的水形成的水层的垂直深度，常用单位为毫米(mm)，是重要的积雪参数之一。与积雪覆盖面积、积雪深度等其他积雪属性不同，雪水当量能够提供雪的物质和能量信息，更能体现空气温度、降水等因素的综合影响，也是地表水文模型和气候模式中的一个重要参数。

雪水当量是积雪深度与积雪密度的乘积，故而积雪深度和积雪密度的影响因素也是影响雪水当量的因素，比如降水、空气温度、风速、植被覆盖、海拔、坡度、坡向等。对雪水当量及其影响因素进行回归建模可以探究雪水当量影响因素对于雪水当量的影响。但地理数据普遍存在空间自相关性，传统回归模型比如普通最小二乘线性回归模型往往会忽略数据的空间自相关性导致模型精度较低，而考虑空间效应的空间计量模型可以解决上述问题。空间计量模型分为空间滞后模型与空间误差模型，利用空间权重矩阵，将空间关系引入模型中，对线性回归模型进行修正。然而，空间计量模型难以直接在可进行大规模积雪监测的遥感数据集上进行应用，因为它所基于的参数的极大似然估计难以在大数据集上使用，且由于遥感数据产品的栅格尺寸一般较大，其空间权重矩阵的尺寸也会变大，计算耗时长且效率低下。因此，亟需开发一种能够消除空间自相关的影响的空间计量模型。

发明内容

本发明的目的在于针对现有技术的不足之处，提供一种基于雪水当量栅格数据的空间回归建模方法，，该方法能够消除空间自相关的影响，提高模型的精度，探究不同因素对雪水当量的影响。

为解决上述技术问题，本发明采用如下技术方案：

一种基于雪水当量栅格数据的空间回归建模方法，其特征在于，包括如下步骤：

步骤1：获取雪水当量栅格数据产品，并从各数据源获取影响雪水当量变化的气象因素以及地形因素的栅格数据产品，再对所有数据进行预处理；

步骤2：对所有预处理后的数据进行划分，得到若干个大小为N×N的正方形区域用于后续建模；

步骤3：针对每个划分得到的正方形区域，以雪水当量为因变量，影响雪水当量变化的因素为自变量，分别进行普通最小二乘线性回归模型的构建：

步骤4：判断每个分块的普通最小二乘线性回归模型残差是否具有显著的空间自相关性，若存在显著的空间自相关性，则执行步骤5，否则，将步骤3计算得到的普通最小二乘线性回归模型的参数作为该分块的最终模型参数；

步骤5：构建分块残差模型；对分块的普通最小二乘线性回归模型残差构建空间误差模型或空间滞后模型；

步骤6：将分块残差模型结果并入到步骤3构建的普通最小二乘线性回归模型中，得到最终的空间回归模型；若该分块未执行步骤5，则将步骤3构建的普通最小二乘线性回归模型作为该分块的最终模型；

步骤7：对模型进行评价和分析。

进一步地，步骤1中，获取的影响雪水当量变化的因素包括但不限于降水、空气温度、风速、植被覆盖、海拔、坡度、坡向等；预处理包括统一数据投影坐标系、地理配准、统一数据空间分辨率、统一数据空间范围、异常值处理、数据标准化。

进一步地，步骤3中构件的普通最小二乘线性回归模型为：

y＝α+β

式中，y表示该分块的雪水当量观测值，α表示该分块普通最小二乘线性回归模型的常数项，x

进一步地，步骤4中，通过计算分块模型残差的莫兰指数的p值(p值为当原假设为真时出现偏离原假设值的观测值以及比观测值更极端的值的概率，该值越小说明结果越显著)判断分块模型残差是否具有显著的空间自相关性，若p值大于阈值，则说明该分块的模型残差仍具有显著的空间自相关，进行步骤5和步骤6；反之则说明分块的模型残差不具有显著的空间自相关，则将步骤3计算得到的普通最小二乘线性回归模型的参数作为该分块的最终模型参数。

进一步地，步骤5具体包括：

步骤5.1：按照分块模型残差的地理分布构建空间权重矩阵W；

步骤5.2：进行拉格朗日乘子检验以及稳健性检验，分别计算LM-Error(空间误差最大似然)、Robust LM-Error(稳健性空间误差最大似然)、LM-Lag(空间滞后最大似然)和Robust LM-Lag(稳健性空间滞后最大似然)统计量，若LM-Error的p值大于阈值则构建空间误差模型，若LM-Lag的p值大于阈值则构建空间滞后模型，若LM-Error和LM-Lag的p值均大于阈值，则对比Robust LM-Error和Robust LM-Lag的p值，若Robust LM-Error的p值大于Robust LM-Lag的p值，则构建空间误差模型，反之则构建空间滞后模型；

步骤5.3：使用步骤5.1生成的权重矩阵W，按照步骤5.2的结果构建空间误差模型或空间滞后模型，得到的空间滞后模型如下：

ε＝ρWε+β′

式中，ε表示该分块的普通最小二乘线性回归模型残差，ρ表示空间效应系数，W表示该分块残差的空间权重矩阵，β′

得到的空间误差模型如下：

ε＝λWδ+β′

式中，ε表示该分块的普通最小二乘线性回归模型残差，λ表示空间误差相关系数，W表示该分块残差的空间权重矩阵，δ表示ε的普通最小二乘线性回归模型的残差，β′

进一步地，步骤6中，将分块残差模型结果并入到步骤3构建的普通最小二乘线性回归模型中，得到最终的空间回归模型表达式为：

y＝α+λWε+(β

式中，(β

进一步地，步骤7具体包括：

步骤7.1：模型精度评价；计算最终空间回归模型的R

式中，y

式中，y′

步骤7.2：检验空间自相关；计算空间回归模型残差的莫兰指数及其p值，用以衡量模型残差的空间自相关性，其p值大于阈值说明模型残差已经不存在显著的空间自相关，证明模型已经消除了空间自相关的影响；其中，

MC＝e

式中，MC表示残差莫兰指数，e表示残差向量；

步骤7.3：交叉验证；交叉验证采用m折交叉验证的方法计算雪水当量空间回归模型的RMSE，将雪水当量及其影响因素数据集随机分成m份，将其中的m-1份数据作为训练集，剩余的1份数据则为验证集；每一份数据都需要作为1次验证集，每一份数据都做过验证集后，计算这m次模型结果的RMSE均值作为交叉验证的结果，结果越小说明模型的鲁棒性越高。

与现有技术相比，本发明的有益效果为：

1)本发明通过对栅格数据产品进行划分，并对每个分块构建普通线性回归模型，再针对分块残差构建空间计量模型，从而消除了普通线性回归模型中残差空间自相关的影响，提高了模型的精度，探究不同因素对雪水当量的影响；

2)本发明考虑了地理数据空间自相关的影响，使得空间统计模型能够应用于大尺寸的积雪遥感数据产品，一定程度上缓解了在探究雪水当量影响因素对雪水当量的影响时进行回归建模的计算压力，从而提高了模型精度，更为准确地描述影响因素与雪水当量之间的定量关系。

附图说明

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

具体实施方式

下面将结合本发明实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面结合具体实施例对本发明作进一步说明，但不作为本发明的限定。

雪水当量的分布受空间效应的影响，传统回归模型往往会忽略数据的空间自相关性以及空间异质性等空间效应，导致模型精度较低。考虑空间效应的空间计量模型虽然可以解决上述问题，却难以直接在大尺寸栅格数据集上进行应用，计算耗时长且效率低下。本发明通过对数据集进行分块并针对分块的普通最小二乘线性回归模型的残差构建空间计量模型，消除空间自相关的影响，提高模型精度，从而便于探究影响因素对于雪水当量的影响。

参见图1，本发明提供的一种基于雪水当量栅格数据的空间回归建模方法，，包括如下步骤：

步骤1：数据获取与预处理。包括以下子步骤：

步骤1.1：在科学数据银行获取高亚洲地区雪水当量数据集，并从中国科学院计算机网络信息中心地理空间数据云平台、美国国家航空航天局官网、国家青藏高原科学数据中心等平台获取影响雪水当量变化的气象因素以及地形因素的栅格数据产品，包括降水、空气温度、风速、植被覆盖、海拔、坡度、坡向等；

步骤1.2：对照数据本身的分类描述对数据进行筛选和清洗，剔除数据集中的无效值或使用插值方法对数据进行插补；

步骤1.3：将所有数据进行投影转换，并以雪水当量数据为基准对其余自变量进行地理配准。针对空间分辨率不一致的数据，对数据进行重采样，统一数据的空间分辨率。针对空间范围不一致的数据，对数据进行裁剪，统一数据的空间范围。针对时间分辨率不一致的数据，使用地图代数将数据统一处理为月值数据，统一数据的时间分辨率；

步骤1.4：采用z-score标准化对自变量数据进行标准化处理，消除因量纲不同或部分数值过大产生的误差。

步骤2：对预处理后的数据进行分块处理。根据裁剪后的栅格尺寸和以及实验精度的要求，确定所需划分的子集大小，将原始栅格精确分割为若干个N×N的正方形区域，便于后续建模。

步骤3：将步骤2划分得到的雪水当量数据和影响因素数据按照地理位置一一对应，针对分块形成对应的数据表。对每一个分块的数据表，以雪水当量为因变量，影响因素为自变量，分别构建普通最小二乘线性回归模型，得到每一个分块的模型残差ε。其中，残差的表达式为：

式中，y为雪水当量的观测值，

步骤4：通过计算步骤3得到的分块残差的莫兰指数判断该残差是否具有显著的空间自相关性。若分块的残差莫兰指数的p值小于阈值(p值为当原假设为真时出现偏离原假设值的观测值以及比观测值更极端的值的概率，该值越小说明结果越显著)，在本实施例中，阈值为0.05，则说明分块的模型残差具有显著的空间自相关性。针对具有显著空间自相关的模型残差将进一步构建空间计量模型，执行步骤5。若分块的残差莫兰指数的p值大于阈值，在本实施例中即当p值大于0.05，则将步骤3计算得到的模型结果作为该分块的最终模型结果；其中，残差莫兰指数的计算公式为：

MC＝e

式中，MC表示残差莫兰指数，e表示残差向量。

步骤5：针对步骤4中模型残差仍具有显著空间自相关性的分块构建空间计量模型，包括以下步骤：

步骤5.1：按照该分块残差的地理分布构建空间权重矩阵W，可以依据分块中残差的邻接关系构建0-1邻接矩阵代替，邻接关系包含只考虑公共边的Rook邻接、只考虑公共点的Bishop邻接和同时考虑公共边和公共点的Queen邻接，两单元邻接时，矩阵对应元素值为1，反之为0；也可以根据距离阈值构建空间邻接矩阵代替。栅格数据通常考虑使用Queen邻接构建相应的空间邻接矩阵；

步骤5.2：构建以残差为因变量，影响因素为自变量的普通最小二乘线性回归模型，分别计算LM-Error(空间误差最大似然)、Robust LM-Error(稳健性空间误差最大似然)、LM-Lag(空间滞后最大似然)和Robust LM-Lag(稳健性空间滞后最大似然)统计量，若LM-Error的p值大于阈值，在本实施例中，阈值为0.05，即若LM-Error的p值大于0.05，则构建空间误差模型，若LM-Lag的p值大于0.05则构建空间滞后模型，若LM-Error和LM-Lag的p值均大于0.05，则对比Robust LM-Error和Robust LM-Lag的p值，若Robust LM-Error的p值大于Robust LM-Lag的p值，则构建空间误差模型，反之则构建空间滞后模型；

步骤5.3：根据步骤5.2的结果，构建以残差为因变量，影响因素为自变量的空间误差模型或空间滞后模型，得到的空间滞后模型如下：

ε＝ρWε+β′

式中，ε表示该分块的普通最小二乘线性回归模型残差，ρ表示空间效应系数，W表示该分块残差的空间权重矩阵，β′

得到的空间误差模型如下：

ε＝λWδ+β′

式中，ε表示该分块的普通最小二乘线性回归模型残差，λ表示空间误差相关系数，W表示该分块残差的空间权重矩阵，δ表示ε的普通最小二乘线性回归模型的残差，β′

步骤6：将步骤5生成的空间计量模型结果加入到步骤3构建的普通最小二乘线性回归模型中，得到该分块的空间回归模型。若该分块未执行步骤5，则将步骤3构建的普通最小二乘线性回归模型作为该分块的最终模型。

步骤7：模型评价和分析，包括以下步骤：

步骤7.1：计算最终空间回归模型的R

式中，y

式中，y′

步骤7.2：计算空间回归模型残差的莫兰指数及其p值，用以衡量模型残差的空间自相关性。其p值大于阈值说明模型残差已经不存在显著的空间自相关，证明模型已经消除了空间自相关的影响；其中，

式中，MC表示残差莫兰指数，n是要素总数，x

步骤7.3：交叉验证；采用交叉验证的方法(在本实施例中，采用10折交叉验证)计算雪水当量空间回归模型的RMSE，将雪水当量及其影响因素数据集随机分成10份，将其中的9份数据作为训练集，剩余的1份数据则为验证集。每一份数据都需要作为1次验证集，每一份数据都做过验证集后，计算这10次模型结果的RMSE均值作为交叉验证的结果。结果越小说明模型的鲁棒性越高。

以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，均应当包含在本发明的保护范围内。