基于支持向量机二元数据同化的地下水水位变化中长期预报方法与流程

本发明属于水文水资源学科地下水预报领域，尤其涉及一种基于支持向量机二元数据同化的地下水水位变化中长期预报方法。

背景技术：

地下水资源作为全球最大的淡水资源[shiklomanovandrodda,2003]，在全球水循环中占据着重要的位置[tayloretal.,2013]。地下水作为生活饮用水的主要来源，维系着生态环境的可持续发展[ipcc,2014]。随着人口增加和气候变化，地下水的需求量急剧增加[konikow,2011；wadaetal.,2011；wadaetal.,2012]。随着地下水的开采，地下水含水层面临着紧缩的危机[sun,2013]。了解认知地下水的时空分布，对地下水进行中长期预报将有利于地下水资源的合理分配，为应对气候变化，保障水资源的可持续利用提供决策依据。

地下水水位动态研究和中长期预报在地下水水资源管理应用、环境污染检测管理等方面有着广泛的应用。由于地下水埋藏于地下，加上岩层结构的复杂和影响因素的多变，给地下水动态研究带来了困难。基于地下水运动的物理机理构建模型如水均衡模型，物理概念清晰，但有些参数在实际应用中包含有大量的物理参数，有些参数在实际应用中难以准确的获得。随着计算机技术的迅速发展，机器学习算法也开始用于地下水水位动态预报。与传统的物理模型相比，机器学习在预报的时候不需要研究区域的具体的地质结构等相关参数。支持向量机技术在水文水资源学科预报中应用较广[gillandmckee,2007；gilletal.,2006；dliuetal.,2010；dliuetal.,2016；yliuetal.,2012；yoonetal.,2011；dliuetal.,2017]。数据同化技术可以有效的提高模型的模拟效率。目前，基于支持向量机的二元数据同化技术在地下水水位变化的中长期预报中暂无相关研究。

技术实现要素：

发明目的：为预报地下水水位动态变化，解决现有技术的不足，本发明首次将基于支持向量机的二元数据同化技术应用于地下水水位变化的中长期预报中，提出一种基于支持向量机二元数据同化的地下水水位变化中长期预报方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于支持向量机二元数据同化的地下水水位变化中长期预报方法，包括如下步骤：

步骤1，获取n个时刻研究区域有效地下水水位和气象要素观测数据集，根据地下水水位变化与气象要素的相关性选取支持向量机的气象驱动要素；

步骤2，采用支持向量机非线性回归方法，构建地下水水位变化预报模型；模型的输入为某一时刻的驱动要素，输出为该时刻的地下水水位变化值；

步骤3，将数据集样本中前t时刻的驱动要素和地下水水位变化值作为训练集，训练预报模型，获得训练好的地下水水位变化预报模型；

步骤4，当预报时刻t+m存在地下水水位观测数据，更新训练集，采用数据同化技术同化校正该集合，采用校正的训练集重新训练预报模型，获取一元数据同化预报模型；

步骤5，将校正后的t+m时刻的驱动要素，输入至一元数据同化预报模型，获取该预报时刻的地下水水位变化一元数据同化预报值；

步骤6，对预报时刻t+m的地下水水位变化一元数据同化预报值再次进行数据同化校正，获取二元数据同化预报值，即为该时刻最终的地下水水位变化预报值。

进一步的，所述步骤1，具体包括：

(1)获取地下水水位和气象要素的观测数据时间序列，对观测数据进行有效性筛选，去除异常数据和缺失数据，选用有效时段的数据序列；

(2)将选定时段内的数据序列转化为逐月数据序列，进行去趋势处理，获取逐月地下水水位变化时间序列；

(3)采用皮尔逊相关分析方法分析地下水水位变化与气象要素在不同时段的相关性，根据相关性选取支持向量机的驱动要素。

进一步的，所述步骤2，具体包括：

获取svm模型识别与回归软件包ls-svm；根据驱动要素调用ls-svm软件包中的定向结构函数initlssvm，在该函数中设置核函数，此处采用径向基核函数；调用trainlssvm命令，构建地下水水位变化预报模型。

进一步的，所述步骤3，根据构建的地下水水位变化预报模型，调用simlssvm命令，输入预报时刻的驱动要素，获取该预报时刻的地下水水位变化预报值。

进一步的，所述步骤4，具体包括：

若在预报的t+m时刻存在地下水水位观测数据，则将训练集向前更新m-1步，即，移除1至m-1时段的驱动要素和地下水水位变化值，将t+1至t+m-1时段的驱动要素和地下水水位变化的预报值添加至训练集，构建新的训练集；

采用数据同化技术，融合t+m时刻地下水水位观测数据，同化校正训练集；使用更新校正的训练集，创建支持向量机，即为支持向量机数据同化模型svm-da(supportvectormachine-dataassimilation)；将校正后的t+m时刻的驱动要素输入svm-da模型，输出该时刻的地下水水位预报值，即为一元数据同化预报值。

进一步的，所述数据同化技术采用集合卡尔曼滤波，或者集合粒子滤波，或者无迹粒子滤波。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

本发明采用支持向量机技术，能够同时控制模型的复杂度和经验误差，具有很强的泛化能力，在解决小样本、非线性及高维数问题方面表现出优异的特性，可以有效避免维数灾难和过学习问题；与传统的物理模型相比，机器学习在预报的时候不需要输入研究区域的具体的地质结构等相关参数。

本发明采用的数据同化技术通过集合估计真实统计值的方法提高了估计精度，同时也大幅度地减少了计算量，可以有效的提高模型的模拟效率。本发明提出的一种基于支持向量机二元数据同化的地下水水位变化中长期预报方法，基于机器学习理论，采用数据同化技术融合大数据预报地下水水位变化，可直接服务于地下水水资源的中长期预报及优化配置，应用前景广阔。

附图说明

图1是支持向量机svm-da二元数据同化模型流程图；

图2是选取的各站点的逐月地下水水位变化及气象驱动要素时间序列；

图3是选取的各站点的地下水水位变化在次一月的预报值与观测值的逐月时间序列；

图4是选取的各站点的地下水水位变化在次二月的预报值与观测值的逐月时间序列；

图5是选取的各站点的地下水水位变化在次三月的预报值与观测值的逐月时间序列。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本实施例现有2007年1月1日至2016年12月31日的逐日地下水水位及气温、降水、太阳辐射和地表温度气象要素的观测数据序列。地下水水位观测站站点所在的经纬度和观测井井位深度如表1所示。依据图1所示的本发明方法，预报各观测站点的地下水水位在次一月至次三月的动态变化。过程如下：

(1)根据研究区域内的逐日地下水水位和气温、降水、太阳辐射、地表温度数据序列，对所述数据进行有效性筛选，去除地下水水位和气象要素数据序列中的异常数据和缺失数据，选用有效时段的数据序列。

表1

(2)将选定时段内的逐日数据序列通过累加求平均方法转化为逐月数据序列。根据逐月地下水水位数据序列，采用matlab工具箱中的hpfilter函数计算逐月地下水水位变化数据序列。

图2是选取的各站点s1至s4在2007年1月至2016年12月期间的逐月地下水水位变化(groundwaterlevelchanges(gw))及气温(2mairtemperature(t2m))，降水(precipitation(p))，太阳辐射(solarradiation(sr))和地表温度(groundtemperature(gt))数据序列。由图2可见，地下水水位变化在各站点均有年季周期性变化的特征，其中以夏秋季变化幅度较大，冬春季变化较小。与此同时，气温，太阳辐射和地表温度的年季变化也呈现周期性，四季较为分明，夏季量值最大，冬季量值最小。降水的年季变化较地下水水位变化及气温、太阳辐射和地表温度较弱，量值在夏季较大。

(3)采用皮尔逊相关分析方法分析地下水水位变化与气温、降水、太阳辐射、地表温度在不同时段的相关性，根据相关关系强弱选取支持向量机的气象驱动要素。

地下水水位变化与各气象要素的相关关系计算结果如表2所示，其中s1至s4表示选取的4个站点，lag0表示当前月地下水水位变化与当前月气温，降水，太阳辐射和地表温度的相关系数。lag1至lag3表示当前月地下水水位变化与前1个月至前3个月的气温，降水，太阳辐射和地表温度的相关系数。由表2可得，地下水水位变化与当前月及前1个月至前3个月的气温的相关系数在各站点s1至s4上均呈正相关，且在大部分站点上呈显著正相关。通过对比发现，地下水水位变化与前1个月和前2个月的气温的相关系数量值较气温在当前月及前3个月的略大。由此推断地下水水位变化受前1个月至前2个月的气温的影响较大。同理，地下水水位变化与当前月及前1个月至前3个月的地表温度的相关系数在各站点s1至s4上也呈正相关且在大部分站点上呈显著正相关。其中地下水水位变化与前1个月至前2个月的地表温度的相关系数的量值较大。地下水水位变化与当前月及前1个月至前3个月的太阳辐射的相关关系在各站点上呈正相关性。与气温和地表温度一致，地下水水位变化与前1个月至前3个月的太阳辐射的相关关系的量值较大，在s1至s3上呈显著正相关。与此同时，地下水水位变化与当前月及前1个月至前3个月的降水的相关系数在各站点上较弱，且大部分呈负相关。由相关关系强弱判断地下水水位变化受前1个月至2个月的气温、地表温度和太阳辐射的影响较大，而受降水的影响较小。

表2

(4)计算地下水水位变化在月季尺度的自相关系数，计算结果如表3所示。由表3可得，各站点的地下水水位变化在月尺度的自相关性显著正相关，随着间隔时间的增长，自相关性减弱。

表3

(5)选取支持向量机的驱动要素；采用支持向量机非线性回归方法，建立地下水水位变化预报模型svm。地下水水位变化预报模型参数设置如表4所示。

由表4可得，t+1，t+2，t+3为地下水水位变化的预报时刻，分别为对次一月(t+1)至次三月(t+3)地下水水位变化的预报。其中t+1预报时刻的输入的驱动要素包括t时刻月降水总量(precipitation(p)，单位：mm)，月平均太阳辐射(solarradiation(sr)，单位：wj/m²)；t-1和t时刻月平均气温(2mairtemperature(t2m)，单位：摄氏度)及月平均地表温度(surfacegroundtemperature(gt)，单位：摄氏度)。t+2预报时刻的输入的驱动要素包括t时刻月降水总量(p)，月平均太阳辐射(sr)和t+1时刻的地下水水位预报值(groundwaterlevelchanges(gw)，单位：mm)；t-1和t时刻月平均气温(t2m)及月平均地表温度(gt)。t+3预报时刻的输入的驱动要素包括t时刻月降水总量(p)，月平均太阳辐射(sr)；t+1和t+2时刻的地下水水位预报值(groundwaterlevelchanges(gw))；t-1和t时刻月平均气温(t2m)及月平均地表温度(gt)。输出要素为逐月平均的地下水水位变化(单位：mm)。

其中，2007年1月至2011年12月逐月数据序列作为模型训练期，通过模拟训练，得到支持向量机预报模型。2012年1月至2016年12月数据序列作为模型预报期，根据构建的支持向量机预报模型，预报次一月至次三月的地下水水位变化。

表4

(7)若在预报的t+m时刻存在地下水水位观测数据，则将训练集向前更新m-1步，即，移除1至m-1时段的驱动要素和地下水水位变化值，将t+1至t+m-1时段的驱动要素和地下水水位变化的预报值添加至训练集，构建新的训练集。

采用enkf数据同化技术，融合t+m时刻地下水水位观测数据，同化校正训练集。使用更新校正的训练集，创建支持向量机，即为支持向量机数据同化模型svm-da(supportvectormachine-dataassimilation)。将校正后的t+m时刻的驱动要素输入一元数据同化模型svm-da，输出该时刻的地下水水位变化，即为一元数据同化预报值。

(8)再次通过enkf数据同化技术融合地下水水位观测数据，更新校正步骤(7)中输出的地下水水位变化的一元数据同化预报值，获得地下水水位变化二元数据同化预报值，即为t+m时刻最终的同化校正预报值。

图3是2012年1月至2016年12月次一月地下水水位变化的预报值与观测值的逐月时间序列对比图。s1至s4为选取的研究站点。其中(a)为观测的地下水水位变化，(b)-(d)分别为svm初始预报值，svm-da一元数据同化预报和二元数据同化预报值。如图3所示，与观测的地下水水位变化时间序列相比，svm的初始预报值和svm-da的一元和二元数据同化预报值都较好的拟合了次一月地下水水位变化。其中，svm-da二元数据同化校正的预报值拟合效果最好，与观测的地下水水位变化的偏差最小。而svm的初始预报值和svm-da的一元数据同化预报值在一些预报时刻与观测的地下水水位变化的偏差较大。图3表明svm-da的二元数据同化预报值有效提高次一月地下水水位预报的准确性。

图4是2012年1月至2016年12月次二月地下水水位变化的预报值与观测值的逐月时间序列对比图。s1至s4为选取的研究站点。其中(a)为观测的地下水水位变化，(b)-(d)分别为svm初始预报值，svm-da一元数据同化预报和二元数据同化预报值。如图4所示，构建的svm和svm-da模型对次二月地下水水位变化的预报在大部分时段效果较好：预报值与观测值的变化趋势一致，偏差较小。其中，svm-da的二元数据同化预报值的效果最优，预报的结果与观测值拟合效果最好，偏差最小。

图5是2012年1月至2016年12月次三月地下水水位变化的预报值与观测值的逐月时间序列对比图。s1至s4为选取的研究站点。其中(a)为观测的地下水水位变化，(b)-(d)分别为svm初始预报值，svm-da一元数据同化预报和二元数据同化预报值。如图5所示，构建的svm和svm-da模型对次三月地下水水位变化的预报也较好，各预报值与观测值拟合性较好，偏差较小。其中，svm-da的二元数据同化预报值的效果最优。

(9)采用相关系数(r)，平均绝对误差(mae)，和均方根误差(rmse)，纳什系数(nse)四种统计指标进一步检验svm及svm-da对次月季地下水水位变化的预报能力。各统计指标的计算公式如下：

其中，yi是第i步观测的地下水水位变化，i从1到n；n为预报的时段长；yi^*是预报的第i步地下水水位变化；和分别是yi和yi^*在预报期的均值。

通常情况下，rmse和mae量值越小且r和nse量值越大，预报效果越好。一般认为，当r量值超过0.5时，预报效果达到满意水平，若量值超过0.75，则预报效果非常好；若mae和rmse量值低于0.6时，则预报效果非常好；若nse大于0，则预报效果非常好，越接近1效果越好。

各站点的预报结果如表5所示。由表5可得，svm及svm-da对各站点的次一月至次三月的地下水水位变化的预报效果均较好。各预报值的相关系数(r)均为显著正相关。与此同时，平均绝对误差(mae)和均方根误差(rmse)的量值较小，且在大部分站点及各预报时刻量值小于1。而纳什系数(nse)大于零，且在大部分站点及各预报时刻量值大于0.5。在svm初始预报和svm-da一元数据同化预报和二元数据同化预报值中，svm-da二元数据同化的预报值最接近观测值，预报效果最好，准确性最高。在各站点的次一月至次三月的预报中，svm-da二元数据同化的预报值与观测值的相关系数(r)量值最大，平均绝对误差(mae)和均方根误差(rmse)的量值最小，而纳什系数(nse)量值最大。

表5

综上所述，本发明提出的一种基于支持向量机二元数据同化的地下水水位变化中长期预报方法基于机器学习理论，根据研究区域内观测的地下水水位和气象要素时间序列，采用支持向量机非线性回归方法和数据同化技术，构建支持向量机地下水水位变化预报模型，适应于对次月季尺度的地下水水位变化进行实时预报。本发明可直接服务于地下水水资源的中长期预报及优化配置，应用前景广阔。

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