一种基于零膨胀时空回归模型的雪水当量建模和预测方法

1.本发明属于地学统计与空间分析应用领域，特别涉及一种基于零膨胀时空回 归模型的雪水当量建模和预测方法。


背景技术：

2.积雪是地球表面水循环的重要环节，随季节发生周期性变化，是能够反映全 球气候变化的敏感指标。雪水当量被定义为积雪完全融化后所得到的水形成水层 的垂直深度，是能够将积雪以统一量纲标准进行度量的重要参数。积雪在可见光 波段具有高反射性，在部分中红外波段具有更强的吸收性，雪的光谱属性以及其 电磁辐射的交互方式，以及积雪的分布和动态特征，使得遥感成为监测积雪分布 和变化的有效手段。被动微波能够穿透云层、积雪和大气，能够全天候、全天时 地工作，在估算积雪参数上具有很大优势。目前覆盖北半球的基于格网的雪水当 量产品，包括重分析类产品、融合雪深观测数据的被动微波遥感产品和被动微波 直接获取产品。其中，融合了被动微波遥感和地面站点观测数据的globsnow数 据集具有较低的整体误差和偏差，成为估算雪水当量的最佳产品。2020年发布 的globsnow 3.0产品延续了上一代产品的数据同化算法，融合了被动微波遥感 和气象站点观测资料。
3.积雪的空间分布会随季节发生周期性变化，不仅受到气候、下垫面等地表环 境因素影响，同时还与环境相互作用从而影响地表环境因素变化，因此分析雪水 当量与环境因素的交互效应、构建雪水当量与环境因素之间的定量模型，是研究 雪水当量与环境因素的关系和预测雪水当量未来变化的基础。根据地理学第一定 律，雪水当量和气候、下垫面、热辐射等环境因素的空间中连续分布，具有空间 自相关性，不满足传统回归模型独立同分布的前提假设。以空间自回归模型、地 统计模型、地理加权回归模型、空间滤值模型等为代表的空间回归模型，考虑了 观测对象的空间关系，将空间效应的影响加入到回归模型中，考虑了空间自相关 性导致的方差膨胀效应和回归系数的偏移效应，从而降低空间效应对模型的影响、 提高模型精度。
4.从时间变化尺度上看，雪水当量不仅与同期环境因素有交互作用，也受到过 去时期雪水当量甚至各种环境因素的影响，具有时间自相关性。在建模过程中， 这种时间滞后效应需要通过构建时间滞后变量来代表，使静态分析的问题有可能 成为动态分析。以arima为代表的传统的时间序列模型能够简单清晰的定量描 述时间滞后效应，但需要满足平稳性的前提要求；以lstm、rnn等为代表的 机器学习模型虽然能够提高时间序列模拟和预测的精度，但其复杂的内在原理降 低了模型的可解释性，同时，这些时间序列模型仅考虑时间效应的影响，未同时 考虑空间效应的影响。以star为代表的时空自回归模型，能够考虑时间效应和 空间效应对模型的影响，但需要构建时空权重矩阵，对于基于格网的雪水当量和 环境因素变量，由于其数据量非常庞大，会面临严重的计算瓶颈。
5.积雪的时空分布随季节变化呈现出明显的变化，在很多纬度和海拔较低的地 区会有很长时间在较大空间范围内没有被积雪覆盖，在雪水当量的时空数据中表 现为大量
的零观测值，也就是零膨胀效应。对于连续性变量，以hurdle模型、 two-part模型为代表的零膨胀模型能够将零值和非零值的分别考虑，混合模型甚 至可以将零值的来源进行细致的区分，从而避免了大量零观测值造成数据分布不 满足模型假设的前提。但对于受大量零观测值影响雪水当量与环境因素时空数据 的建模分析，尚未有相应的建模方法。
6.综上所述，对于雪水当量与环境因素建模过程中，缺乏对时间效应、空间效 应以及零膨胀效应的同时考虑，因此，亟待提供一种基于零膨胀时空回归模型的 雪水当量建模方法，为雪水当量与环境因素间时空定量关系描述、及雪水当量变 化预测提供基础支撑。


技术实现要素：

7.本发明的目的在于提供一种基于零膨胀时空回归模型的雪水当量建模和预 测方法，从而同时解决在雪水当量与环境因素建模过程中空间效应、时间效应及 零膨胀效应对模型造成的影响，更加准确的定量描述雪水当量与环境因素的交互 效应，提高雪水当量预测精度。
8.本发明所采用的技术方案包括以下步骤：
9.步骤1：数据获取与处理；获取雪水当量栅格数据产品及与积雪聚积和消融 相关的环境因素遥感数据产品，将所有的栅格数据产品进行投影转换成与雪水当 量相同的投影坐标系，将投影后的栅格数据重采样成与雪水当量栅格数据相同的 空间分辨率；
10.步骤2：构建空间邻接列表；根据栅格单元的空间邻接关系，构建空间邻接 列表；在空间邻接列表中，每个栅格单元都对应的空间邻接关系表，表中是对应 的邻接栅格单元的id号和权重值；
11.步骤3：时间序列数据重组与拼接；为每一期的雪水当量生成其时间滞后变 量，将每一期的雪水当量和环境因素数据按照行或列的形式重新排列，构建时间 序列数据集；同时，将空间邻接列表按照相应的时间序列拼接，构建时间序列空 间邻接列表；
12.步骤4：生成雪水当量二元分类变量；为重组后的时间序列数据集生成雪水 当量二元分类变量π以区分零值和非零值，新生成的二元分类变量π取值为0或1， 其中0代表雪水当量为零的单元，π取值为1代表雪水当量不为零的单元；
13.步骤5：二元部分时空建模；以步骤4中雪水当量二元分类变量为因变量， 环境因素数据及时间滞后变量作为自变量，采用广义线性回归模型对二元部分雪 水当量进行时空建模；
14.步骤6：连续部分时空建模；提取出步骤4中雪水当量二元分类变量的值为 1的记录并存储为非负子数据集，以子数据集中的雪水当量为因变量，环境因素 数据及时间滞后变量作为自变量，进行连续部分雪水当量的时空建模；
15.步骤7：构建零膨胀时空回归模型；将步骤5中二元部分得到的概率值和步 骤6连续部分得到的预测值进行匹配相乘，得到零膨胀时空回归模型拟合的雪水 当量值；
16.步骤8：预测分析；对步骤7中的模型进行精度评价，当模型精度较高时， 将最近的环境因素遥感数据产品预处理后，作为自变量x输入到步骤7构建的 零膨胀时空回归模型中，与步骤7中的建模参数进行运算，得到最近的雪水当量 的预测值。
17.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤1 中，雪水当量栅格数据产品通常采用被动微波遥感和地面观测的数据融合产品， 相关的环境
因素遥感数据包括但不限于海拔、纬度、降水、温度、地表热通量、 水蒸气含量、云水含量、植被覆盖、风速。
18.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤2 中，对于栅格数据结构，常见的邻接关系有：rock邻接，bishop邻接，queen 邻接等；除了邻接关系，也可采用基于距离的空间权重方式。
19.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤3 中，将该像元当前雪水当量对应的前期雪水当量整合进数据集中，再将所有时期 的数据按照时间序列拼接起来，构建时间序列数据集；同时按照整合的时间序列 数据的新id将空间邻接列表按照相应的时间序列拼接，构建时间序列空间邻接 列表。
20.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤5 中，二元部分时间效应由时间滞后变量π
t-1
及相应的线性系数γ
t
来代表，二元部 分空间效应fb可采用广义线性空间回归模型，其中x
b，k，t
和γ
k，t
分别代表二元部分 的因变量和对应的系数，w为步骤3拼接得到的空间邻接列表生成的空间邻接矩 阵；二元部分模型可表示如下：
21.π
t
～γ
t
π
t-1
+fb(x
b，k，t
，w，γ
k，t
)。
22.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤6 中，连续部分时间效应由时间滞后变量y
t-1
及相应的线性系数β
t
来代表，连续部 分空间效应fc可采用空间回归模型，其中x
c，k，t
和β
k，t
分别代表连续部分的因变量 和对应的系数，wc为提取出的雪水当量二元分类变量的值为1的记录对应的空 间邻接列表生成的空间邻接矩阵；连续部分模型可表示如下：
23.y
t
～β
tyt-1
+fc(x
c，k，t
，wc，β
k，t
)。
24.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤7 中，零膨胀模型拟合值计算可表示如下：
25.e(swe
t
)＝y
t
×
π
t
。
26.在上述的一种基于零膨胀时空回归模型的雪水当量建模和预测方法，步骤8 中，雪水当量预测步骤包括：
27.(1).模型精度评价，选择统计学常用的均方根误差、拟合优度、残差的莫兰 指数对步骤7中的模型精度进行评价，若拟合优度较高且均方根误差和残差的莫 兰指数小，则用步骤7中的模型进行后续预测；若模型精度低，则返回步骤1 重新选择环境因素；
28.(2).采集预测所需环境数据，采集步骤1中最终确定的环境因素最近遥感产 品并按照步骤1中的预处理方式进行数据预处理；
29.(3).雪水当量预测，将数据输入到步骤7构建的模型中进行运算，得到最近 雪水当量的预测值。
30.本发明所提供的一种基于零膨胀时空回归模型的雪水当量建模和预测方法， 符合雪水当量和环境因素的时空分布特征，同时考虑了时间滞后效应、空间自相 关效应和零膨胀效应，从而更准确地构建雪水当量时空分布与环境因素间的定量 关系，提高雪水当量建模和预测精度。
附图说明
31.图1为本发明技术方案的流程图。
32.图2为本发明实施例数据处理流程图。
33.图3为本发明实施例步骤2中空间邻接关系表示意图。
34.具体实施方法
35.为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对 本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解 释本发明，并不用于限定本发明。
36.本发明要解决的核心问题是：积雪时空分布随季节发生周期性变化，与气候、 下垫面等地表环境因素间具有很高的关联交互作用。在积雪与地表环境因素时空 建模中，其空间分布、积雪时间滞后效应以及大量的零观测值都不满足现有模型 的前提假设。本发明通过构建雪水当量与环境因素的零膨胀时空回归模型，消除 零膨胀效应、空间自相关效应及时间滞后效应对建模的影响，从而提高雪水当量 建模和预测的精度，为雪水当量预测和全球气候变化趋势分析提供重要参考依据。
37.参见图1，本发明提供的一种基于零膨胀时空回归模型的雪水当量建模和预 测方法，包括如下步骤：
38.步骤1：数据获取与处理。参见图2，包括以下步骤：
39.步骤1.1：获取globsnow3.0雪水当量栅格数据产品和环境因素遥感数据产 品，包括海拔、空气温度、地表温度(日间/夜间)、地表热通量、云层水含量、 空气水蒸气含量、降水量、植被覆盖、风速等。其中部分产品可以进一步处理得 到其他关联的环境因子，海拔高程dem可以计算得到坡度和坡向变量，温度可 以通过统计得到每月最低温度低于0℃天数。针对北半球地区，选择ease网格 北极方位角等面积投影作为投影坐标系，将以上变量全部进行投影转换。
40.步骤1.2：将所有投影变换后的数据集进行重采样，使之具有统一的空间分 辨率。对不同时间分辨率的数据，采用计算和、均值、最大/最小值等方式进行 处理统一成月值数据，使之具有相同的时间分辨率和空间分辨率。
41.步骤1.3：数据检查和探索性分析。对数据集中所有变量进行描述性统计分 析，对异常值、空值等进行处理，根据具体情况，可选择插值方法插补或直接剔 除。
42.步骤2：构建空间邻接列表。根据栅格单元的空间邻接关系，构建空间邻接 列表。在空间邻接列表中，每个栅格单元都对应的空间邻接关系表，表中是对应 的邻接栅格单元的id号和权重值。参见图3，对于栅格单元通常选择queen相 邻(8方向)方式构建空间邻接关系表，即每个栅格单元都和邻近的8个栅格单 元具有邻接关系，且权重值为1。
43.步骤3：时间序列数据重组与拼接。具体步骤如下：
44.步骤3.1：为每一期的雪水当量生成其时间滞后变量，此例中仅考虑月值数 据中前一个月的雪水当量值作为时间滞后变量，并将每一期的雪水当量和环境因 素数据以列形式重新排列，并按照倒序时间顺序构建时间序列数据集。
45.步骤3.2：将空间邻接列表按照相应的时间序列数据集重新生成id，并将具 有新id的空间邻接列表按照时间序列数据集的顺序对应拼接，构建时间序列空 间邻接列表。
46.步骤4：生成雪水当量二元分类变量。为重组后的时间序列数据集生成雪水 当量二元分类变量π以区分零值和非零值，新生成的二元分类变量π取值为0或1， 其中0代表雪水当量为零的单元，π取值为1代表雪水当量不为零的单元。
47.步骤5：二元部分时空建模。以步骤4中雪水当量二元分类变量为因变量， 环境因素数据及时间滞后变量作为自变量，采用广义线性回归模型对二元部分雪 水当量进行时空建模。二元部分模型可表示如下：
48.π
t
～γ
t
π
t-1
+fb(x
b，k，t
，w，γ
k，t
)
49.其中，二元部分时间效应由时间滞后变量π
t-1
及相应的线性系数γ
t
来代表， 二元部分空间效应fb可采用广义线性空间回归模型，其中x
b，k，t
和γ
k，t
分别代表二 元部分的因变量和对应的系数，w为步骤3拼接得到的空间邻接列表生成的空间 邻接矩阵。
50.步骤6：连续部分时空建模。提取出步骤4中雪水当量二元分类变量的值为 1的记录并存储为非负子数据集，以子数据集中的雪水当量为因变量，环境因素 数据及时间滞后变量作为自变量，进行连续部分雪水当量的时空建模。连续部分 模型可表示如下：
51.y
t
～β
tyt-1
+fc(x
c，k，t
，wc，β
k，t
)
52.其中，连续部分时间效应由时间滞后变量y
t-1
及相应的线性系数β
t
来代表， 连续部分空间效应fc可采用空间回归模型(空间滞后模型/空间误差模型/空间滤 值模型)，其中x
c，k，t
和β
k，t
分别代表连续部分的因变量和对应的系数，wc为提取 出的雪水当量二元分类变量的值为1的记录对应的空间邻接列表生成的空间邻 接矩阵。
53.步骤7：构建零膨胀时空回归模型。将步骤5中二元部分得到的概率值和步 骤6连续部分得到的预测值进行匹配相乘，得到零膨胀时空回归模型拟合的雪水 当量值，零膨胀模型可表示如下：
54.e(swe
t
)＝y
t
×
π
t
55.其中π
t
为二元部分建模计算得到的概率值，y
t
是连续部分计算得到的雪水当 量拟合值。由于y
t
连续部分的数据集是整个时间序列数据集的子集计算得到，需 要将对应未提取的部分用0值补齐后与二元部分相乘。
56.步骤8：模型评价与预测分析。具体步骤如下：
57.步骤8.1：模型精度评价。计算模型步骤7得到模型的均方根误差、拟合优 度、残差的莫兰指数，拟合优度r2越高、均方根误差rmse越低代表模型精度 越好，若拟合优度较高且均方根误差和残差的莫兰指数较小，则进行后续预测； 若模型精度较低，则返回步骤1重新选择环境因素，直到精度较好，通常精度需 要高于0.8。莫兰指数计算公式如下，其中e代表残差，是残差平均值，c
ij
是根 据步骤3中的时间序列空间邻接列表构建的观测点i和j之间的空间权重。
[0058][0059]
步骤8.2：采集预测所需环境数据，按照步骤1中最终确定的环境因素采集 其最近遥感产品，并对环境因素数据按照步骤1中的投影转换、重采样方式进行 预测数据预处理。
[0060]
步骤8.3：雪水当量预测。将步骤8.2中预处理后的环境数据输入到步骤7 构建的模型中进行运算，得到最近雪水当量的预测值。如有新的雪水当量产品或 观测值，可进一步计算预测误差。
[0061]
应当理解的是，上述针对本发明中较佳实施例的表述较为详细，但不能因此 认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下， 在本发明权利要求保护范围内，可以做出替换或变形，均落入本发明的保护范围 之内，本发明的请求
保护范围以所附权利要求书为准。