CN116541681A - 一种基于协同克里金插值的复合灾害空间变异性识别方法 - Google Patents
一种基于协同克里金插值的复合灾害空间变异性识别方法 Download PDFInfo
- Publication number
- CN116541681A CN116541681A CN202310437686.XA CN202310437686A CN116541681A CN 116541681 A CN116541681 A CN 116541681A CN 202310437686 A CN202310437686 A CN 202310437686A CN 116541681 A CN116541681 A CN 116541681A
- Authority
- CN
- China
- Prior art keywords
- variable
- sequence
- variation function
- variable sequence
- value
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 239000002131 composite material Substances 0.000 title claims abstract description 51
- 238000009826 distribution Methods 0.000 claims abstract description 63
- 238000011160 research Methods 0.000 claims abstract description 24
- 238000004519 manufacturing process Methods 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 claims description 180
- 238000002790 cross-validation Methods 0.000 claims description 19
- 238000004364 calculation method Methods 0.000 claims description 15
- 238000010219 correlation analysis Methods 0.000 claims description 13
- 238000004590 computer program Methods 0.000 claims description 9
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 claims description 9
- 239000010931 gold Substances 0.000 claims description 9
- 229910052737 gold Inorganic materials 0.000 claims description 9
- 238000012360 testing method Methods 0.000 claims description 9
- 230000000694 effects Effects 0.000 claims description 7
- 238000012795 verification Methods 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 238000012549 training Methods 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 14
- 238000004458 analytical method Methods 0.000 description 8
- 238000006243 chemical reaction Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 150000001875 compounds Chemical class 0.000 description 2
- 238000011835 investigation Methods 0.000 description 2
- 238000010998 test method Methods 0.000 description 2
- 238000006424 Flood reaction Methods 0.000 description 1
- 230000007488 abnormal function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000003653 coastal water Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000001939 inductive effect Effects 0.000 description 1
- 238000012502 risk assessment Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000010792 warming Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/24—Querying
- G06F16/245—Query processing
- G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
- G06F16/2462—Approximate or statistical queries
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/29—Geographical information databases
-
- 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Databases & Information Systems (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- Probability & Statistics with Applications (AREA)
- Evolutionary Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Algebra (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computational Linguistics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Operations Research (AREA)
- Remote Sensing (AREA)
- Fuzzy Systems (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于协同克里金插值的复合灾害空间变异性识别方法，包括：S1、采集研究区域内多种相关灾害数据，制作相关灾害数据集；S2、确定研究区域的复合灾害事件的主变量和次变量；S3、采用协同克里金插值法，对符合正态分布的主变量序列和次变量序列进行分析评估；S4、根据协同半变异函数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。本发明通过引入协同克里金方法，结合复合灾害是由多种因素共同引发的这一特征，利用多变量对其主要因素进行分析，其空间插值结果经过验证，更加精确，所得到的半变异函数参数可表征空间变异的程度。
Description
技术领域
本发明涉及水灾害及其衍生地质灾害分析技术，特别是涉及一种基于协同克里金插值的复合灾害空间变异性识别方法。
背景技术
由于全球变暖带来的世界范围气候变化加剧，沿海地区自然气候条件发生了巨大改变，极端灾害性事件频发。多个致灾因相叠加形成的极端事件被称为复合灾害。复合灾害事件具有影响范围广、灾害损失大等特点，如何准确、客观且高效地评估复合灾害事件的灾害风险，对有效防范沿海地区的极端灾害有重要指导意义。
但是，在高强度人类活动的多年影响下，沿海地区水文条件及河流地貌特征已发生显著变异，且因沿海地区河网密度及水文站布设位置所限，现有的水文统计资料远不能满足对复合灾害的风险评估的需求。
发明内容
发明目的：本发明的目的是提供一种基于协同克里金插值法的复合灾害空间变异性识别方法。
技术方案：本发明的一种基于协同克里金插值的复合灾害空间变异性识别方法，包括以下步骤：
S1、采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
S2、根据相关灾害数据集确定复合灾害事件的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，使用Spearman和/或Kendall相关系数对两两组合进行相关性分析，选择相关性高且显著的数据序列作为次变量序列；对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
S3、采用协同克里金插值法，对步骤S2得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，采用均方根误差RMSE对各模型拟合结果进行评价，选择RMSE值最小的模型为最优模型，得到最佳主变量半变异函数、最佳次变量半变异函数与最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
S4、根据步骤S3得到的最佳协同半变异函数的参数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。
进一步的，步骤S2中Spearman相关系数γs计算公式为：
其中，n为进行相关性分析的相关灾害数据集的长度；Ri和Si分别为主变量序列的秩和与主变量序列两两组合的其他数据序列的秩，和/>分别为主变量序列的平均值和与主变量序列两两组合的其他数据序列的平均值；
Kendall相关系数τ计算公式为：
其中，xi和xj分别为第i个和第j个主变量序列的数据，yi和yj分别为第i个和第j个与主变量序列两两组合的其他数据序列的数据。
进一步的，步骤S3中对于主变量序列或次变量序列的单变量，其半变异函数公式为：
其中，γ(h)为计算出的主变量或次变量半变异函数值序列，χi为主变量序列或次变量序列中的数据点，χi+h为主变量序列或次变量序列中与数据点χi相距直线距离为h的数据点，h为数据点χi与数据点χi+h之间的直线距离；N(h)为主变量序列或次变量序列中相距直线距离为h的数据点的对数；Z(χi)和Z(χi+h)分别为主变量序列或次变量序列中数据点χi与数据点χi+h处的变量值；
对于主变量序列与次变量序列，两者组成的协同半变异函数公式为：
其中，γ12(h)为计算出的协同半变异函数值序列，Z1(xi)和Z1(xi+h)分别为主变量序列中数据点xi与数据点xi+h处的变量值；Z2(y′j)和Z2(y′j+h)分别为次变量序列中数据点y′j与数据点y′j+h处的变量值。
进一步的，步骤S3中拟合模型包括：
(1)球状模型：
其中，γ1(h)为球状模型拟合出的半变异函数或协同半变异函数；C0为块金值；C为拱高，即偏基台值；(C0+C)为基台值；a为变程；h为主变量序列或次变量序列中数据点χi与数据点χi+h之间的直线距离；
(2)高斯模型：
其中，γ2(h)为高斯模型拟合出的半变异函数或协同半变异函数；
(3)指数模型：
其中，γ3(h)为指数模型拟合出的半变异函数或协同半变异函数；
(4)幂指数模型：
γ4(h)＝Ahθ0<θ<2
其中，γ4(h)为幂指数模型拟合出的半变异函数或协同半变异函数；A为常数；θ为幂指数。
进一步的，步骤S3中根据最佳半变异函数和最佳协同半变异函数，对主变量序列和次变量序列进行空间插值的具体过程包括：
计算协同克里金法的点位插值，计算公式为：
其中，Z*(a0)为待估点a0处主变量预测值；n1和n2分别为主变量序列和次变量序列的样本长度，Z1(xi)为主变量序列中数据点xi处的变量值，Z2(y’j)为次变量序列中数据点y’j处的变量值；λ1i为计算插值时，赋予Z1(xi)的权重系数，λ2j则为赋予Z2(y’j)的权重系数，且
利用一阶平稳条件与权重系数的假设，结合拉格朗日函数，得到如下矩阵：
其中，C1为主变量的协方差，C2为次变量的协方差，C12为两个变量间的互协方差，μ1和μ2为拉格朗日系数；
通过求解权重系数λ1i与λ2j，便得到待估点a0处的主变量预测值Z*(a0)，求出研究区域所有待估点处的预测值，便得到主变量空间插值结果。
进一步的，步骤S3中对于主变量和次变量的空间插值结果，分别运用K折交叉验证的方法分别对其结果进行评价；K折交叉验证具体方法为：将主变量序列或次变量序列平均分为K份，其中K-1份作为训练样本，剩余1份作为验证样本，重复K次，即得到交叉验证后的预测变量值；随后采用均方根误差RMSE来分别对主变量和次变量的残差进行评估，RMSE值越小越好；
RMSE的计算公式为：
其中，Zpred(χi)为交叉验证后的主变量序列或次变量序列中数据点χi处的预测变量值，Z(χi)为主变量序列或次变量序列中数据点χi处的变量值，N为主变量序列或次变量序列的长度。
进一步的，步骤S4具体为：
协同半变异函数拟合得到的参数包括块金值、偏基台值和变程，块金值C0的大小即块金效应，表示随机部分的空间异质性；基台值C+C0为当样点距离h增大到一定值后，半变异函数值将达到的稳定常数值，表示变量的最大变异程度；而半变异函数值达到基台值时的样点间距离称为变程a，因此在变程范围内，具有空间相关性，变程范围外，不具有空间相关性；区域内空间相关度用C/(C+C0)表示，比值越小，空间相关性越强。
本发明的一种基于协同克里金插值的复合灾害空间变异性识别系统，包括：
数据采集及处理模块，用于采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
相关性分析模块，用于根据相关灾害数据集确定复合灾害事件的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，并对两两组合进行相关性分析，选择相关性高且显著的数据序列作为次变量序列；
正态性检验模块，用于对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
协同克里金插值模块，用于采用协同克里金插值法，对得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，采用均方根误差RMSE对各模型拟合结果进行评价，选择RMSE值最小的模型为最优模型，得到最佳主变量半变异函数、最佳次变量半变异函数与最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
识别模块，用于根据得到的最佳协同半变异函数的参数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。
本发明的一种装置设备，包括存储器和处理器，其中：
存储器，用于存储能够在处理器上运行的计算机程序；
处理器，用于在运行所述计算机程序时，执行如上述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤。
本发明的一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如上所述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤。
有益效果：与现有技术相比，本发明技术方案的显著技术效果为：通过引入协同克里金方法，结合复合灾害是由多种因素共同引发的这一特征，利用多变量对其主要因素进行分析，其空间插值结果经过验证，更加精确，所得到的半变异函数参数可表征空间变异的程度，可对沿海地区水文统计资料进行插值模拟，并对复合灾害特征进行空间变异性识别与分析。
附图说明
图1是本发明方法流程图；
图2是研究区水文站点分布图；
图3是降雨量、水位值半变异函数以及协同半变异函数散点图；
图4是协同克里金插值法计算的降雨量空间分布图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
由多种致灾因子相叠加形成的灾害称为复合灾害，如在沿海地区强降雨碰上天文大潮以及台风风暴潮引发的复合洪水，同时极端灾害事件也会带来一系列衍生灾害。复合灾害发生时间内的各水文特性之间具有相关性，因此可以采用协同克里金对复合灾害进行空间插值并识别其空间变异性。本发明的一种基于协同克里金插值的复合灾害空间变异识别方法，首先，收集整理研究区域内相关灾害数据(如降雨事件强度、沿海地区潮位、相关衍生灾害及其损失等)，选定主要灾害事件后，计算主变量数据与其他事件数据间的相关性，选取相关性较强的两类数据表征复合灾害事件，将它们分别作为主/次变量。随后，选定主/次变量后进行正态性检验，采用协同克里金法，拟合出主/次变量的半变异函数与协同半变异函数，通过交叉验证，对空间插值结果进行评价分析。最后，根据拟合出的最佳协同半变异函数相关参数以及主变量的空间分布特征，可对研究区域内的复合灾害特征进行空间变异性识别与分析。
如图1所示，本发明的一种基于协同克里金插值的复合灾害空间变异性识别方法，包括以下步骤：
S1、采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
本发明实施例中，采集了来自珠江三角洲河网区域内的多种水文数据并进行筛选剔除，以下仅列出剔除了多种相关灾害数据中整体分布距离较远的数据后，制作成的相关灾害数据集。相关灾害数据集选用来自珠江三角洲河网区域内共78个水文站点的数据，78个水文站点中降雨测站35个，水位测站43个，相关灾害数据包括35个降雨测站的降雨数据与43个水位测站的水位数据。相关灾害数据集包括水文站的经纬度以及对应位置上的降雨或水位数据。各水文站点位置如图2所示。
S2、选定主变量和次变量；
根据相关灾害数据集确定需要研究的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，使用Spearman和/或Kendall相关系数对两两组合进行相关性分析，选择相关性高且显著的一组作为研究样本，即研究样本包括主变量序列样本和次变量序列样本；对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
本发明实施例中，根据收集到的数据进行相关性分析，以下仅列出相关性高且显著的研究样本作为本例分析所用的主变量序列样本和次变量序列样本。本例选择了降雨与水位数据作为引发复合洪水的重要因素，以降雨量作为主变量，降雨与海岸水位在一定程度上由共同的大尺度天气控制，因此可以选择水位作为次变量辅助计算。以下为使用Spearman与Kendall相关系数对其进行相关性分析的结果。
两类相关系数的计算公式如下：
Spearman相关系数：
其中，n为进行相关性分析的相关灾害数据集的长度；Ri和Si分别为主变量序列的秩和与主变量序列两两组合的其他数据序列的秩，和/>分别为主变量序列的平均值和与主变量序列两两组合的其他数据序列的平均值；
Kendall相关系数：
其中，xi和xj分别为第i和第j个主变量序列的数据，yi和yj分别为第i和第j个与主变量序列两两组合的其他数据序列的数据。
通过两变量不存在相关性的概率值(p值)来确定相关性是否显著，在本发明中，当p<0.05时，即表示数据通过了95％的显著性检验。
计算得到二者的相关系数如表1所示。
表1降雨量与水位的相关系数及p值
由表1可知，降雨量与水位之间存在着正相关性，且p值均小于0.05，具有95％的显著相关，因此可以用水位作为次变量来辅助得到降雨量的空间分布。
对主变量序列与次变量序列分别进行正态性检验，根据所选主变量序列/次变量序列的长度选择不同的检验方法，如Shapiro-Wilk(SW)检验与Kolmogorov-Smirnov(KS)检验，若数据序列为非正态分布，则需要将其转换为正态分布，常见的转换方式有对数转换、倒数转换和平方根转换等，利用转换后的主/次变量序列进行后续计算。
本发明实施例中，采用Q-Q图和夏皮罗-威尔克(Shapiro-Wilk)检验法，分别对降雨数据序列与水位数据序列进行检验，发现降雨数据序列需要通过平方根转换为正态分布，而水位数据序列符合正态分布，无需转换。
S3、采用协同克里金Cokriging插值法，对步骤S2得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，采用多种模型分别对主变量序列、次变量序列以及两者的组合进行拟合，得到拟合后的多种主变量半变异函数、次变量半变异函数以及两者的协同半变异函数，采用均方根误差RMSE对各拟合模型结果进行评价，选择RMSE值最小的模型为最优模型，最优模型拟合出的半变异函数和协同半变异函数分别为最佳半变异函数和最佳协同半变异函数，利用最佳协同半变异函数参数(块金值、偏基台值、变程及空间相关度)进行变异情况分析；其次，根据最佳半变异函数和最佳协同半变异函数，采用协同克里金插值方法对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
通过协同克里金法，可对研究区域内的已经过正态分布转换的主/次变量进行分析评估；本步骤分为三部分，首先根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，得到拟合后的多种主变量半变异函数、次变量半变异函数以及两者的协同半变异函数，采用均方根误差RMSE对各模型拟合结果进行评价，选择RMSE值最小的模型为最优模型，最优模型拟合出的半变异函数和协同半变异函数分别为最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证。
S31、最佳半变异函数拟合
对转换为正态分布后的主变量/次变量的半变异函数以及协同半变异函数计算，通过R软件计算并分别绘制出函数的散点图，如图3所示。
对于单变量，其半变异函数公式为：
其中，γ(h)为计算出的主变量或次变量半变异函数值序列，χi为主变量序列或次变量序列中的数据点，χi+h为主变量序列或次变量序列中与数据点χi相距直线距离为h的数据点，h为数据点χi与数据点χi+h之间的直线距离；N(h)为主变量序列或次变量序列中相距直线距离为h的数据点的对数；Z(χi)和Z(χi+h)分别为主变量序列或次变量序列中数据点χi与数据点χi+h处的变量值；
对于主变量序列与次变量序列，两者组成的协同半变异函数公式为：
其中，γ12(h)为计算出的协同半变异函数值序列，Z1(xi)和Z1(xi+h)分别为主变量序列中数据点xi与数据点xi+h处的变量值；Z2(y′j)和Z2(y′j+h)分别为次变量序列中数据点y′j与数据点y′j+h处的变量值。
在对半变异函数的拟合中，根据散点图的分布及趋势选用合适的连续分布形函数来表征散点图的特征，本例选取了4个模型来对数据进行拟合，分别为：
(1)球状模型
其中，γ1(h)为球状模型拟合出的半变异函数或协同半变异函数；C0为块金值；C为拱高，即偏基台值；(C0+C)为基台值；a为变程；h为主变量序列或次变量序列中数据点χi与数据点χi+h之间的直线距离；
(2)高斯模型
其中，γ2(h)为高斯模型拟合出的半变异函数或协同半变异函数；
(3)指数模型
其中，γ3(h)为指数模型拟合出的半变异函数或协同半变异函数；
(4)幂指数模型
γ4(h)＝Ahθ0<θ<2
其中，γ4(h)为幂指数模型拟合出的半变异函数或协同半变异函数；A为常数，θ为幂指数。
为了从半变异函数模型中选出最佳拟合模型，采用均方根误差(RMSE)来对拟合模型结果进行评价，选择RMSE值最小的模型为最优模型，利用其模型参数进行变异情况分析。
RMSE的计算公式为：
其中，Zpred(χi)为交叉验证后的主变量序列或次变量序列中数据点χi处的预测变量值，Z(χi)为主变量序列或次变量序列中数据点χi处的变量值，N为主变量序列或次变量序列的长度。
4种拟合模型的均方根误差如表2所示。
表2不同半变异函数模型均方根误差(RMSE)拟合精度
由表2可知，应选择高斯模型来对主/次变量半变异函数以及协同半变异函数进行拟合。
S32、空间插值计算
通过协同克里金法，对研究区域内的变量的空间分布进行分析评估。
协同克里金法的点位插值计算公式为：
其中，Z*(a0)为待估点a0处主变量预测值；n1和n2分别为主变量序列和次变量序列的样本长度，Z1(xi)为主变量序列中数据点xi处的变量值，Z2(y’j)为次变量序列中数据点y’j处的变量值；；λ1i为计算插值时，赋予Z1(xi)的权重系数，λ2j则为赋予Z2(y’j)的权重系数，且
利用一阶平稳条件与权重系数的假设，结合拉格朗日函数，得到如下矩阵：
其中，C1为主变量的协方差，C2为次变量的协方差，C12为两个变量间的互协方差，μ1和μ2为拉格朗日系数；
通过求解权重系数λ1i与λ2j，便得到待估点a0处的主变量预测值Z*(a0)。首先，在研究区域内，建立400×400的网格，随后通过R软件，用最佳半变异函数和最佳协同半变异函数来进行协同克里金插值计算，得到插值后的降雨量空间分布(图4)。
由图4可知，降雨量的空间分布整体呈现出自东北向西南逐渐升高的趋势，但是在西南方向的台山市东北部，降雨量呈现出由四周向中心逐渐降低的趋势。
S33、交叉验证
对于降雨量的空间插值结果，运用K折交叉验证的方法对结果进行评价。K折交叉验证具体方法为：将样本平均分为K份，其中K-1份作为训练样本，剩余1份作为验证样本，重复K次，即得到交叉验证后的预测降雨量。随后采用均方根误差(RMSE)来分别对降雨量和水位值的残差进行评估，RMSE值越小越好。
RMSE的计算公式为：
其中，Zpred(χi)为交叉验证的预测结果。本研究中采用5折交叉验证法，主变量与次变量交叉验证的结果如表3所示。
表3交叉验证结果
由表3可知，相对于次变量水位，主变量降雨的预测值与观测值之间的均方根误差较小，说明协同克里金法能够对主变量进行更加准确的空间插值。
S4、复合灾害空间变异性识别；
根据步骤S3得到的协同半变异函数的参数以及主变量插值结果，可对复合灾害进行空间变异性识别；
协同半变异函数拟合得到的参数在空间上具有不同的意义。块金值C0的大小即块金效应，表示随机部分的空间异质性；基台值C+C0为，当样点距离h增大到一定值后，半变异函数值将达到的稳定常数值，表示变量的最大变异程度；而半变异函数值达到基台值时的样点间距离称为变程a，因此在变程范围内，具有空间相关性，变程范围外，不具有空间相关性；区域内空间相关度用C/(C+C0)表示，比值越小，空间相关性越强。
主变量空间插值结果可在研究区域内呈现，通过其空间分布特征，结合协同半变异函数参数，进一步对灾害空间变异性进行识别。
本发明实施例中，最优模型为高斯模型，高斯模型拟合的协同半变异函数参数为：
偏基台值C＝1.82，块金值C0＝1.93，变程a＝65.84km。因此可对研究区域进行空间变异性分析。区域内的随机部分具有空间异质性，在最小取样尺度下，误差为1.93；区域内最大变异程度为3.75km2，当两点之间的直线距离超过变程65.84km时，降雨量分布的空间相关性将逐渐降低；空间相关度C/(C+C0)＝0.49，表明所研究的两变量具有中等空间相关性。由降雨量空间插值分布可知，在台山市东北部，发生空间变异，相比于附近区域，台山市东北部的降雨量有降低趋势。
本发明的一种基于协同克里金插值的复合灾害空间变异性识别系统，包括：
数据采集及处理模块，用于采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
相关性分析模块，用于根据相关灾害数据集确定复合灾害事件的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，并对两两组合进行相关性分析，选择相关性高且显著的数据序列作为次变量序列；
正态性检验模块，用于对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
协同克里金插值模块，用于采用协同克里金插值法，对得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，得到最佳主变量半变异函数、最佳次变量半变异函数与最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
识别模块，用于根据得到的最佳协同半变异函数的参数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。
本发明的一种装置设备，包括存储器和处理器，其中：
存储器，用于存储能够在处理器上运行的计算机程序；
处理器，用于在运行所述计算机程序时，执行如上所述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤，并能达到上述方法所述的技术效果。
本发明的一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如上所述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤，并能达到上述方法所述的技术效果。
Claims (10)
1.一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，包括以下步骤：
S1、采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
S2、根据相关灾害数据集确定复合灾害事件的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，使用Spearman和/或Kendall相关系数对两两组合进行相关性分析，选择相关性高且显著的数据序列作为次变量序列；对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
S3、采用协同克里金插值法，对步骤S2得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，采用均方根误差RMSE对各模型拟合结果进行评价，选择RMSE值最小的模型为最优模型，得到最佳主变量半变异函数、最佳次变量半变异函数与最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
S4、根据步骤S3得到的最佳协同半变异函数的参数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。
2.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S2中Spearman相关系数γs计算公式为：
其中，n为进行相关性分析的相关灾害数据集的长度；Ri和Si分别为主变量序列的秩和与主变量序列两两组合的其他数据序列的秩，R和S分别为主变量序列的平均值和与主变量序列两两组合的其他数据序列的平均值；
Kendall相关系数τ计算公式为：
其中，xi和xj分别为第i个和第j个主变量序列的数据，yi和yj分别为第i个和第j个与主变量序列两两组合的其他数据序列的数据。
3.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S3中对于主变量序列或次变量序列的单变量，其半变异函数公式为：
其中，γ(h)为计算出的主变量或次变量半变异函数值序列，χi为主变量序列或次变量序列中的数据点，χi+h为主变量序列或次变量序列中与数据点χi相距直线距离为h的数据点，h为数据点χi与数据点χi+h之间的直线距离；N(h)为主变量序列或次变量序列中相距直线距离为h的数据点的对数；Z(χi)和Z(χi+h)分别为主变量序列或次变量序列中数据点χi与数据点χi+h处的变量值；
对于主变量序列与次变量序列，两者组成的协同半变异函数公式为：
其中，γ12(h)为计算出的协同半变异函数值序列，Z1(xi)和Z1(xi+h)分别为主变量序列中数据点xi与数据点xi+h处的变量值；Z2(y′j)和Z2(y′j+h)分别为次变量序列中数据点y′j与数据点y′j+h处的变量值。
4.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S3中拟合模型包括：
(1)球状模型：
其中，γ1(h)为球状模型拟合出的半变异函数或协同半变异函数；C0为块金值；C为拱高，即偏基台值；(C0+C)为基台值；a为变程；h为主变量序列或次变量序列中数据点χi与数据点χi+h之间的直线距离；
(2)高斯模型：
其中，γ2(h)为高斯模型拟合出的半变异函数或协同半变异函数；
(3)指数模型：
其中，γ3(h)为指数模型拟合出的半变异函数或协同半变异函数；
(4)幂指数模型：
γ4(h)＝Ahθ 0<θ<2
其中，γ4(h)为幂指数模型拟合出的半变异函数或协同半变异函数；A为常数；θ为幂指数。
5.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S3中根据最佳半变异函数和最佳协同半变异函数，对主变量序列和次变量序列进行空间插值的具体过程包括：
计算协同克里金法的点位插值，计算公式为：
其中，Z*(a0)为待估点a0处主变量预测值；n1和n2分别为主变量序列和次变量序列的样本长度，Z1(xi)为主变量序列中数据点xi处的变量值，Z2(y’j)为次变量序列中数据点y’j处的变量值；λ1i为计算插值时，赋予Z1(xi)的权重系数，λ2j则为赋予Z2(y’j)的权重系数，且
利用一阶平稳条件与权重系数的假设，结合拉格朗日函数，得到如下矩阵：
其中，C1为主变量的协方差，C2为次变量的协方差，C12为两个变量间的互协方差，μ1和μ2为拉格朗日系数；
通过求解权重系数λ1i与λ2j，便得到待估点a0处的主变量预测值Z*(a0)，求出研究区域所有待估点处的预测值，便得到主变量空间插值结果。
6.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S3中对于主变量和次变量的空间插值结果，分别运用K折交叉验证的方法分别对其结果进行评价；K折交叉验证具体方法为：将主变量序列或次变量序列平均分为K份，其中K-1份作为训练样本，剩余1份作为验证样本，重复K次，即得到交叉验证后的预测变量值；随后采用均方根误差RMSE来分别对主变量和次变量的残差进行评估，RMSE值越小越好；
RMSE的计算公式为：
其中，Zpred(χi)为交叉验证后的主变量序列或次变量序列中数据点χi处的预测变量值，Z(χi)为主变量序列或次变量序列中数据点χi处的变量值，N为主变量序列或次变量序列的长度。
7.根据权利要求1所述的一种基于协同克里金插值的复合灾害空间变异性识别方法，其特征在于，步骤S4具体为：
协同半变异函数拟合得到的参数包括块金值、偏基台值和变程，块金值C0的大小即块金效应，表示随机部分的空间异质性；基台值C+C0为当样点距离h增大到一定值后，半变异函数值将达到的稳定常数值，表示变量的最大变异程度；而半变异函数值达到基台值时的样点间距离称为变程a，因此在变程范围内，具有空间相关性，变程范围外，不具有空间相关性；区域内空间相关度用C/(C+C0)表示，比值越小，空间相关性越强。
8.一种基于协同克里金插值的复合灾害空间变异性识别系统，其特征在于，包括：
数据采集及处理模块，用于采集研究区域内多种相关灾害数据，剔除多种相关灾害数据中整体分布距离较远的数据，制作相关灾害数据集，相关灾害数据集包括剔除多种相关灾害数据中整体分布距离较远的数据后的多种相关灾害数据及其经纬度；
相关性分析模块，用于根据相关灾害数据集确定复合灾害事件的主变量序列，并将主变量序列与相关灾害数据集中的其他数据序列两两组合，并对两两组合进行相关性分析，选择相关性高且显著的数据序列作为次变量序列；
正态性检验模块，用于对主变量序列与次变量序列分别进行正态性检验，若主变量序列和次变量序列为非正态分布，则将主变量序列和次变量序列分别转换为正态分布；
协同克里金插值模块，用于采用协同克里金插值法，对得到的符合正态分布的主变量序列和次变量序列进行分析评估；首先，根据符合正态分布的主变量序列或次变量序列，分别通过半变异函数公式计算出对应的主变量半变异函数值序列和次变量半变异函数值序列，通过协同半变异函数公式计算出对应的主变量和次变量的协同半变异函数值序列；其次，对得到的主变量半变异函数值序列、次变量半变异函数值序列、协同半变异函数值序列分别进行多模型拟合，采用均方根误差RMSE对各模型拟合结果进行评价，选择RMSE值最小的模型为最优模型，得到最佳主变量半变异函数、最佳次变量半变异函数与最佳协同半变异函数；随后根据最佳主变量半变异函数、最佳次变量半变异函数和最佳协同半变异函数，采用协同克里金插值法来对主变量序列和次变量序列进行空间插值，得到主变量在研究区域内的空间分布；然后，采用均方根误差对主变量序列和次变量序列的空间插值进行交叉验证；
识别模块，用于根据得到的最佳协同半变异函数的参数以及主变量序列的空间插值，对复合灾害进行空间变异性识别。
9.一种装置设备，其特征在于，包括存储器和处理器，其中：
存储器，用于存储能够在处理器上运行的计算机程序；
处理器，用于在运行所述计算机程序时，执行如权利要求1-7任一项所述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤。
10.一种存储介质，其特征在于，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如权利要求1-7任一项所述一种基于协同克里金插值的复合灾害空间变异性识别方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310437686.XA CN116541681A (zh) | 2023-04-21 | 2023-04-21 | 一种基于协同克里金插值的复合灾害空间变异性识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310437686.XA CN116541681A (zh) | 2023-04-21 | 2023-04-21 | 一种基于协同克里金插值的复合灾害空间变异性识别方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116541681A true CN116541681A (zh) | 2023-08-04 |
Family
ID=87456956
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310437686.XA Pending CN116541681A (zh) | 2023-04-21 | 2023-04-21 | 一种基于协同克里金插值的复合灾害空间变异性识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116541681A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117293826A (zh) * | 2023-11-27 | 2023-12-26 | 山东大学 | 一种分布式光伏缺失功率实时预测方法、系统、介质及设备 |
-
2023
- 2023-04-21 CN CN202310437686.XA patent/CN116541681A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117293826A (zh) * | 2023-11-27 | 2023-12-26 | 山东大学 | 一种分布式光伏缺失功率实时预测方法、系统、介质及设备 |
| CN117293826B (zh) * | 2023-11-27 | 2024-04-05 | 山东大学 | 一种分布式光伏缺失功率实时预测方法、系统、介质及设备 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108761574B (zh) | 基于多源信息融合的降雨量估算方法 | |
| CN111651941B (zh) | 一种全球电离层电子总含量预测的算法 | |
| CN105740991B (zh) | 基于改进bp神经网络拟合多种气候模式的气候变化预测方法及系统 | |
| CN107885951B (zh) | 一种基于组合模型的水文时间序列预测方法 | |
| CN111242404B (zh) | 一种强降雨诱发洪灾事件的极端性评估方法及系统 | |
| CN111027175A (zh) | 基于耦合模型集成模拟的洪水对社会经济影响的评估方法 | |
| CN107545103A (zh) | 煤矿区土壤重金属含量空间模型建立方法 | |
| CN111581803B (zh) | 一种全球电离层电子含量的克里金代理模型构建方法 | |
| CN115423163A (zh) | 一种流域短期洪水事件预测方法、装置及终端设备 | |
| CN112131731B (zh) | 一种基于空间特征向量滤波的城市生长元胞模拟方法 | |
| CN106372277A (zh) | 森林立地指数时空估测中的变异函数模型优化方法 | |
| CN114254802B (zh) | 气候变化驱动下植被覆盖时空变化的预测方法 | |
| CN111401602A (zh) | 基于神经网络的卫星以及地面降水测量值同化方法 | |
| CN113704693B (zh) | 一种高精度的有效波高数据估计方法 | |
| CN111639803A (zh) | 一种应用于气候变化情景下区域未来植被指数的预估方法 | |
| CN116541681A (zh) | 一种基于协同克里金插值的复合灾害空间变异性识别方法 | |
| CN116910041B (zh) | 一种基于尺度分析的遥感降水产品的逐日订正方法 | |
| CN111879915B (zh) | 一种滨海湿地高分辨率的逐月土壤盐度监测方法及系统 | |
| CN105974495A (zh) | 利用分类拟合法预判目标区域未来平均云量的方法 | |
| CN111915158A (zh) | 一种基于Flood Area模型的暴雨灾害天气风险评估方法、装置及设备 | |
| CN114840616A (zh) | 一种基于时空插值的动态大气自然环境建模方法 | |
| CN113486295B (zh) | 基于傅里叶级数的臭氧总量变化预测方法 | |
| Wang et al. | Spatial Variation of Extreme Rainfall Observed From Two Century‐Long Datasets | |
| CN111597692B (zh) | 一种地表净辐射估算方法、系统、电子设备和存储介质 | |
| CN112949944A (zh) | 一种基于时空特征的地下水位智能预测方法及系统 |
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 |