CN112257021A - 一种可选的克里金空间插值降雨量估算方法 - Google Patents

一种可选的克里金空间插值降雨量估算方法 Download PDF

Info

Publication number
CN112257021A
CN112257021A CN202011112171.5A CN202011112171A CN112257021A CN 112257021 A CN112257021 A CN 112257021A CN 202011112171 A CN202011112171 A CN 202011112171A CN 112257021 A CN112257021 A CN 112257021A
Authority
CN
China
Prior art keywords
scale
rainfall
function
kriging
variation function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202011112171.5A
Other languages
English (en)
Inventor
车磊
王海起
陈奇东
甄卫民
梁庆娜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Radio Wave Propagation CETC 22 Research Institute
Original Assignee
China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Institute of Radio Wave Propagation CETC 22 Research Institute filed Critical China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority to CN202011112171.5A priority Critical patent/CN112257021A/zh
Publication of CN112257021A publication Critical patent/CN112257021A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/14—Rainfall or precipitation gauges

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Hydrology & Water Resources (AREA)
  • Algebra (AREA)
  • Atmospheric Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种可选的克里金空间插值降雨量估算方法，包括如下步骤：步骤1，通过离散变异函数公式计算所有降雨量样本点对的实验变异函数值：步骤2，采用两尺度小波最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型：步骤3，建立克里金空间插值方程组，求解克里金权重系数；步骤4，根据克里金权重系数计算待估位置点的降雨量估计值；步骤5，定量化精度评价。本发明公开的估算方法，克服了普通克里金空间插值降雨量估算方法中传统理论变异函数模型形状固定且未考虑空间变化趋势及多尺度特征的不足，使得降雨量估算的结果符合本身的空间变化趋势且体现空间变化的多尺度特征。

Description

一种可选的克里金空间插值降雨量估算方法
技术领域
本发明属于空间插值降雨量估算领域，特别涉及该领域中的一种可选的克里金空间插值降雨量估算方法。
背景技术
降雨量作为区域水资源的输入项，对于区域水资源的评估至关重要。当前，对于区域降雨量的估算主要通过站点观测的手段，再利用空间插值模型估算区域范围内降雨量。这里提及的空间插值模型是根据已知观测样本点的降雨量真实信息来估计待估位置点的降雨量信息，其原理是通过已知观测点样本降雨量数据构建函数关系，综合空间位置关系以及空间相关性，从而估算其他任意点或任意分区的降雨量。
常用的空间插值法可以归结为两类，一类是确定性方法，另一类是空间统计方法。确定性方法中最具代表性的是反距离权重法和泰森多边形法。这种传统的概念性模型对降雨的空间变化一般采用均化处理方式，用以方便产流计算，但随之带来的问题就是难以客观反映降雨的空间分布，一定程度上忽略了降雨在空间上的随机性，也未能有效考虑不同观测点的降雨信息在空间上的相关性。空间统计方法中最具代表性的是克里金插值方法，该方法以区域化变量独有的随机性和结构性性质为基础，使得探索降雨量空间结构和空间变异规律变得有迹可循，同时以变异函数为基本工具对区域化变量的连续性、相关性、尺度性等要素进行空间描述。降雨量正是具有这种随机性与结构性双重特征的区域化变量，应用克里金空间插值实现对降雨量的估算其实质在于通过已知位置点的降雨量内插或外推的方式，对待估位置点降雨量的取值进行无偏、最优估计。
克里金空间插值降雨量估算过程中，插值模型的精度取决于模型对空间变异性和空间相关性的反映程度。传统方法中插值模型需要用理论变异函数拟合实验变异函数，而理论变异函数模型往往是根据人的经验来选定，因此会导致不同的理论变异函数模型对降雨量估算结果的优劣程度产生较大的影响。此外，理论变异函数形状固定，有限的已知降雨量站点数据进行拟合时，无法反映实际降雨量在空间上的相关性和差异性，降雨量的空间变化趋势被淹没，而被淹没的这种空间变化趋势往往还具有多尺度的特征，对于尺度的选择和尺度效应的处理也是至关重要的。故如何根据现有降雨量观测站的数据生成空间计算单元内的降雨量，且尽可能地客观反映降雨量空间变化特征，是本专利的研究重点。
发明内容
本发明所要解决的技术问题就是提供一种可选的克里金空间插值降雨量估算方法。
本发明采用如下技术方案：
一种可选的克里金空间插值降雨量估算方法，其改进之处在于，包括如下步骤：
步骤1，通过离散变异函数公式计算所有降雨量样本点对的实验变异函数值：
将降雨量z(x)在点x处和点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure BDA0002728954110000021
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差，
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的降雨量数学期望，
整理得：
Figure BDA0002728954110000022
这时，变异函数γ(x，h)依赖于两个变量：站点位置x和降雨量样本点对的距离h，若变异函数仅仅依赖于距离h而与位置x无关时，则γ(x，h)可以写为γ(h)：
Figure BDA0002728954110000023
离散样本的实验变异函数的计算公式为：
Figure BDA0002728954110000024
其中，i＝1，...，N(h)，h为降雨量样本点对的距离，N(h)代表降雨量样本点对距离为h所有样本点对的个数，z(xi)和z(xi+h)分别是降雨量z(x)在空间位置点xi和xi+h处的降雨量真实值；
步骤2，采用两尺度小波最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型：
步骤21，两尺度小波最小二乘支持向量机通过两尺度的小波核函数来表现空间不同细节的变化，反映空间两尺度特征，两尺度小波最小二乘支持向量机简称两尺度小波LS-SVM；
两尺度小波LS-SVM模型在回归问题上，输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure BDA0002728954110000025
其中，hi∈Rd，在这里d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；
在两尺度问题中，假定尺度1上小波核函数尺度因子偏大，尺度2上小波核函数尺度因子偏小；
在尺度1因子上，f1(x)用集合
Figure BDA0002728954110000031
来回归，尺度1函数模型f1(h)可以用公式(6)表示：
Figure BDA0002728954110000032
其中，h为降雨量样本点对的距离，ω1为尺度1下的权系数向量，ω1 T代表尺度1下权系数的转置向量，b1为常数，
Figure BDA0002728954110000033
为尺度1下输入空间到特征空间的映射函数；
在尺度2因子上，f2(x)用集合
Figure BDA0002728954110000034
来回归，尺度2函数模型f2(h)可以用公式(7)表示：
Figure BDA0002728954110000035
其中，ω2为尺度2下的权系数向量，ω2 T代表尺度2下权系数的转置向量，b2为常数，
Figure BDA0002728954110000036
为尺度2下输入空间到特征空间的映射函数；
将两尺度的模型融合为最终函数模型f(h)，可以表示为：
Figure BDA0002728954110000037
其中，b＝b1+b2，b为常数，
f1(x)是对原始降雨量样本数据回归，f2(x)则是在f1(x)的基础上对残差进行回归，故又可以表示为：
Figure BDA0002728954110000038
Figure BDA0002728954110000039
步骤22，将支持向量机模型转换为优化函数
Figure BDA00027289541100000310
Figure BDA00027289541100000311
其中，i＝1，...，n，n为分组后需要拟合的实验变异函数值总个数，
Figure BDA00027289541100000312
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω1||2为ω1的2-范数平方，||ω2||2为ω2的2-范数平方，
Figure BDA00027289541100000313
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，e1i和e2i分别表示尺度1和尺度2下的误差项，γ1和γ2分别表示尺度1和尺度2下的正则化参数；
在两尺度LS-SVM方法中，误差项e1i等于观测站降雨量真实值yi与回归模型计算的降雨量
Figure BDA0002728954110000041
之差，而误差项e2i等于yi-f1(hi)与回归模型计算值
Figure BDA0002728954110000042
之差，因此优化函数须满足约束条件：
Figure BDA0002728954110000043
Figure BDA0002728954110000044
步骤23，利用拉格朗日乘数法将步骤22中含约束条件公式(12)和(13)的优化函数(11)转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure BDA00027289541100000410
为：
Figure BDA0002728954110000045
其中，α1i和α2i分别为尺度1和尺度2下的拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure BDA0002728954110000046
经拉格朗日方程求解，矩阵形式为：
Figure BDA0002728954110000047
其中，α＝[α11，α12，…，α1n]T，β＝[α21，α22，…，α2n]T，
Figure BDA0002728954110000048
y＝[y1，y2，…，yn]T，I为单位矩阵，Ω1和Ω2为小波核函数；
Figure BDA0002728954110000049
Figure BDA0002728954110000051
得到两尺度LS-SVM模型为：
Figure BDA0002728954110000052
其中，i＝1，...，n，j＝1，...，n，hi和hj分别表示第i个或第j个距离，n为分组后需要拟合的实验变异函数值总个数，α1i表示尺度1下的拉格朗日乘子，α2i和α2j表示尺度2下的拉格朗日乘子，b1和b2为常数，K1(hi，h)和K2(hj，h)分别为尺度1和尺度2下的小波核函数，小波核函数的低频信息作为空间变化的大尺度框架，而高频信息反映空间变化的细节部分，故可以区别不同尺度的空间特征；
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi：
克里金插值方法需要满足两个条件：无偏性和估计方差最小：
E[z(x0)-z*(x0)]＝0 (22)
Var[z(x0)-z*(x0)]＝min (23)
其中，z(x0)为待估位置点x0处的降雨量真实值，z*(x0)为待估位置点x0处的降雨量估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure BDA0002728954110000053
即：
Figure BDA0002728954110000054
估计方差最小表现为：
Figure BDA0002728954110000055
根据拉格朗日乘数原理求条件极值：
Figure BDA0002728954110000056
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即：
Figure BDA0002728954110000061
其中，k为降雨量样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的降雨量z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值，
上述公式用矩阵表示，如下：
Figure BDA0002728954110000062
其中，λ为列向量，代表λi，i＝1，...，k，
普通克里金方程组为：K*λ＝D (30)
解得：λ＝K-1D； (31)
步骤4，根据克里金权重系数λi计算待估位置点的降雨量估计值：
根据克里金权重系数λi计算待估位置点x0处的降雨量估计值z*(x0)，即实现待估位置点降雨量无偏、最优的估计，可用公式(32)进行计算：
Figure BDA0002728954110000063
步骤5，根据所有待估位置点的降雨量估计值和真实值，利用平均绝对误差MAE和均方根误差RMSE进行克里金空间插值降雨量估算方法定量化精度评价：
Figure BDA0002728954110000064
Figure BDA0002728954110000065
其中，S代表待估位置点的总个数，待估位置点的降雨量真实值为yi，i＝1，...，S，待估位置点的降雨量估计值为
Figure BDA0002728954110000066
i＝1，...，S。
进一步的，在步骤2实验变异函数值拟合之前，若需要拟合的实验变异函数值较多，则先进行分组操作；
进一步的，步骤2采用多尺度小波最小二乘支持向量机拟合实验变异函数值，多尺度取2，3，…，N。
本发明的有益效果是：
本发明公开了一种可选的克里金空间插值降雨量估算方法，克服了普通克里金空间插值降雨量估算方法中传统理论变异函数模型形状固定且未考虑空间变化趋势及多尺度特征的不足，以多尺度小波最小二乘支持向量机拟合实验变异函数值，使得降雨量估算的结果符合本身的空间变化趋势且体现空间变化的多尺度特征。
本发明方法从已知降雨量样本数据点的空间变化趋势出发，通过多尺度小波LS-SVM拟合实验变异函数，这种方式拟合的结果符合降雨量本身的空间变化趋势，解决了传统方法无法反映实际空间变化趋势的问题；同时小波核作为LS-SVM核函数，能根据不同部分调整核参数，灵活多样；最后，多尺度小波核利用小波多分辨率的特性，可以体现降雨量空间变化的多尺度特征，体现空间变化的尺度效应。
本发明方法直观的效果体现为高精度的插值估算结果，即空间插值降雨量估算问题中对待估位置点降雨量的准确估计，为政府主管部门及时评估洪涝灾害发生等级和完善防洪布防提供参考依据，深层次的技术效果体现为提高地理数据解释地理现象的能力，这是极具规律性的。
附图说明
图1是本发明实施例1所公开方法的流程示意图；
图2是SIC97数据集的观测点与待估位置点的地理分布图；
图3是所有降雨量样本点对分组前的变异函数云图；
图4是所有降雨量样本点对分组后的变异函数云图；
图5是需要拟合的实验变异函数值；
图6是传统理论变异模型(指数、球状、高斯)拟合曲线图；
图7是单尺度LS-SVM(RBF核、Morlet小波核)和多尺度LS-SVM(Morlet小波核)拟合曲线图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例1，如图1所示，本实施例公开了一种可选的克里金空间插值降雨量估算方法，包括如下步骤：
步骤1，通过离散变异函数公式计算所有降雨量样本点对的实验变异函数值：
将降雨量z(x)在点x处和点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure BDA0002728954110000081
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差，
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的降雨量数学期望，
整理得：
Figure BDA0002728954110000082
这时，变异函数γ(x，h)依赖于两个变量：站点位置x和降雨量样本点对的距离h，若变异函数仅仅依赖于距离h而与位置x无关时，则γ(x，h)可以写为γ(h)：
Figure BDA0002728954110000083
离散变异函数γ(h)是克里金空间插值独有的参量，通过对所有降雨量样本点对的计算，得到该参量值。
离散样本的实验变异函数的计算公式为：
Figure BDA0002728954110000084
其中，i＝1，...，N(h)，h为降雨量样本点对的距离，N(h)代表降雨量样本点对距离为h所有样本点对的个数，z(xi)和z(xi+h)分别是降雨量z(x)在空间位置点xi和xi+h处的降雨量真实值；
步骤2，采用多尺度小波最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型，若需要拟合的实验变异函数值较多，则在拟合之前进行分组操作，以便降低计算过程中的复杂度。
步骤21，多尺度小波最小二乘支持向量机(简称多尺度小波LS-SVM)通过多尺度的小波核函数来表现空间不同细节的变化，反映空间多尺度特征，多尺度可取2，3，…，N。一般情况下，为方便计算取两尺度小波LS-SVM，对于更多尺度小波LS-SVM可以用类似的方法进行推导。
两尺度小波LS-SVM模型在回归问题上，输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure BDA0002728954110000091
其中，hi∈Rd，在这里d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；
在两尺度问题中，假定尺度1上小波核函数尺度因子偏大，尺度2上小波核函数尺度因子偏小；
在尺度1因子上，f1(x)用集合
Figure BDA0002728954110000092
来回归，尺度1函数模型f1(h)可以用公式(6)表示：
Figure BDA0002728954110000093
其中，h为降雨量样本点对的距离，ω1为尺度1下的权系数向量，ω1 T代表尺度1下权系数的转置向量，b1为常数，
Figure BDA0002728954110000094
为尺度1下输入空间到特征空间的映射函数；
在尺度2因子上，f2(x)用集合
Figure BDA0002728954110000095
来回归，尺度2函数模型f2(h)可以用公式(7)表示：
Figure BDA0002728954110000096
其中，ω2为尺度2下的权系数向量，ω2 T代表尺度2下权系数的转置向量，b2为常数，
Figure BDA0002728954110000097
为尺度2下输入空间到特征空间的映射函数；
将两尺度的模型融合为最终函数模型f(h)，可以表示为：
Figure BDA0002728954110000098
其中，b＝b1+b2，b为常数，
需要说明的是，f1(x)是对原始降雨量样本数据回归，f2(x)则是在f1(x)的基础上对残差进行回归，故又可以表示为：
Figure BDA0002728954110000099
Figure BDA00027289541100000910
步骤22，依据统计学习理论，支持向量机模型的目的是使结构风险和经验风险同时达到最小，将支持向量机模型转换为优化函数
Figure BDA00027289541100000911
Figure BDA00027289541100000912
其中，i＝1，...，n，n为分组后需要拟合的实验变异函数值总个数，
Figure BDA00027289541100000913
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω1||2为ω1的2-范数平方，||ω2||2为ω2的2-范数平方，
Figure BDA0002728954110000101
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，e1i和e2i分别表示尺度1和尺度2下的误差项，γ1和γ2分别表示尺度1和尺度2下的正则化参数；
在两尺度LS-SVM方法中，误差项e1i等于观测站降雨量真实值yi与回归模型计算的降雨量
Figure BDA0002728954110000102
之差，而误差项e2i等于yi-f1(hi)与回归模型计算值
Figure BDA0002728954110000103
之差，因此优化函数须满足约束条件：
Figure BDA0002728954110000104
Figure BDA0002728954110000105
步骤23，利用拉格朗日乘数法将步骤22中含约束条件公式(12)和(13)的优化函数(11)转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure BDA0002728954110000109
为：
Figure BDA0002728954110000106
其中，α1i和α2i分别为尺度1和尺度2下的拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure BDA0002728954110000107
经拉格朗日方程求解，矩阵形式为：
Figure BDA0002728954110000108
其中，α＝[α11，α12，…，α1n]T，β＝[α21，α22，…，α2n]T，
Figure BDA0002728954110000111
y＝[y1，y2，…，yn]T，I为单位矩阵，Ω1和Ω2为小波核函数；
Figure BDA0002728954110000112
Figure BDA0002728954110000113
小波函数是一种时频局部化的工具，是把数据、函数或算子分割成不同频率的成分，然后再用分解的方法研究不同尺度下的成分，故非常适合做空间多尺度分析。小波函数的低频信息可作为空间变化的大尺度框架，而高频信息反映空间变化的细节部分，区别不同尺度的空间特征。目前应用较多的实数和复数小波函数中，以Morlet小波、Gaussian小波、Mexican hat小波等为主要应用，各形式如下：
Morlet小波函数形式为：
Figure BDA0002728954110000114
Gaussian小波函数形式为：
Figure BDA0002728954110000115
Mexican hat小波函数形式为：
Figure BDA0002728954110000116
以Morlet小波函数为例，验证其满足Mercer条件和平移不变定理，核函数形式为：
Figure BDA0002728954110000117
其中，d代表维数，在这里d＝1，ω0和ai为待定参数。
最后，得到两尺度LS-SVM模型为：
Figure BDA0002728954110000118
其中，i＝1，....，n，j＝1，...，n，hi和hj分别表示第i个或第j个距离，n为分组后需要拟合的实验变异函数值总个数，α1i表示尺度1下的拉格朗日乘子，α2i和α2j表示尺度2下的拉格朗日乘子，b1和b2为常数，K1(hi，h)和K2(hj，h)分别为尺度1和尺度2下的小波核函数，小波核函数的低频信息作为空间变化的大尺度框架，而高频信息反映空间变化的细节部分，故可以区别不同尺度的空间特征；
采用多尺度最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型，通过多尺度小波支持向量机与克里金空间插值相结合，创新性思维解决了降雨量估算过程中对待估位置点降雨量无偏、最优的估计问题，也是本发明区别于普通克里金插值技术最明显的特征。
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi：
以克里金插值方法的原理进行说明，克里金插值方法需要满足两个条件：无偏性和估计方差最小：
E[z(x0)-z*(x0)]＝0 (22)
Var[z(x0)-z*(x0)]＝min (23)
其中，z(x0)为待估位置点x0处的降雨量真实值，z*(x0)为待估位置点x0处的降雨量估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure BDA0002728954110000121
即：
Figure BDA0002728954110000122
估计方差最小表现为：
Figure BDA0002728954110000123
克里金插值方法在满足无偏性和估计方差最小这两个原则的条件下，建立了含有约束条件的拉格朗日函数。其中，无偏性体现在约束条件上，而方差最小体现在求极值问题上。根据拉格朗日乘数原理求条件极值：
Figure BDA0002728954110000124
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即：
Figure BDA0002728954110000131
其中，k为降雨量样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的降雨量z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值，
上述公式用矩阵表示，如下：
Figure BDA0002728954110000132
其中，λ为列向量，代表λi，i＝1，...，k，通过拟合得到的理论变异函数模型，解方程求解克里金权重系数λi。
普通克里金方程组为：K*λ＝D (30)
解得：λ＝K-1D； (31)
步骤4，根据克里金权重系数λi计算待估位置点的降雨量估计值：
根据克里金权重系数λi计算待估位置点x0处的降雨量估计值z*(x0)，即实现待估位置点降雨量无偏、最优的估计，可用公式(32)进行计算：
Figure BDA0002728954110000133
步骤5，根据所有待估位置点的降雨量估计值和真实值，利用平均绝对误差MAE和均方根误差RMSE进行克里金空间插值降雨量估算方法定量化精度评价，实现对该发明优劣性能评价的参考。评价指标采用平均绝对误差MAE(Mean Absolute Error)和均方根误差RMSE(Root Mean Square Error)。MAE反映估测的总体误差，RMSE反映样本数据的估测灵敏度和极值效应，二者数值越小，表明插值效果越好。其定义如下：
Figure BDA0002728954110000134
Figure BDA0002728954110000135
其中，S代表待估位置点的总个数，待估位置点的降雨量真实值为yi，i＝1，...，S，待估位置点的降雨量估计值为
Figure BDA0002728954110000141
i＝1，...，S。
下面以SIC97数据集为具体实例进行说明：
Spatial Interpolation Comparison 97数据集(简称SIC97数据集)为1986年5月8日瑞士467个降雨量观测点的降雨量数据。其中，样本观测点总个数为100，即k＝100，待估位置点的总个数为367，即S＝367。数据来源：https：//www.researchgate.net/profile/Gregoire_Dubois/publication/281292076_Spatial_Interpolation_Comparison_97_(SIC97)_dataset，该数据集是用来测试空间插值降雨量估算方法的典型数据集。该数据集的观测点与待估位置点地理分布如图2所示，由于样本点数目较多，拟合实验变异函数之前，进行分组操作，分组后需要拟合的实验变异函数值总个数为20，即n＝20，图3和图4分别表示分组前的变异函数云图和分组后的变异函数云图，图5表示需要拟合的实验变异函数值。
多尺度小波LS-SVM中核函数采用Morlet小波核函数，用本实施例所公开两尺度小波LS-SVM优化的克里金空间插值方法进行估计，同时采用了球状模型、指数模型、高斯模型三种传统理论变异函数模型和单尺度RBF核LS-SVM、单尺度Morlet小波核LS-SVM。LS-SVM的惩罚参数和核参数优化都是通过交叉验证的方法获得，传统方法拟合曲线如图6所示，单尺度LS-SVM(RBF核、Morlet小波核)和多尺度小波LS-SVM(两尺度Morlet小波LS-SVM)模型拟合曲线如图7所示，最终得到单尺度RBF核LS-SVM的惩罚参数γ＝168，核参数σ2＝8×103；单尺度Morlet小波核LS-SVM的惩罚参数γ＝681，核参数a＝328；两尺度Morlet小波核LS-SVM的惩罚参数γ1＝3.7×103γ2＝378，尺度1下的核参数ω01＝27a1＝0.03，尺度2下的核参数ω02＝299a2＝14。
通过MAE和RMSE得到的最终插值结果如下表所示：
Figure BDA0002728954110000142
Figure BDA0002728954110000151
从上表中对比发现，在传统方法中，以球状模型精度最优；单尺度RBF核函数LS-SVM模型精度优于高斯模型和指数模型，但略低于球状模型；单尺度Mor l et核函数LS-SVM模型略优于所有的传统方法；两尺度小波LS-SVM模型插值精度是最好的。综合上述评价指标MAE、RMSE结果，可以认为，基于两尺度小波LS-SVM优化的克里金空间插值降雨量估算方法效果最优，为降雨量估算提供了一种可选的估算方法。
为详细具体阐述本发明方法的具体过程，以SIC97数据集为具体实例进行说明，通过克里金空间插值实现降雨量估算是本发明要解决的典型问题，问题本身契合具体的技术发明领域，具备了良好的代表性，因此选择此数据集作为本发明技术方案的具体实例。同时，本发明的技术方案不能脱离技术流程和数学模型而独立存在，各个步骤出现的计算和公式是实现本发明方案的必不可少的技术手段，其技术特征已经体现在数学函数参数所代表的物理意义及具体的应用领域。

Claims (3)

1.一种可选的克里金空间插值降雨量估算方法，其特征在于，包括如下步骤：
步骤1，通过离散变异函数公式计算所有降雨量样本点对的实验变异函数值：
将降雨量z(x)在点x处和点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure FDA0002728954100000011
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差，
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的降雨量数学期望，
整理得：
Figure FDA0002728954100000012
这时，变异函数γ(x，h)依赖于两个变量：站点位置x和降雨量样本点对的距离h，若变异函数仅仅依赖于距离h而与位置x无关时，则γ(x，h)可以写为γ(h)：
Figure FDA0002728954100000013
离散样本的实验变异函数的计算公式为：
Figure FDA0002728954100000014
其中，i＝1，...，N(h)，h为降雨量样本点对的距离，N(h)代表降雨量样本点对距离为h所有样本点对的个数，z(xi)和z(xi+h)分别是降雨量z(x)在空间位置点xi和xi+h处的降雨量真实值；
步骤2，采用两尺度小波最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型：
步骤21，两尺度小波最小二乘支持向量机通过两尺度的小波核函数来表现空间不同细节的变化，反映空间两尺度特征，两尺度小波最小二乘支持向量机简称两尺度小波LS-SVM；
两尺度小波LS-SVM模型在回归问题上，输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure FDA0002728954100000015
其中，hi∈Rd，在这里d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；
在两尺度问题中，假定尺度1上小波核函数尺度因子偏大，尺度2上小波核函数尺度因子偏小；
在尺度1因子上，f1(x)用集合
Figure FDA0002728954100000021
来回归，尺度1函数模型f1(h)可以用公式(6)表示：
Figure FDA0002728954100000022
其中，h为降雨量样本点对的距离，ω1为尺度1下的权系数向量，ω1 T代表尺度1下权系数的转置向量，b1为常数，
Figure FDA0002728954100000023
为尺度1下输入空间到特征空间的映射函数；
在尺度2因子上，f2(x)用集合
Figure FDA0002728954100000024
来回归，尺度2函数模型f2(h)可以用公式(7)表示：
Figure FDA0002728954100000025
其中，ω2为尺度2下的权系数向量，ω2 T代表尺度2下权系数的转置向量，b2为常数，
Figure FDA0002728954100000026
为尺度2下输入空间到特征空间的映射函数；
将两尺度的模型融合为最终函数模型f(h)，可以表示为：
Figure FDA0002728954100000027
其中，b＝b1+b2，b为常数，
f1(x)是对原始降雨量样本数据回归，f2(x)则是在f1(x)的基础上对残差进行回归，故又可以表示为：
Figure FDA0002728954100000028
Figure FDA0002728954100000029
步骤22，将支持向量机模型转换为优化函数
Figure FDA00027289541000000210
Figure FDA00027289541000000211
其中，i＝1，...，n，n为分组后需要拟合的实验变异函数值总个数，
Figure FDA00027289541000000212
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω1||2为ω1的2-范数平方，||ω2||2为ω2的2-范数平方，
Figure FDA00027289541000000213
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，e1i和e2i分别表示尺度1和尺度2下的误差项，γ1和γ2分别表示尺度1和尺度2下的正则化参数；
在两尺度LS-SVM方法中，误差项e1i等于观测站降雨量真实值yi与回归模型计算的降雨量
Figure FDA0002728954100000031
之差，而误差项e2i等于yi-f1(hi)与回归模型计算值
Figure FDA0002728954100000032
之差，因此优化函数须满足约束条件：
Figure FDA0002728954100000033
Figure FDA0002728954100000034
步骤23，利用拉格朗日乘数法将步骤22中含约束条件公式(12)和(13)的优化函数(11)转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure FDA0002728954100000035
为：
Figure FDA0002728954100000036
其中，a1i和α2i分别为尺度1和尺度2下的拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure FDA0002728954100000037
经拉格朗日方程求解，矩阵形式为：
Figure FDA0002728954100000038
其中，α＝[α11，α12，…，α1n]T，β＝[α21，α22，…，α2n]T，
Figure FDA0002728954100000039
y＝[y1，y2，…，yn]T，I为单位矩阵，Ω1和Ω2为小波核函数；
Figure FDA00027289541000000310
Figure FDA0002728954100000041
得到两尺度LS-SVM模型为：
Figure FDA0002728954100000042
其中，i＝1，...，n，j＝1，...，n，hi和hj分别表示第i个或第j个距离，n为分组后需要拟合的实验变异函数值总个数，α1i表示尺度1下的拉格朗日乘子，α2i和α2j表示尺度2下的拉格朗日乘子，b1和b2为常数，K1(hi，h)和K2(hj，h)分别为尺度1和尺度2下的小波核函数，小波核函数的低频信息作为空间变化的大尺度框架，而高频信息反映空间变化的细节部分，故可以区别不同尺度的空间特征；
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi：
克里金插值方法需要满足两个条件：无偏性和估计方差最小：
E[z(x0)-z*(x0)]＝0 (22)
Var[z(x0)-z*(x0)]＝min (23)
其中，z(x0)为待估位置点x0处的降雨量真实值，z*(x0)为待估位置点x0处的降雨量估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure FDA0002728954100000043
即：
Figure FDA0002728954100000044
估计方差最小表现为：
Figure FDA0002728954100000045
根据拉格朗日乘数原理求条件极值：
Figure FDA0002728954100000046
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即：
Figure FDA0002728954100000051
其中，k为降雨量样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的降雨量z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值，
上述公式用矩阵表示，如下：
Figure FDA0002728954100000052
其中，λ为列向量，代表λi，i＝1，...，k，
普通克里金方程组为：K*λ＝D (30)
解得：λ＝K-1D； (31)
步骤4，根据克里金权重系数λi计算待估位置点的降雨量估计值：
根据克里金权重系数λi计算待估位置点x0处的降雨量估计值z*(x0)，即实现待估位置点降雨量无偏、最优的估计，可用公式(32)进行计算：
Figure FDA0002728954100000053
步骤5，根据所有待估位置点的降雨量估计值和真实值，利用平均绝对误差MAE和均方根误差RMSE进行克里金空间插值降雨量估算方法定量化精度评价：
Figure FDA0002728954100000054
Figure FDA0002728954100000055
其中，S代表待估位置点的总个数，待估位置点的降雨量真实值为yi，i＝1，...，S，待估位置点的降雨量估计值为
Figure FDA0002728954100000056
i＝1，...，S。
2.根据权利要求1所述可选的克里金空间插值降雨量估算方法，其特征在于：在步骤2实验变异函数值拟合之前，若需要拟合的实验变异函数值较多，则先进行分组操作。
3.根据权利要求1所述可选的克里金空间插值降雨量估算方法，其特征在于：步骤2采用多尺度小波最小二乘支持向量机拟合实验变异函数值，多尺度取2，3，…，N。
CN202011112171.5A 2020-10-16 2020-10-16 一种可选的克里金空间插值降雨量估算方法 Pending CN112257021A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011112171.5A CN112257021A (zh) 2020-10-16 2020-10-16 一种可选的克里金空间插值降雨量估算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011112171.5A CN112257021A (zh) 2020-10-16 2020-10-16 一种可选的克里金空间插值降雨量估算方法

Publications (1)

Publication Number Publication Date
CN112257021A true CN112257021A (zh) 2021-01-22

Family

ID=74245670

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011112171.5A Pending CN112257021A (zh) 2020-10-16 2020-10-16 一种可选的克里金空间插值降雨量估算方法

Country Status (1)

Country Link
CN (1) CN112257021A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112612995A (zh) * 2021-03-08 2021-04-06 武汉理工大学 一种基于贝叶斯回归的多来源降雨数据融合算法及装置
CN112632472A (zh) * 2021-03-09 2021-04-09 武汉理工大学 一种基于空间统计学的多来源降雨数据融合算法及装置
CN116580123A (zh) * 2023-03-29 2023-08-11 云南省生态环境科学研究院 一种区域大气环境光化学网格模型图像生成系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6961719B1 (en) * 2002-01-07 2005-11-01 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Hybrid neural network and support vector machine method for optimization
CN105043390A (zh) * 2015-06-29 2015-11-11 中国船舶重工集团公司第七0七研究所 基于泛克里金法的重力场插值方法
CN106600534A (zh) * 2016-12-12 2017-04-26 中国石油大学(华东) 一种基于多尺度小波支持向量机优化的克里金空间插值方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6961719B1 (en) * 2002-01-07 2005-11-01 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Hybrid neural network and support vector machine method for optimization
CN105043390A (zh) * 2015-06-29 2015-11-11 中国船舶重工集团公司第七0七研究所 基于泛克里金法的重力场插值方法
CN106600534A (zh) * 2016-12-12 2017-04-26 中国石油大学(华东) 一种基于多尺度小波支持向量机优化的克里金空间插值方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112612995A (zh) * 2021-03-08 2021-04-06 武汉理工大学 一种基于贝叶斯回归的多来源降雨数据融合算法及装置
CN112612995B (zh) * 2021-03-08 2021-07-09 武汉理工大学 一种基于贝叶斯回归的多来源降雨数据融合算法及装置
CN112632472A (zh) * 2021-03-09 2021-04-09 武汉理工大学 一种基于空间统计学的多来源降雨数据融合算法及装置
CN116580123A (zh) * 2023-03-29 2023-08-11 云南省生态环境科学研究院 一种区域大气环境光化学网格模型图像生成系统

Similar Documents

Publication Publication Date Title
CN112257021A (zh) 一种可选的克里金空间插值降雨量估算方法
Guénot et al. Adaptive sampling strategies for non‐intrusive POD‐based surrogates
US6208982B1 (en) Method and apparatus for solving complex and computationally intensive inverse problems in real-time
Ma et al. Phased microphone array for sound source localization with deep learning
CN106600534A (zh) 一种基于多尺度小波支持向量机优化的克里金空间插值方法
Hodyss et al. The error of representation: Basic understanding
Constantinescu et al. Physics-based covariance models for Gaussian processes with multiple outputs
Park et al. Unstructured grid adaptation and solver technology for turbulent flows
Sasaki et al. Transfer functions for flow predictions in wall-bounded turbulence
CN115203865B (zh) 一种基于数字孪生的产品装配过程机械性能在线预测方法
CN110837921A (zh) 基于梯度提升决策树混合模型的房地产价格预测研究方法
Solonen et al. On dimension reduction in Gaussian filters
Cheng et al. Background error covariance iterative updating with invariant observation measures for data assimilation
Cheng et al. A graph clustering approach to localization for adaptive covariance tuning in data assimilation based on state-observation mapping
CN115983029A (zh) 一种航空装备可靠性仿真数字孪生模型构建方法、设备、介质
Drakoulas et al. FastSVD-ML–ROM: A reduced-order modeling framework based on machine learning for real-time applications
CN109540089B (zh) 一种基于贝叶斯-克里金模型的桥面高程拟合方法
Gálvez et al. Immunological-based approach for accurate fitting of 3D noisy data points with Bézier surfaces
CN113435089B (zh) 一种基于高斯过程的板材折弯回弹预测方法
Regazzoni et al. A physics-informed multi-fidelity approach for the estimation of differential equations parameters in low-data or large-noise regimes
Mooers et al. Generative modeling of atmospheric convection
Mehra et al. An adaptive spectral graph wavelet method for PDEs on networks
CN109960146A (zh) 提高软测量仪表模型预测精度的方法
CN112949944A (zh) 一种基于时空特征的地下水位智能预测方法及系统
Cortés et al. Geometry simplification of open-cell porous materials for elastic deformation FEA

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20210122

RJ01 Rejection of invention patent application after publication