CN114970187B - 一种实现水文气候时间序列趋势无偏估计的方法 - Google Patents
一种实现水文气候时间序列趋势无偏估计的方法 Download PDFInfo
- Publication number
- CN114970187B CN114970187B CN202210650692.9A CN202210650692A CN114970187B CN 114970187 B CN114970187 B CN 114970187B CN 202210650692 A CN202210650692 A CN 202210650692A CN 114970187 B CN114970187 B CN 114970187B
- Authority
- CN
- China
- Prior art keywords
- trend
- time sequence
- beta
- time series
- natural evolution
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 175
- 230000008569 process Effects 0.000 claims abstract description 82
- 238000012360 testing method Methods 0.000 claims abstract description 14
- 230000035772 mutation Effects 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims description 15
- 238000010586 diagram Methods 0.000 claims description 11
- 230000005923 long-lasting effect Effects 0.000 claims description 11
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 3
- 238000012353 t test Methods 0.000 claims description 3
- 230000001932 seasonal effect Effects 0.000 abstract 1
- 230000008859 change Effects 0.000 description 5
- 238000010998 test method Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003121 nonmonotonic effect Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000010485 coping Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000000528 statistical test Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/08—Probabilistic or stochastic CAD
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种实现水文气候时间序列趋势无偏估计的方法，包括：利用蒙特卡罗试验获取五种自然演变类型的时间序列差分后各自对应的一阶和二阶自相关系数的95%置信区间；剔除待分析时间序列的突变成分和季节成分，对剩余成分进行差分后求解一阶和二阶自相关系数，与上述各种类型对应的95%置信区间对比，确定其自然演变类型：若为白噪声，利用广义最小二乘法直接评估趋势显著性；若为AR过程，利用Newey‑West方法消除异方差和自相关性后评估趋势显著性；若为长持续过程，拟合并去除长持续成分后再评估趋势显著性。本发明针对不同自然演变类型采取针对性的方法检测趋势类型并评估显著性，消除自然演变特征的影响，实现对水文气候时间序列趋势的无偏估计。
Description
技术领域
本发明属于水文气候科学技术领域，具体指代一种实现水文气候时间序列趋势无偏估计的方法。
背景技术
准确揭示水文气候过程的演变特征、掌握其未来演变情势，是科学评估气候变化以及合理应对气候变化影响的基本依据和必要前提。趋势作为描述变量变化的一个重要指标，被广泛用于描述水文气候过程的长期演变特征。目前时间序列线性趋势识别方法主要有线性回归方法和Mann-Kendall(MK)趋势检验法。随着对水文气候过程演变特征认识的不断提高，研究者认识到水文气候过程受多种复杂因素的共同作用和影响，其演变特征往往表现出复杂的非线性趋势。然而，广泛应用的线性趋势只能以恒定速率描述其变化的平均速率，缺乏表征序列固有非线性趋势的能力。
为更好地描述水文气象要素演变特征的时变性，识别非线性和非单调趋势的方法被越来越多的应用到水文气候领域，如滑动平均法和基于数据拟合的指数拟合法、双曲线拟合法等。由于这些方法通常设置了特定的参数，导致获得的趋势具有主观性。近年来，基于时频分析的新技术和新方法被不断引入水文气候领域，如离散小波变换(Discretewavelet transform，DWT)能更有效地识别时间序列的非单调趋势，经验模态分解(Empiricmode decomposition，EMD)方法可依据时间序列本身的特性对序列进行分解，具有直接和自适应等突出特点。
在大量的实际分析计算过程中，线性趋势因其直观、易实现等优势被广泛应用于描述区域和全球水文气候要素的演变特征，因此仍然应用最为广泛。MK趋势检验法是非参数检验方法，相比于其他参数检验法，其数据不需要遵循一定的分布，受少数异常值的干扰较小，适用于分析水文气象要素等非正态分布数据。MK方法的前提假设为数据的独立性，但在实际分析中水文气候过程自然演变特征存在不同程度的自相关特性。理论上，当存在正的自相关特性时，出现较大值的时刻之后更有可能出现较大值，而小值后紧随小值，最终形成峰-谷交替结构。在这种情况下，水文气候过程的自然演变特征在某一时间段内会呈现趋势变化现象，对真实趋势识别的结果产生影响，且观测到的趋势也不能简单地归为人类活动等因素导致的外部强迫趋势。因此，评估水文气候过程的自然演变特征对趋势的影响，进而实现时间序列趋势的无偏估计十分重要。
目前，针对水文气候过程自然演变特征对趋势的影响，研究者提出了多种方法。当存在异方差或自相关性时，MK方法的方差估计是不准确的，从而影响对趋势显著性的统计检验。为了解决这一问题，White提出了Heteroskedasticity Consistent Covariances方法，当存在异方差时能够对协方差矩阵进行一致性估计，而无需知道异方差的形式，但该方法假定序列的残差是不存在自相关的。为此，Newey-West提出了一个更为普适性的估计量，当存在异方差和自相关性时仍然能对协方差矩阵进行一致性估计。同时，研究者进一步考虑了时间序列的长持续特性，首先消除时间序列中的长持续过程，并利用MK趋势检验法对剩余成分进行趋势识别。然而，由于不同的自然演变特征对趋势的影响方式和程度存在很大差异，导致目前仍然无法对水文气候过程的长期演变特征进行合理有效的评估。因此，准确判别水文气候过程的自然演变类型，是实现对时间序列趋势进行无偏估计的重要前提。
发明内容
针对于上述现有技术的不足，本发明的目的在于提供一种实现水文气候时间序列趋势无偏估计的方法，以解决现有技术中忽略时间序列自然演变特征对趋势识别的影响，导致无法准确识别水文气候时间序列趋势的类型与显著性的难题。
为达到上述目的，本发明采用的技术方案如下：
本发明的一种实现水文气候时间序列趋势无偏估计的方法，步骤如下：
1)分别生成与待分析时间序列S(t)长度相同的白噪声、AR(1)过程、AR(2)过程、单位根过程、长持续过程五种类型的时间序列，对生成的各时间序列进行差分处理后求解各自对应的一阶自相关系数和二阶自相关系数；
2)重复上述步骤1)，直至各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数的统计特征趋于稳定，进而获取各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间；
3)识别出时间序列S(t)中的突变成分B0，求解多年平均的季节成分S0，剔除时间序列S(t)的突变成分B0和季节成分S0，将剩余成分作为新时间序列S’(t)＝S(t)-B0-S0；
4)对新时间序列S’(t)做差分处理后，求解其一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)；
5)将一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)与步骤2)中得到的各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间进行对比，来判定时间序列S(t)的具体自然演变类型；
6)当时间序列S(t)的自然演变类型为白噪声过程时，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β1，并利用t检验方法评估其显著性：若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β1；
7)当时间序列S(t)的自然演变类型为AR(1)或AR(2)过程时，利用Newey-West方法处理该时间序列的异方差和自相关性后求得趋势斜率为β2；若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β2；
8)当时间序列S(t)的自然演变类型为单位根过程时，则时间序列S(t)呈现明显的随机性趋势，不存在确定性线性趋势，β＝0；
9)当时间序列S(t)的自然演变类型为长持续过程时，依次取不同的长持续特性di值，再利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)，得到剩余成分Si(t)＝S’(t)-Mi(t)；分别估计剩余成分Si(t)的趋势斜率βi与对应的残差Ri(t)；当残差最小时对应的剩余成分Si(t)趋势斜率记为β3，进而得到最终的趋势斜率的无偏估计结果β＝β3；
10)得到时间序列S(t)的趋势斜率无偏估计结果β，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β0；对比β0与β的差别，得到时间序列S(t)的自然演变特征对其趋势无偏估计结果的影响程度λ＝|β0-β|。
进一步地，所述的步骤1)具体包括：
11)利用蒙特卡罗方法生成白噪声的时间序列y1(t)；
12)利用一阶自回归模型生成AR(1)过程的时间序列y2(t)如下：
y2(t)＝ρ×y2(t-1)+u(t)
式中，t表示时序；ρ为一阶自相关系数，且|ρ|<1，u(t)是均值为0的符合独立同分布的白噪声序列；
13)利用二阶自回归模型生成AR(2)过程的时间序列y3(t)如下：
y3(t)＝ρ1×y3(t-1)+ρ2×y3(t-2)+u(t)
式中，ρ1和ρ2分别为一阶和二阶自相关系数，ρ1+ρ2<1，ρ2-ρ1<1，-1<ρ2<1；
14)生成单位根过程的时间序列y4(t)如下：
y4(t)＝y4(t-1)+u(t)；
15)利用ARFIMA模型生成长持续过程的时间序列y5(t)。
进一步地，所述步骤5)具体包括：
51)当AC_diff(1)和AC_diff(2)属于白噪声的95％置信区间内，则时间序列S(t)的自然演变类型判定为白噪声过程；
52)当AC_diff(1)和AC_diff(2)属于单位根过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为单位根过程；
53)当AC_diff(1)和AC_diff(2)属于AR(2)过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为AR(2)过程；
54)当AC_diff(1)和AC_diff(2)属于长持续过程或AR(1)过程的95％置信区间内，利用DFA方法求解新时间序列S’(t)的标度指数α；进一步判别S(t)的自然演变类型。
进一步地，所述步骤54)具体包括：
541)利用DFA方法获取新时间序列S’(t)的波动函数F(s)和时间尺度s的双对数散点图(ln(F(s))、ln(s))；
542)识别双对数散点图的结构突变点B1；
543)利用最小二乘法对区间B1<s<L/4的ln(F(s))和ln(s)进行线性拟合，线性趋势为标度指数α，L为时间序列S(t)的序列长度；
544)若α＝0.5，则时间序列S(t)的自然演变类型判定为AR(1)过程；
545)若α>0.5，则时间序列S(t)的自然演变类型判定为长持续过程。
进一步地，所述步骤9)具体包括：
91)设置不同的长持续特性di值，即di＝[-0.5:0.01:0.5]；
92)依次取不同的di值，利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)：
Mi(t)＝u(t)/(1-N)di
式中，N为自回归模型AR(n)的滞后算子，n为模型阶数；
93)得到剩余成分Si(t)＝S’(t)-Mi(t)；
94)利用广义最小二乘方法估计剩余成分Si(t)的线性趋势斜率βi；若通过显著性检验，则保持βi＝βi；若未通过显著性检验，则βi＝0；
95)得到趋势成分Ti(t)＝ci+βi×t，ci为趋势成分的截距项；
96)新时间序列S’(t)减去长持续成分Mi(t)和趋势成分Ti(t)后得到残差Ri(t)：
Ri(t)＝S’(t)-Mi(t)-Ti(t)；
97)当残差Ri(t)取最小值时，对应的斜率值为时间序列S(t)趋势斜率的无偏估计结果β＝β3。
本发明的有益效果：
本发明的方法首先考虑到水文气候过程自然演变特征对趋势的影响，通过利用蒙特卡罗试验确定各种自然演变类型统计特征的置信区间，以此为依据准确区分白噪声、单位根过程、AR1)过程、AR(2)过程和长持续过程，实现对时间序列自然演变类型的准确判别。其次，针对不同的自然演变类型采取针对性的方法对时间序列的趋势类型和显著性进行检测与评估，进而消除不同自然演变特征对趋势识别的影响，实现对水文气候时间序列趋势的无偏估计，从而为水文气候时间序列的模拟预测等工作提供科学依据。
附图说明
图1为本发明方法的流程图；
图2a为服从AR(1)过程的人工生成序列S11的示意图；
图2b为服从AR(1)过程的人工生成序列S12的示意图；
图2c为服从AR(1)过程的人工生成序列S13的示意图；
图2d为服从AR(1)过程的人工生成序列S14的示意图；
图3a为服从长持续过程的人工生成序列S21的示意图；
图3b为服从长持续过程的人工生成序列S22的示意图；
图3c为服从长持续过程的人工生成序列S23的示意图；
图3d为服从长持续过程的人工生成序列S24的示意图。
具体实施方式
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
参照图1所示，本发明的一种实现水文气候时间序列趋势无偏估计的方法，步骤如下：
1)分别生成与待分析时间序列S(t)长度相同的白噪声、AR(1)过程、AR(2)过程、单位根过程、长持续过程五种类型的时间序列，对生成的各时间序列进行差分处理后求解各自对应的一阶自相关系数和二阶自相关系数；具体包括：
11)利用蒙特卡罗方法生成白噪声的时间序列y1(t)；
12)利用一阶自回归模型生成AR(1)过程的时间序列y2(t)如下：
y2(t)＝ρ×y2(t-1)+u(t)
式中，t表示时序；ρ为一阶自相关系数，且|ρ|<1，u(t)是均值为0的符合独立同分布的白噪声序列；
13)利用二阶自回归模型生成AR(2)过程的时间序列y3(t)如下：
y3(t)＝ρ1×y3(t-1)+ρ2×y3(t-2)+u(t)
式中，ρ1和ρ2分别为一阶和二阶自相关系数，ρ1+ρ2<1，ρ2-ρ1<1，-1<ρ2<1；
14)生成单位根过程的时间序列y4(t)如下：
y4(t)＝y4(t-1)+u(t)；
15)利用ARFIMA模型生成长持续过程的时间序列y5(t)。
2)重复上述步骤1)，直至各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数的统计特征趋于稳定，进而获取各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间；
3)识别出时间序列S(t)中的突变成分B0，求解多年平均的季节成分S0，剔除时间序列S(t)的突变成分B0和季节成分S0，将剩余成分作为新时间序列S’(t)＝S(t)-B0-S0；
4)对新时间序列S’(t)做差分处理后，求解其一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)；
5)将一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)与步骤2)中得到的各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间进行对比，来判定时间序列S(t)的具体自然演变类型；具体包括：
51)当AC_diff(1)和AC_diff(2)属于白噪声的95％置信区间内，则时间序列S(t)的自然演变类型判定为白噪声过程；
52)当AC_diff(1)和AC_diff(2)属于单位根过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为单位根过程；
53)当AC_diff(1)和AC_diff(2)属于AR(2)过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为AR(2)过程；
54)当AC_diff(1)和AC_diff(2)属于长持续过程或AR(1)过程的95％置信区间内，利用DFA方法求解新时间序列S’(t)的标度指数α；进一步判别S(t)的自然演变类型。
其中，所述步骤54)具体包括：
541)利用DFA方法获取新时间序列S’(t)的波动函数F(s)和时间尺度s的双对数散点图(ln(F(s))、ln(s))；
542)识别双对数散点图的结构突变点B1；
543)利用最小二乘法对区间B1<s<L/4的ln(F(s))和ln(s)进行线性拟合，线性趋势为标度指数α，L为时间序列S(t)的序列长度；
544)若α＝0.5，则时间序列S(t)的自然演变类型判定为AR(1)过程；
545)若α>0.5，则时间序列S(t)的自然演变类型判定为长持续过程。
6)当时间序列S(t)的自然演变类型为白噪声过程时，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β1，并利用t检验方法评估其显著性：若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β1；
7)当时间序列S(t)的自然演变类型为AR(1)或AR(2)过程时，利用Newey-West方法处理该时间序列的异方差和自相关性后求得趋势斜率为β2；若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β2；
8)当时间序列S(t)的自然演变类型为单位根过程时，则时间序列S(t)呈现明显的随机性趋势，不存在确定性线性趋势，β＝0；
9)当时间序列S(t)的自然演变类型为长持续过程时，依次取不同的长持续特性di值，再利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)，得到剩余成分Si(t)＝S’(t)-Mi(t)；分别估计剩余成分Si(t)的趋势斜率βi与对应的残差Ri(t)；当残差最小时对应的剩余成分Si(t)趋势斜率记为β3，进而得到最终的趋势斜率的无偏估计结果β＝β3；具体包括：
91)设置不同的长持续特性di值，即di＝[-0.5:0.01:0.5]；
92)依次取不同的di值，利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)：
Mi(t)＝u(t)/(1-N)di
式中，N为自回归模型AR(n)的滞后算子，n为模型阶数；
93)得到剩余成分Si(t)＝S’(t)-Mi(t)；
94)利用广义最小二乘方法估计剩余成分Si(t)的线性趋势斜率βi；若通过显著性检验，则保持βi＝βi；若未通过显著性检验，则βi＝0；
95)得到趋势成分Ti(t)＝ci+βi×t，ci为趋势成分的截距项；
96)新时间序列S’(t)减去长持续成分Mi(t)和趋势成分Ti(t)后得到残差Ri(t)：
Ri(t)＝S’(t)-Mi(t)-Ti(t)；
97)当残差Ri(t)取最小值时，对应的斜率值为时间序列S(t)趋势斜率的无偏估计结果β＝β3。
10)得到时间序列S(t)的趋势斜率无偏估计结果β，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β0；对比β0与β的差别，得到时间序列S(t)的自然演变特征对其趋势无偏估计结果的影响程度λ＝|β0-β|。
示例中：
由于人工生成序列的趋势大小等情况已知，利用人工生成序列，有利于检验本发明方法的有效性。而实测水文气候时间序列的趋势显著性等情况往往未知，无法准确判断本发明方法识别时间序列趋势结果的准确性。为证明本发明方法识别时间序列趋势结果的准确性，设计方案时生成两类人工序列，分别用于验证自然演变特征对常规趋势识别方法的影响，以及本发明方法提高时间序列趋势识别结果的有效性。
第一类时间序列的自然演变特征为AR(1)过程，序列长度相同，但时间序列的自回归系数不同，分别为0.1、0.3、0.5、0.9，记为S11、S12、S13和S14(图2a，图2b，图2c，图2d)。第二类时间序列的自然演变特征为长持续过程，也具有相同的序列长度，但时间序列的长持续特性d的大小不同，分别为0.1、0.2、0.3、0.4，记为S21、S22、S23和S24(图3a，图3b，图3c，图3d)。选用MK方法、Newey-West方法和Local whittle方法对时间序列的趋势分别进行识别和评估，结果见表1(不同方法对人工生成时间序列的趋势识别结果)：
表1
趋势识别结果显示：对于AR(1)过程，MK方法无法准确识别时间序列的趋势；Localwhittle方法可以准确识别自相关性较弱时间序列(S11、S12、S13)的趋势，而无法识别自相关性较强时间序列(S14)的趋势；Newey-West方法可以消除时间序列中的异方差，进而准确识别时间序列的趋势。对于长持续过程，MK方法无法准确识别时间序列的趋势；Newey-West方法可以准确识别长持续特性较弱时间序列的趋势(S21、S22、S23)，而无法识别长持续特性较强时间序列(S24)的趋势；Local whittle方法可以消除序列中的长持续特性对趋势的影响，进而准确识别时间序列的趋势。上述结果表明，判别水文气候过程的自然演变特征类型，是准确识别时间序列趋势的重要前提。相比MK方法、Newey-West方法和Local whittle方法，本发明方法首先准确判别时间序列的自然演变类型，然后综合多种方法对时间序列的趋势类型进行准确检测与显著性评估，可以消除不同自然演变特征对趋势识别的影响，从而得到准确的时间序列趋势结果。
对比上述时间序列趋势识别结果，可以得到以下几点重要结论：(1)时间序列的自然演变特征对常规方法的趋势识别结果有显著影响；(2)较MK方法，Newey-West方法可以准确识别AR(1)过程时间序列的趋势，Local whittle方法可以准确识别长持续特性时间序列的趋势；(3)本发明方法通过准确判别待分析时间序列的自然演变类型，并综合多种方法对时间序列的趋势进行检测和识别，可以准确消除时间序列中的自然演变特征对趋势识别的影响，因此趋势识别结果更为可靠，可为准确揭示水文气候过程演变特征以及科学评估气候变化等提供方法支持。
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。
Claims (5)
1.一种实现水文气候时间序列趋势无偏估计的方法，其特征在于，步骤如下：
1)分别生成与待分析时间序列S(t)长度相同的白噪声、AR(1)过程、AR(2)过程、单位根过程、长持续过程五种类型的时间序列，对生成的各时间序列进行差分处理后求解各自对应的一阶自相关系数和二阶自相关系数；
2)重复上述步骤1)，直至各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数的统计特征趋于稳定，进而获取各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间；
3)识别出时间序列S(t)中的突变成分B0，求解多年平均的季节成分S0，剔除时间序列S(t)的突变成分B0和季节成分S0，将剩余成分作为新时间序列S’(t)＝S(t)-B0-S0；
4)对新时间序列S’(t)做差分处理后，求解其一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)；
5)将一阶自相关系数AC_diff(1)和二阶自相关系数AC_diff(2)与步骤2)中得到的各类型的时间序列差分处理后的一阶自相关系数和二阶自相关系数对应的95％置信区间进行对比，来判定时间序列S(t)的具体自然演变类型；
6)当时间序列S(t)的自然演变类型为白噪声过程时，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β1，并利用t检验方法评估其显著性：若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β1；
7)当时间序列S(t)的自然演变类型为AR(1)或AR(2)过程时，利用Newey-West方法处理该时间序列的异方差和自相关性后求得趋势斜率为β2；若通过显著性检验，则判定时间序列S(t)存在确定性趋势，趋势斜率的无偏估计结果β＝β2；
8)当时间序列S(t)的自然演变类型为单位根过程时，则时间序列S(t)呈现明显的随机性趋势，不存在确定性线性趋势，β＝0；
9)当时间序列S(t)的自然演变类型为长持续过程时，依次取不同的长持续特性di值，再利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)，得到剩余成分Si(t)＝S’(t)-Mi(t)；分别估计剩余成分Si(t)的趋势斜率βi与对应的残差Ri(t)；当残差最小时对应的剩余成分Si(t)趋势斜率记为β3，进而得到最终的趋势斜率的无偏估计结果β＝β3；
10)得到时间序列S(t)的趋势斜率无偏估计结果β，利用广义最小二乘方法估计时间序列S(t)的线性趋势斜率β0；对比β0与β的差别，得到时间序列S(t)的自然演变特征对其趋势无偏估计结果的影响程度λ＝|β0-β|。
2.根据权利要求1所述的实现水文气候时间序列趋势无偏估计的方法，其特征在于，所述的步骤1)具体包括：
11)利用蒙特卡罗方法生成白噪声的时间序列y1(t)；
12)利用一阶自回归模型生成AR(1)过程的时间序列y2(t)如下：
y2(t)＝ρ×y2(t-1)+u(t)
式中，t表示时序；ρ为一阶自相关系数，且|ρ|<1，u(t)是均值为0的符合独立同分布的白噪声序列；
13)利用二阶自回归模型生成AR(2)过程的时间序列y3(t)如下：
y3(t)＝ρ1×y3(t-1)+ρ2×y3(t-2)+u(t)
式中，ρ1和ρ2分别为一阶和二阶自相关系数，ρ1+ρ2<1，ρ2-ρ1<1，-1<ρ2<1；
14)生成单位根过程的时间序列y4(t)如下：
y4(t)＝y4(t-1)+u(t)；
15)利用ARFIMA模型生成长持续过程的时间序列y5(t)。
3.根据权利要求1所述的实现水文气候时间序列趋势无偏估计的方法，其特征在于，所述步骤5)具体包括：
51)当AC_diff(1)和AC_diff(2)属于白噪声的95％置信区间内，则时间序列S(t)的自然演变类型判定为白噪声过程；
52)当AC_diff(1)和AC_diff(2)属于单位根过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为单位根过程；
53)当AC_diff(1)和AC_diff(2)属于AR(2)过程的95％置信区间内，则时间序列S(t)的自然演变类型判定为AR(2)过程；
54)当AC_diff(1)和AC_diff(2)属于长持续过程或AR(1)过程的95％置信区间内，利用DFA方法求解新时间序列S’(t)的标度指数α；进一步判别S(t)的自然演变类型。
4.根据权利要求3所述的实现水文气候时间序列趋势无偏估计的方法，其特征在于，所述步骤54)具体包括：
541)利用DFA方法获取新时间序列S’(t)的波动函数F(s)和时间尺度s的双对数散点图(ln(F(s))、ln(s))；
542)识别双对数散点图的结构突变点B1；
543)利用最小二乘法对区间B1<s<L/4的ln(F(s))和ln(s)进行线性拟合，线性趋势为标度指数α，L为时间序列S(t)的序列长度；
544)若α＝0.5，则时间序列S(t)的自然演变类型判定为AR(1)过程；
545)若α>0.5，则时间序列S(t)的自然演变类型判定为长持续过程。
5.根据权利要求1所述的实现水文气候时间序列趋势无偏估计的方法，其特征在于，所述步骤9)具体包括：
91)设置不同的长持续特性di值，即di＝[-0.5:0.01:0.5]；
92)依次取不同的di值，利用长持续特性模型拟合新时间序列S’(t)的长持续成分Mi(t)：
Mi(t)＝u(t)/(1-N)di
式中，N为自回归模型AR(n)的滞后算子，n为模型阶数，u(t)是均值为0的符合独立同分布的白噪声序列；
93)得到剩余成分Si(t)＝S’(t)-Mi(t)；
94)利用广义最小二乘方法估计剩余成分Si(t)的线性趋势斜率βi；若通过显著性检验，则保持βi＝βi；若未通过显著性检验，则βi＝0；
95)得到趋势成分Ti(t)＝ci+βi×t，ci为趋势成分的截距项；
96)新时间序列S’(t)减去长持续成分Mi(t)和趋势成分Ti(t)后得到残差Ri(t)：
Ri(t)＝S’(t)-Mi(t)-Ti(t)；
97)当残差Ri(t)取最小值时，对应的斜率值为时间序列S(t)趋势斜率的无偏估计结果β＝β3。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210650692.9A CN114970187B (zh) | 2022-06-09 | 2022-06-09 | 一种实现水文气候时间序列趋势无偏估计的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210650692.9A CN114970187B (zh) | 2022-06-09 | 2022-06-09 | 一种实现水文气候时间序列趋势无偏估计的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114970187A CN114970187A (zh) | 2022-08-30 |
| CN114970187B true CN114970187B (zh) | 2024-04-16 |
Family
ID=82960890
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210650692.9A Active CN114970187B (zh) | 2022-06-09 | 2022-06-09 | 一种实现水文气候时间序列趋势无偏估计的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114970187B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009023660A1 (en) * | 2007-08-14 | 2009-02-19 | Nuance Communication, Inc. | Synthesis by generation and concatenation of multi-form segments |
| CN105069309A (zh) * | 2015-08-21 | 2015-11-18 | 中国科学院地理科学与资源研究所 | 一种识别水文时间序列非线性趋势的方法 |
| CN107066425A (zh) * | 2017-03-17 | 2017-08-18 | 中山大学 | 一种变化环境下超定量洪水非一致性分析方法 |
| CN110321518A (zh) * | 2019-06-14 | 2019-10-11 | 中国科学院地理科学与资源研究所 | 一种判定水文时间序列趋势类型的方法 |
| CN111611692A (zh) * | 2020-04-26 | 2020-09-01 | 武汉大学 | 气候变化情景下基于等可靠度的设计洪水推求方法及系统 |
| US10852439B1 (en) * | 2020-04-30 | 2020-12-01 | Beihang University | Global ionospheric total electron content prediction system |
| CN112149296A (zh) * | 2020-09-17 | 2020-12-29 | 中国科学院地理科学与资源研究所 | 一种判定水文时间序列平稳性类型的方法 |
-
2022
- 2022-06-09 CN CN202210650692.9A patent/CN114970187B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009023660A1 (en) * | 2007-08-14 | 2009-02-19 | Nuance Communication, Inc. | Synthesis by generation and concatenation of multi-form segments |
| CN105069309A (zh) * | 2015-08-21 | 2015-11-18 | 中国科学院地理科学与资源研究所 | 一种识别水文时间序列非线性趋势的方法 |
| CN107066425A (zh) * | 2017-03-17 | 2017-08-18 | 中山大学 | 一种变化环境下超定量洪水非一致性分析方法 |
| CN110321518A (zh) * | 2019-06-14 | 2019-10-11 | 中国科学院地理科学与资源研究所 | 一种判定水文时间序列趋势类型的方法 |
| CN111611692A (zh) * | 2020-04-26 | 2020-09-01 | 武汉大学 | 气候变化情景下基于等可靠度的设计洪水推求方法及系统 |
| US10852439B1 (en) * | 2020-04-30 | 2020-12-01 | Beihang University | Global ionospheric total electron content prediction system |
| CN112149296A (zh) * | 2020-09-17 | 2020-12-29 | 中国科学院地理科学与资源研究所 | 一种判定水文时间序列平稳性类型的方法 |
Non-Patent Citations (1)
| Title |
|---|
| 水文时间序列趋势和跳跃分析的再抽样方法研究;刘攀;郭生练;肖义;李玮;郭富强;;水文;20070425(第02期);52-56 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114970187A (zh) | 2022-08-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gerstengarbe et al. | Estimation of the beginning and end of recurrent events within a climate regime | |
| Wu et al. | A prediction method using the grey model GMC (1, n) combined with the grey relational analysis: a case study on Internet access population forecast | |
| CN110321518B (zh) | 一种判定水文时间序列趋势类型的方法 | |
| CN104280526A (zh) | 水质自动在线监测设备测量误差的分析和估计方法 | |
| CN112149296B (zh) | 一种判定水文时间序列平稳性类型的方法 | |
| CN102488517A (zh) | 一种检测脑电信号中爆发抑制状态的方法以及装置 | |
| Syrjala | Critique on the use of the delta distribution for the analysis of trawl survey data | |
| EP3649568B1 (en) | Method for automatic detection of physical modes in a modal analysis model | |
| Huang et al. | Time series modeling and filtering method of electric power load stochastic noise | |
| CN117434153A (zh) | 基于超声波技术的道路无损检测方法及系统 | |
| CN110852906A (zh) | 一种基于高维随机矩阵进行窃电嫌疑识别的方法和系统 | |
| Han et al. | An online approach for intracranial pressure forecasting based on signal decomposition and robust statistics | |
| CN108664452B (zh) | 一种室内人员数量的确定方法及确定系统 | |
| CN105046372B (zh) | 一种逐日蔬菜价格的预测方法和装置 | |
| CN114970187B (zh) | 一种实现水文气候时间序列趋势无偏估计的方法 | |
| CN110837088B (zh) | 一种星载激光测高仪数据去噪方法 | |
| Varshavskiy et al. | Efficiency estimation of the noise digital filtering algorithms | |
| Farrelly et al. | Determination of uncertainty in environmental noise measurements by bootstrap method | |
| Stadnytska et al. | Analyzing fractal dynamics employing R | |
| Warsza | Evaluation of the type A uncertainty in measurements with autocorrelated observations | |
| CN114840802B (zh) | 一种判别水文气候过程自然演变类型的方法 | |
| CN113537156A (zh) | 一种基于区间标准差结合频谱分析的振动数据清洗方法 | |
| CN111210278A (zh) | 一种基于时间序列的煤炭行业股价预测方法 | |
| Kabović et al. | The influence of the input parameters variation of the non-seasonal ARIMAX model on the accuracy of meteorological parameters forecasting | |
| CN111915858A (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |