CN111310109A - 一种基于vmd-arma-garch模型的非良态风速建模方法 - Google Patents

一种基于vmd-arma-garch模型的非良态风速建模方法 Download PDF

Info

Publication number
CN111310109A
CN111310109A CN202010174623.6A CN202010174623A CN111310109A CN 111310109 A CN111310109 A CN 111310109A CN 202010174623 A CN202010174623 A CN 202010174623A CN 111310109 A CN111310109 A CN 111310109A
Authority
CN
China
Prior art keywords
wind speed
model
arma
garch
vmd
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.)
Granted
Application number
CN202010174623.6A
Other languages
English (en)
Other versions
CN111310109B (zh
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 Railway Eryuan Engineering Group Co Ltd CREEC
Original Assignee
China Railway Eryuan Engineering Group Co Ltd CREEC
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 Railway Eryuan Engineering Group Co Ltd CREEC filed Critical China Railway Eryuan Engineering Group Co Ltd CREEC
Priority to CN202010174623.6A priority Critical patent/CN111310109B/zh
Publication of CN111310109A publication Critical patent/CN111310109A/zh
Application granted granted Critical
Publication of CN111310109B publication Critical patent/CN111310109B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明公开了一种基于VMD‑ARMA‑GARCH模型的非良态风速建模方法，包括以下步骤：步骤一：建立非良态风速方程，脉动风速方程，脉动风速的功率谱函数等；步骤二：采用判别准则并结合VMD分解法获取时变平均风速
Figure DDA0002410365770000011
和脉动风速u(t)；采用ARMA‑GARCH模型获取脉动风速强度包络函数σu(t)；采用ARMA模型/AR模型获取归一化平稳风速功率谱函数Sα(ω)。本发明提供了一种普遍适用于各种不同类型非良态风速的统一建模方法，通过VMD分解法，ARMA‑GARCH模型和ARMA模型/AR模型，解决了非良态风速精确建模的三个关键问题，即：时变平均风速的提取、脉动风速强度包络函数的计算、脉动风速时变功率谱的估算，从而能够准确建立具有明确物理意义的非良态风速模型，且三种方法协同使用，能大幅减少计算工作量，提高建模精确性。

Description

一种基于VMD-ARMA-GARCH模型的非良态风速建模方法
技术领域
本发明涉及一种风速建模方法，特别是一种基于VMD-ARMA-GARCH模型的非良态风速建模方法。
背景技术
大量实测风速资料表明：非良态风荷载具有很强的非平稳(瞬态)特性，将会对结构物造成严重的破坏，例如，雷暴风、台风、龙卷风、飓风等。一般来说，假定当阵风中平稳部分的持时远大于大型结构的基本振动周期时，那么将这段风速简化处理为平稳随机过程后进行经典风致响应分析是合理的，否则将会引起计算结果的偏差。为了更加真实地反映结构的风致振动性能，对非良态风速进行精确的数学建模是开展结构风致响应分析的前提。
大气边界层(atmospheric-boundary-layer,简称ABL)的良态风速模型可表示为常值平均风速与平稳脉动风速的叠加。其中，常值平均风速可由10min或1h时间窗宽对原始风速进行总体平均得到。与ABL截然不同，非良态风速模型则表示为确定性的时变平均风速与非平稳脉动风速的叠加。精确的非平稳风场模型需要从三个方面进行考虑，(1)时变平均风速的提取；(2)脉动风速强度包络函数的计算；(3)脉动风速时变功率谱的估算。
在工程中，获得多个非良态随机样本并非易事，往往仅能够获得单个观测样本。更重要的是，非良态风速随机过程具有非各态历经的特点，不能够完全照搬良态平稳风速的统计建模思路，而且也很难通过统一的数学表达式对各种不同类型的非良态风速进行拟合。然而，近年来，自然界非良态极端灾害性大风事件频发，严重威胁着建筑结构(例如，大跨桥梁结构、超高层建筑、空间复杂结构物)的安全性，而且现阶段尚未得到极端非良态风速精确的数学模型。因此，为了掌握极端风速的变化规律以及避免其对结构物产生的不利影响，发展一种普遍适用于各种不同类型的非良态风速建模方法具有重大的工程价值。
发明内容
本发明的目的在于：针对现有技术存在的问题，提供一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，解决非良态风速精确建模的三个关键问题，即：时变平均风速的提取、脉动风速强度包络函数的计算、脉动风速时变功率谱的估算。
为了实现上述目的，本发明采用的技术方案为：
一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，包括以下步骤：
步骤一：建立非良态风速方程
Figure BDA0002410365750000021
式中，U(t)为原始风速，
Figure BDA0002410365750000022
为确定性的时变平均风速，u(t)为均匀调制的随机脉动风速，t为时间；
所述脉动风速u(t)表示为与时间t相关的强度包络函数σu(t)与归一化平稳风速α(t)的乘积，即u(t)＝σu(t)α(t)；
所述脉动风速的功率谱函数Su(ω,t)表示为强度包络函数平方
Figure BDA0002410365750000023
与归一化平稳风速功率谱函数Sα(ω)的乘积，即
Figure BDA0002410365750000024
空间x1处的脉动风速u(x1,t)与x2处的u(x2,t)之间的互谱，表示为：
Figure BDA0002410365750000025
式中：Su(x1,x2,ω,t)为脉动风速u(x1,t)与u(x2,t)之间的互谱，Su(x1,ω,t)为x1处的脉动风速的功率谱函数，Su(x2,ω,t)为x2处的脉动风速的功率谱函数，coh(x1,x2,ω,t)为脉动风速间的相干函数，ω为频率；
步骤二：采用提出判别准则并结合VMD分解法(Variational ModeDecomposition，变分模态分解法)获取所述时变平均风速
Figure BDA0002410365750000031
和所述脉动风速u(t)；
采用ARMA-GARCH模型(Autoregressive moving average-GeneralizedAutoRegressive Conditional Heteroskedasticity model，自回归移动平均—广义自回归条件异方差模型)获取所述强度包络函数σu(t)；
采用ARMA模型(Autoregressive moving average model，自回归移动平均模型)或AR模型(Autoregressive model，自回归模型)获取所述归一化平稳风速功率谱函数Sα(ω)。
本发明提供了一种普遍适用于各种不同类型的非良态风速(例如雷暴风、台风、龙卷风等)建模方法，通过VMD分解法，ARMA-GARCH模型和ARMA模型/AR模型，解决了非良态风速精确建模的三个关键问题，即：时变平均风速的提取、脉动风速强度包络函数的计算、脉动风速时变功率谱的估算，从而能够准确建立具有明确物理意义的非良态风速模型，且三种方法协同使用，能大幅减少计算工作量，提高建模精确性。
作为本发明的优选方案，所述步骤二中提出的判别准则为：当平均风速中最高频率fmax取值为结构基频f1的1/5～1/10时，可忽略结构的动力效应，即，fmax∈[0.1f1,0.2f1]，所以当模态分量的中心频率fc大于0.2f1时，则认为其为脉动风成分，如下式所示：
Figure BDA0002410365750000032
作为本发明的优选方案，所述步骤二中，采用VMD分解法获取所述时变平均风速
Figure BDA0002410365750000033
和所述脉动风速u(t)，包括以下步骤：
步骤A21：初始化VMD分解法的参数，包括保真度系数τ，第一个中心频率更新参数DC，中心频率初始化参数init，收敛准则阈值ε和惩罚因子α，其中，模态分量层数K＝1；
步骤A22:采用所述步骤A21中设定的参数对原始风速信号进行第一层分解，并获得第一层模态分量c1的中心频率；
步骤A23:将所述模态分量层数K值加1，其余参数保持不变，重复所述步骤A22，对所述原始风速信号进行K层模态分解，并获得第K层模态分量ck的中心频率；
步骤A24:判断所述第K层模态分量的中心频率ck是否大于0.2f1，所述f1为结构的基频，
如果“是”，则第1层～第K-1层模态分量的叠加视为时变平均风速，并进入步骤A25；
如果“否”，则返回步骤A23；
步骤A25:将所述原始风速信号减去所述时变平均风速
Figure BDA0002410365750000041
获得所述脉动风速u(t)。
作为本发明的优选方案，所述步骤A21中，初始化VMD分解法的参数时，设定保真度系数τ＝0，第一个中心频率更新参数DC＝0，中心频率初始化参数init＝1，收敛准则阈值ε＝1e-7，惩罚因子α＝2000。
作为本发明的优选方案，所述步骤二中，采用ARMA-GARCH模型获取所述强度包络函数σu(t)时，
所述脉动风速u(t)的均值方程表示为：
φ(B)u(t)＝θ(B)ε(t)，ε(t)＝σε(t)μ(t)；
所述脉动风速u(t)的方差方程表示为：
Figure BDA0002410365750000051
式中，φ(B)为p阶AR(p)多项式，表示为φ(B)＝1+a1B+a2B2+…apBp，a＝[a1,a2,…,ap]T为AR模型的参数向量，p为AR模型的阶数，θ(B)为q阶MA(q)多项式，表示为θ(B)＝1+b1B+b2B2+…bqBq，b＝[b1,b2,…,bq]T为MA模型的参数向量，q为MA模型的阶数，B为向后移位算子，B[u(t)]＝u(t-1)，
u(t)为脉动风速，ε(t)为零均值的异方差序列，σε(t)为ε(t)的时变标准差，μ(t)为均值为0、方差为1的独立同分布随机变量序列，ηi、λj、γ为GARCH模型的参数，i＝1,2,…,m，j＝1,2,…,l，m和l为GARCH模型的阶数。
作为本发明的优选方案，进一步地，通过推导所述强度包络函数σu(t)的解析式表示为：
Figure BDA0002410365750000052
式中，[G1，G2，…]为格林函数，通过待定系数法，可得到格林函数的计算公式。
进一步地，由于脉动风速u(t)的均值近似为0，在GARCH模型公式中所描述的波动特征，可近似认为脉动风速u(t)与其残差的包络函数相一致。
作为本发明的优选方案，采用ARMA-GARCH模型获取所述强度包络函数σu(t)，包括以下步骤：
步骤B21：设置ARMA模型参数p、q的取值范围，其中，p＝1,2,3,…,pmax，q＝0,1,2,…,qmax；
步骤B22：ARMA模型定阶，包括：
步骤B221：遍历ARMA模型参数p、q的取值范围，形成共计pmax·(qmax+1)组ARMA模型参数对；
步骤B222：针对每一组p、q的取值，利用Ljung-Box Q测试方法计算残差序列εpq的自相关性，并计算相应的统计值Qtest；
步骤B223：选择Qtest最小的一组p、q组合作为ARMA模型的最优参数估计；
步骤B23：GARCH模型定阶，包括通过最优原则比较，基于t分布的GARCH(1,1)模型效果最优，采用最大似然估计方法确定GARCH项的参数值；
步骤B24：根据所述的方差方程计算残差序列εpq的时变标准差；
步骤B25：利用所述残差序列εpq的时变标准差，推算脉动风速的所述强度包络函数σu(t)。
作为本发明的优选方案，所述步骤B21中，pmax＝20，qmax＝pmax-1。
作为本发明的优选方案，所述步骤二中，采用ARMA模型或AR模型获取所述归一化平稳风速功率谱函数Sα(ω)时，
归一化平稳风速α(t)表示为：
φ(B)α(t)＝θ(B)e(t)；
式中，φ(B)为p阶AR(p)多项式，表示为φ(B)＝1+a1B+a2B2+…apBp，a＝[a1,a2,…,ap]T为AR模型的参数向量，p为AR模型的阶数，θ(B)为q阶MA(q)多项式，表示为θ(B)＝1+b1B+b2B2+…bqBq，b＝[b1,b2,…,bq]T为MA模型的参数向量，q为MA模型的阶数，B为向后移位算子，B[u(t)]＝u(t-1)，e(t)为高斯白噪声序列且方差为常值
Figure BDA0002410365750000061
根据N个已知的观测数据α(0)，α(1)，…，α(N-1)，归一化的脉动风谱表示为：
Figure BDA0002410365750000071
根据Wold分解定理可知，一个ARMA模型可以用一个阶数足够大的AR模型来代替。通过Cadzow谱估计子线性化处理，可将功率谱表达式转换成只需要AR参数，而不需要MA具体参数值(即参数b1，b2，…，bq全为0)的计算表达式，当采用AR模型时，进一步简化为：
Figure BDA0002410365750000072
式中，△t为采样间隔。
作为本发明的优选方案，所述AR模型阶数采用线性代数法或信息量准则进行确定，所述线性代数法包括奇异值分解法、Gram-Schmidt正交法，所述信息量准则包括FPE准则、AIC准则，所述AR模型参数采用最小二乘方法、Yule-Walker方法进行求解。
作为本发明的优选方案，采用奇异值分解法结合总体最小二乘方法进行求解时，包括以下步骤：
C21：用样本的相关矩阵R代替增广矩阵B，计算B的SVD，并储存奇异值σ11≥σ11…≥σhh≥0和矩阵V；
C22：采用归一化比值法，确定增广矩阵B的有效秩p，得到ARMA模型的AR模型阶数估计；
C23：计算矩阵S(p)，
Figure BDA0002410365750000073
Figure BDA0002410365750000074
为矩阵V第j列的加窗段，上标H表示矩阵共轭转置；
C24：计算S(p)的逆矩阵S-(p)，并由ai＝S-(p)(i+1,1)/S-(P)(1,1),i＝1,…,p计算出待求AR模型参数的总体最小二乘估计值。
综上所述，由于采用了上述技术方案，本发明的有益效果是：
(1)本发明提供一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，从数据自身特征出发，具有普遍性，可适用于不同类型的非良态风速样本，填补了这一方面研究的空白。
(2)基于VMD分解法对原始风速的多个分量分解结果，提出了提取时变平均风速的判别准则，即当分量中心频率>0.2f1时，则认为其为脉动风成分，相反地，分量中心频率≤0.2f1时，则认为其为平均风速成分。应用该判别准则不会导致低估由脉动风速产生的结构动力响应以及总体风致响应，具有较好的工程应用价值。
(3)ARMA-GARCH模型在非平稳风速建模方面具有显著的优势，突破了已有计算方法基于“局部平均”假设的束缚，能够获得精确的强度包络函数。
(4)计算的强度包络函数能够反映非良态风速的非平稳“强、弱”程度，可将弱非平稳风速简化为平稳风速处理，可极大提高风致响应计算效率。
(5)ARMA模型(或简化的AR模型)能够提供高分辨率、低方差的归一化脉动风速时变功率谱估计结果，有助于提高非良态风谱的估算精度。
(6)非良态脉动风速时变功率谱的数学模型可表示为可分离谱的形式，即时变功率谱写为强度包络函数与归一化脉动风速功率谱的乘积。该分离谱的形式不仅具有明确的物理意义，而且符合随机振动理论中关于非平稳谱的相关定义，可便于揭示非良态风速对复杂大跨结构的瞬态作用机理。
(7)所述建模方法中，三种方法协同使用，能大幅减少计算工作量，提高建模精确性，不仅可掌握非良态风速的数学模型、演变规律，而且可为非良态风速与大跨桥梁相互作用的随机振动分析，提供时域和频域两方面的荷载输入。
附图说明
图1是本发明所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法的流程图。
图2是本发明实施例1所述的大跨悬索桥的结构示意图。
图3是本发明实施例1所述的原始雷暴风速及其时变平均风速的示意图。
图4是本发明实施例1所述的脉动风速及其强度包络函数的示意图。
图5是本发明实施例1所述的归一化脉动风速的示意图。
图6是本发明实施例1所述的归一化风谱的示意图。
图7是本发明实施例1所述的时变风谱的示意图。
图8是本发明实施例2所述的香港某大桥的结构示意图。
图9是本发明实施例2所述的原始台风“杜鹃”样本的风速示意图。
图10是本发明实施例2所述的时变平均风速的示意图。
图11是本发明实施例2所述的脉动风速及其强度包络函数的示意图。
图12是本发明实施例2所述的时变风谱的示意图。
图13是本发明实施例3中采用本发明所述的建模方法得到的时变谱图。
图14是本发明实施例3中采用Priestley方法得到的时变谱图。
具体实施方式
下面结合附图，对本发明作详细的说明。
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
以大跨度桥梁结构为例，对原始风速信号进行三维矢量分解，获得顺风向风速U(t)、顺桥向风速v(t)和竖向风速w(t)，以顺风向风速U(t)为例。
如图1所示，一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，包括以下步骤：
步骤一：建立非良态风速方程：
Figure BDA0002410365750000101
式中，U(t)为原始风速，
Figure BDA0002410365750000102
为确定性的时变平均风速，u(t)为均匀调制的随机脉动风速，t为时间；
其中，平均风速随时间连续变化并且能够刻画原始风速的变化趋势；根据已有研究成果表明，风速中内蕴的频率成分随时间的变化较弱，而描述风速非平稳特征主要体现在强度信息的瞬态变化，所述脉动风速u(t)表示为与时间t相关的强度包络函数σu(t)与归一化平稳风速α(t)的乘积，即
u(t)＝σu(t)α(t) (2)
相应地，所述脉动风速的功率谱函数Su(ω,t)表示为强度包络函数平方σu2(t)与归一化平稳风速功率谱函数Sα(ω)的乘积，即
Figure BDA0002410365750000103
不同桥梁主梁跨长坐标x1、x2位置处，x1处的脉动风速u(x1,t)与x2处的u(x2,t)之间的互谱，表示为：
Figure BDA0002410365750000104
式中：Su(x1,x2,ω,t)为脉动风速u(x1,t)与u(x2,t)之间的互谱，Su(x1,ω,t)为x1处的脉动风速的功率谱函数，Su(x2,ω,t)为x2处的脉动风速的功率谱函数，coh(x1,x2,ω,t)为脉动风速间的相干函数，ω为频率。
而顺桥向风速v(t)和竖向风速w(t)可完全参照顺风向风速U(t)的研究方法进行建模，不再赘述。
步骤二：采用提出的判别准则并结合VMD分解法获取所述时变平均风速
Figure BDA0002410365750000111
和所述脉动风速u(t)。
目前，时-频分析方法已成为处理非平稳、非线性信号的有力工具，例如，小波类方法(小波分解法、平稳小波变换法)，模态分解类方法(经验模态分解法、总体经验模态分解法)。这些时-频分析方法均可将任意频带宽度的风速时程分解为一系列子分量的叠加。然而，选择上述不同的方法以及确定各类方法中的未知参数很大程度上依赖于研究者的经验，并非基于数据驱动的完全自适应方法，而且即使采用完全一致的风速数据也会得到不一致的结果。
VMD分解法是一种自适应信号的、准正交的、完全非递归的分解方法。该方法以经典维纳滤波、Hilbert变换及混频问题为核心，通过迭代搜寻变分模型最优解，来确定信号中各个模态分量的中心频率和有限带宽，能够自适应地实现信号中各个模态分量的分离。相比于上述时-频分析方法，VMD分解法具有较好地筛选密集模态以及抑制边界效应的能力，并且该方法是从长周期的低频分量至短周期高频分量依次进行分解，非常适合用于风工程中对于时变平均风速的提取。
原始风速经VMD分解后，哪些分量可归入平均风速？哪些分量可纳入脉动风成分？相比于结构的基频f1，倘若平均风速中最高频率足够小，那么平均风速对结构的动力放大效应可以忽略不计，仅需进行简单的拟静力分析。应用这一判别准则可保证平均风速成分中不含脉动风速，这样则不会导致低估由脉动风速产生的结构动力响应以及总体风致响应(平均风速的静风响应组合脉动风速的抖振响应)。由结构动力学可知，当平均风速中最高频率fmax取值为结构基频f1的1/5～1/10时，可忽略动力效应，即，fmax∈[0.1f1,0.2f1]，当VMD分解后的分量中其中心频率fc大于0.2f1时，则认为其为脉动风成分，如下式所示：
Figure BDA0002410365750000121
采用VMD分解法获取时变所述时变平均风速
Figure BDA0002410365750000123
和所述脉动风速u(t)，包括以下步骤：
步骤A21：初始化VMD分解法的6个参数，包括保真度系数τ＝0，第一个中心频率更新参数DC＝0，中心频率初始化参数init＝1，收敛准则阈值ε＝1e-7和惩罚因子α＝2000，其中，模态分量层数K＝1；
步骤A22:采用所述步骤S21中设定的参数对原始风速信号进行第一层分解，并获得第一层模态分量c1的中心频率；
步骤A23:将所述模态分量层数K值加1，其余5个参数保持不变，重复所述步骤A22，对所述原始风速信号进行K层模态分解，并获得第K层模态分量ck的中心频率；
步骤A24:判断所述第K层模态分量的中心频率ck是否大于0.2f1，所述f1为结构的基频，
如果“是”，则第1层～第K-1层模态分量的叠加视为时变平均风速，并进入步骤A25；
如果“否”，则返回步骤A23；
步骤A25:将所述原始风速信号减去所述时变平均风速
Figure BDA0002410365750000122
获得所述脉动风速u(t)。
步骤三：采用ARMA-GARCH模型获取所述强度包络函数σu(t)。
基于单个风速样本已有计算强度包络函数的方法包括移动窗口加权平均法、Kalman滤波法、Kernel回归法等。目前上述这些计算方法均存在缺陷，其中，移动窗口加权平均法和Kernel回归法会面临选择窗函数以及窗宽的困难。Kalman滤波法则会遇到初始参数选取以及计算耗时的难题。此外，更为重要的是，上述三类方法均束缚于“局部平均”的假设，即样本特征随时间的变化缓慢。显然，已有方法适用于自相关性较弱或者样本能力集中在高频区域的序列样本。由于风速样本的能量聚集在低频区域，所以精确计算风速强度包络函数的方法仍需要深入研究。
本发明提出解决该问题的方案，基于自回归移动平均(ARMA)—广义自回归条件异方差(GARCH)模型在经济学金融领域中分析波动率的出色性能，提出采用ARMA(p,q)-GARCH(m,l)模型计算脉动风速的强度包络函数。
对于脉动风速u(t)，采用ARMA(p,q)-GARCH(m,l)对其进行建模，所述脉动风速u(t)的均值方程表示为：
φ(B)u(t)＝θ(B)ε(t) (6)
ε(t)＝σε(t)μ(t) (7)
所述脉动风速u(t)的方差方程表示为：
Figure BDA0002410365750000131
式中，φ(B)为p阶AR(p)多项式，表示为φ(B)＝1+a1B+a2B2+…apBp，a＝[a1,a2,…,ap]T为AR模型的参数向量，p为AR模型的阶数，θ(B)为q阶MA(q)多项式，表示为θ(B)＝1+b1B+b2B2+…bqBq，b＝[b1,b2,…,bq]T为MA模型的参数向量，q为MA模型的阶数，B为向后移位算子，B[u(t)]＝u(t-1)，
u(t)为脉动风速，ε(t)为零均值的异方差序列，σε(t)为ε(t)的时变标准差，μ(t)为均值为0、方差为1的独立同分布随机变量序列，ηi、λj、γ为GARCH模型的参数，i＝1,2,…,m，j＝1,2,…,l，m和l为GARCH模型的阶数。
本质上，式(6)～式(8)所示ARMA-GARCH模型认为时间序列每个时间点的异方差值是最近l个时间点残差平方的线性组合，与最近m个时间点异方差序列σε(t)的线性组合叠加后得到。
进一步地，对公式(6)进行转换，所述强度包络函数σu(t)的解析式可记为：
Figure BDA0002410365750000141
式中，[G1，G2，…]为格林函数，通过待定系数法，可容易得到格林函数的计算公式，此处不再详述。
进一步地，由于脉动风速u(t)的均值近似为0，在GARCH模型公式(8)中所描述的波动特征，可近似认为脉动风速u(t)与其残差的包络函数相一致。
采用ARMA-GARCH模型获取所述强度包络函数σu(t)，包括以下步骤：
步骤B21：设置ARMA模型参数p、q的取值范围，其中，p＝1,2,3,…,pmax，q＝0,1,2,…,qmax，pmax＝20，qmax＝pmax-1；
步骤B22：ARMA模型定阶，包括：
步骤B221：遍历ARMA模型参数p、q的取值范围，形成共计pmax·(qmax+1)组ARMA模型参数对；
步骤B222：针对每一组p、q的取值，利用Ljung-Box Q测试方法计算残差序列εpq的自相关性，并计算相应的统计值Qtest；
步骤B223：统计值Qtest越小表明残差序列εpq的自相关性越弱，选择Qtest最小的一组p、q组合作为ARMA模型的最优参数估计；
步骤B23：GARCH模型定阶，为了简洁，由于GARCH(1,1)模型(m和l取值为1)描述异方差性简洁并且拟合效果好，通过最优原则比较，基于t分布的GARCH(1,1)模型效果最优，采用最大似然估计方法确定GARCH项的参数值；
步骤B24：根据公式(8)计算残差序列εpq的时变标准差；
步骤B25：由于脉动风速的强度包络函数与其ARMA模型残差序列的时变标准差成比例，利用所述残差序列εpq的时变标准差，推算脉动风速的所述强度包络函数σu(t)。
步骤四：采用ARMA模型或AR模型获取所述归一化平稳风速功率谱函数Sα(ω)。
平稳随机信号功率谱的估计主要包括经典谱估计和现代谱估计两种方法，经典谱估计方法，例如基于Fourier变换的周期图方法、自相关函数法，由于输入信号为有限长，存在谱估计为有偏估计、分辨率低、方差性能差、旁瓣泄漏等缺陷。现代谱估计方法包括AR模型、MA模型、ARMA模型等，可将平稳随机过程看成是白噪声激励线性时不变系统的输出。由于该类谱估计方法不受“不确定性原理”的束缚，频谱具有高分辨率与低方差(光滑)的优点，非常适用于数据点较短的时间序列。
本发明基于强度包络线计算时对公式(6)所述模型的研究成果以及ARMA模型广泛的代表性和实用性，为了便捷，采用ARMA模型方法对归一化的平稳脉动风速α(t)进行功率谱估计。
给予公式(6)，归一化平稳风速α(t)表示为：
φ(B)α(t)＝θ(B)e(t) (10)
式中，e(t)为高斯白噪声序列且方差为常值
Figure BDA0002410365750000151
其余参数含义与上述相同。
根据N个已知的观测数据α(0)，α(1)，…，α(N-1)，归一化的脉动风谱表示为：
Figure BDA0002410365750000152
根据Wold分解定理可知，一个ARMA模型可以用一个阶数足够大的AR模型来代替。通过Cadzow谱估计子线性化处理，可将功率谱表达式转换成只需要AR参数，而不需要MA具体参数值(即参数b1，b2，…，bq全为0)的计算表达式，当采用AR模型时，进一步简化为：
Figure BDA0002410365750000161
式中，△t为采样间隔，其余参数含义同上。
相比于ARMA模型不仅需要估计AR参数和MA参数(MA参数估计必须求解非线性方程组)，AR模型相对简单，故工程上可采用AR模型进行代替。其中，AR模型阶数可采用线性代数法(奇异值分解法、Gram-Schmidt正交法)、信息量准则(FPE准则、AIC准则等)进行确定，模型参数可选用最小二乘方法、Yule-Walker方法等进行求解。在上述方法中，奇异值分解法(SVD)结合总体最小二乘方法(TLS)，可用较少的AR阶数实现高分辨率的风谱估计。
本发明采用奇异值分解法(SVD)结合总体最小二乘方法(TLS)进行求解时，包括以下步骤：
C21：用样本的相关矩阵R代替增广矩阵B，计算B的SVD，并储存奇异值σ11≥σ11…≥σhh≥0和矩阵V；
C22：采用归一化比值法，确定增广矩阵B的有效秩p，得到ARMA模型的AR模型阶数估计；
C23：计算矩阵S(p)，
Figure BDA0002410365750000162
Figure BDA0002410365750000163
为矩阵V第j列的加窗段，上标H表示矩阵共轭转置；
C24：计算S(p)的逆矩阵S-(p)，并由ai＝S-(p)(i+1,1)/S-(P)(1,1),i＝1,…,p计算出待求AR模型参数的总体最小二乘估计值。
本发明提供了一种普遍适用于各种不同类型的非良态风速(例如，雷暴风、台风、龙卷风等)建模方法，通过VMD分解法，ARMA-GARCH模型和ARMA模型/AR模型，解决了非良态风速精确建模的三个关键问题，即：时变平均风速的提取、脉动风速强度包络函数的计算、脉动风速时变功率谱的估算，从而能够准确建立具有明确物理意义的非良态风速模型，且三种方法协同使用，能大幅减少计算工作量，提高建模精确性。
实施例1
以我国西南山区主跨为628m大跨悬索桥(如图2所示)桥址处观测得到的一组山区雷暴风速为例进行说明。其中，该大桥的基频f1＝0.13Hz。图3所示为原始雷暴风速及其时变平均风速、10min常值平均风速，显然常值平均风速不能够刻画雷暴风速样本的变化趋势。VMD分解获取时变平均风速时，分解层数K＝3。可以看出，雷暴风速在10min时间段内存在一次明显的幅值“上升段”和“下降段”。图4所示为脉动风速及其强度包络函数结果，图5～图6所示为归一化脉动风速及其归一化平稳风速功率谱函数，图7所示为雷暴风的时变风谱结果。
实施例2
以主跨为1018m香港某大桥(如图8所示，图8中实心圆点代表风速仪安装位置)桥址处观测得到的近海台风“杜鹃”(如图9所示)样本为例进行说明。其中，该大桥的基频为0.161Hz。选取台风样本中非平稳性最强的1h风速样本进行建模研究，研究得到该组1h风速发生于第21h～22h，VMD分解获取时变平均风速时，分解层数K＝2，如图10所示。
为了便于与实施例1中实测山区雷暴风速结果进行相互比较，图11～图12所示台风“杜鹃”样本的强度包络函数以及时变风谱仅示出前10min的计算结果。通过对比可以看出，台风“杜鹃”的强度包络函数变化的幅度较小，表现在其时变风谱的能量变化并不显著，可视为“弱非平稳”过程；而实施例1中所示山区雷暴风速为“强非平稳”过程。
实施例3
本实施例为对比例，采用本发明所述的建模方法和著名学者Priestly提出的基于单个观测样本的经典时变功率谱估计方法分别进行建模，可以发现Priestley方法得不到可分离形式表示的非平稳风谱(如公式(3)所示)，并不便于后续大跨度桥梁风致振动响应分析，甚至会得到无法根据随机振动理论解释的分析结果
图13～图14分别是两种建模方法的时变谱图，可以观察到：谱图在四小段时间范围内具有较强的能量，在其它时间段能量较小。因此，两种方法能够表征统一的非平稳瞬态特征，表明了本发明所述建模方法不仅具有良好捕捉瞬态特征的效果，而且可分离谱表达式符合随机振动理论中关于非平稳谱的相关定义。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，包括以下步骤：
步骤一：建立非良态风速方程
Figure FDA0002410365740000011
式中，U(t)为原始风速，
Figure FDA0002410365740000012
为确定性的时变平均风速，u(t)为均匀调制的随机脉动风速，t为时间；
所述脉动风速u(t)表示为与时间t相关的强度包络函数σu(t)与归一化平稳风速α(t)的乘积，即u(t)＝σu(t)α(t)；
所述脉动风速的功率谱函数Su(ω,t)表示为强度包络函数平方
Figure FDA0002410365740000013
与归一化平稳风速功率谱函数Sα(ω)的乘积，即
Figure FDA0002410365740000014
空间x1处的脉动风速u(x1,t)与x2处的u(x2,t)之间的互谱，表示为：
Figure FDA0002410365740000015
式中：Su(x1,x2,ω,t)为脉动风速u(x1,t)与u(x2,t)之间的互谱，Su(x1,ω,t)为x1处的脉动风速的功率谱函数，Su(x2,ω,t)为x2处的脉动风速的功率谱函数，coh(x1,x2,ω,t)为脉动风速间的相干函数，ω为圆频率；
步骤二：采用提出的判别准则并结合VMD分解法获取时变所述时变平均风速
Figure FDA0002410365740000016
和所述脉动风速u(t)；
采用ARMA-GARCH模型获取所述强度包络函数σu(t)；
采用ARMA模型或AR模型获取所述归一化平稳风速功率谱函数Sα(ω)。
2.根据权利要求1所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述步骤二中提出的判别准则为：当平均风速中最高频率fmax取值为结构基频f1的1/5～1/10时，忽略结构的动力效应，即，fmax∈[0.1f1,0.2f1]，当模态分量的中心频率fc大于0.2f1时，则认为其为脉动风成分，如下式所示：
Figure FDA0002410365740000021
3.根据权利要求2所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述步骤二中，采用VMD分解法获取时变所述时变平均风速
Figure FDA0002410365740000022
和所述脉动风速u(t)，包括以下步骤：
步骤A21：初始化VMD分解法的参数，包括保真度系数τ，第一个中心频率更新参数DC，中心频率初始化参数init，收敛准则阈值ε和惩罚因子α，其中，模态分量层数K＝1；
步骤A22:采用所述步骤A21中设定的参数对原始风速信号进行第一层分解，并获得第一层模态分量c1的中心频率；
步骤A23:将所述模态分量层数K值加1，其余参数保持不变，重复所述步骤A22，对所述原始风速信号进行K层模态分解，并获得第K层模态分量ck的中心频率；
步骤A24:判断所述第K层模态分量的中心频率ck是否大于0.2f1，所述f1为结构的基频，
如果“是”，则第1层～第K-1层模态分量的叠加视为时变平均风速，并进入步骤A25；
如果“否”，则返回步骤A23；
步骤A25:将所述原始风速信号减去所述时变平均风速
Figure FDA0002410365740000023
获得所述脉动风速u(t)。
4.根据权利要求1所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述步骤二中，采用ARMA-GARCH模型获取所述强度包络函数σu(t)时，
所述脉动风速u(t)的均值方程表示为：
φ(B)u(t)＝θ(B)ε(t)，ε(t)＝σε(t)μ(t)；
所述脉动风速u(t)的方差方程表示为：
Figure FDA0002410365740000031
式中，φ(B)为p阶AR(p)多项式，表示为φ(B)＝1+a1B+a2B2+…apBp，a＝[a1,a2,…,ap]T为AR模型的参数向量，p为AR模型的阶数，θ(B)为q阶MA(q)多项式，表示为θ(B)＝1+b1B+b2B2+…bqBq，b＝[b1,b2,…,bq]T为MA模型的参数向量，q为MA模型的阶数，B为向后移位算子，B[u(t)]＝u(t-1)，
u(t)为脉动风速，ε(t)为零均值的异方差序列，σε(t)为ε(t)的时变标准差，μ(t)为均值为0、方差为1的独立同分布随机变量序列，ηi、λj、γ为GARCH模型的参数，i＝1,2,…,m，j＝1,2,…,l，m和l为GARCH模型的阶数。
5.根据权利要求4所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述强度包络函数σu(t)的解析式表示为：
Figure FDA0002410365740000032
式中，[G1，G2，…]为格林函数。
6.根据权利要求5所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，采用ARMA-GARCH模型获取所述强度包络函数σu(t)，包括以下步骤：
步骤B21：设置ARMA模型参数p、q的取值范围，其中，p＝1,2,3,…,pmax，q＝0,1,2,…,qmax；
步骤B22：ARMA模型定阶，包括：
步骤B221：遍历ARMA模型参数p、q的取值范围，形成共计pmax·(qmax+1)组ARMA模型参数对；
步骤B222：针对每一组p、q的取值，利用Ljung-Box Q测试方法计算残差序列εpq的自相关性，并计算相应的统计值Qtest；
步骤B223：选择Qtest最小的一组p、q组合作为ARMA模型的最优参数估计；
步骤B23：GARCH模型定阶，包括通过最优原则比较，基于t分布的GARCH(1,1)模型效果最优，采用最大似然估计方法确定GARCH项的参数值；
步骤B24：根据所述的方差方程计算残差序列εpq的时变标准差；
步骤B25：利用所述残差序列εpq的时变标准差，推算脉动风速的所述强度包络函数σu(t)。
7.根据权利要求6所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述步骤B21中，pmax＝20，qmax＝pmax-1。
8.根据权利要求1所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述步骤二中，采用ARMA模型或AR模型获取所述归一化平稳风速功率谱函数Sα(ω)时，
归一化平稳风速α(t)表示为：
φ(B)α(t)＝θ(B)e(t)；
式中，φ(B)为p阶AR(p)多项式，表示为φ(B)＝1+a1B+a2B2+…apBp，a＝[a1,a2,…,ap]T为AR模型的参数向量，p为AR模型的阶数，θ(B)为q阶MA(q)多项式，表示为θ(B)＝1+b1B+b2B2+…bqBq，b＝[b1,b2,…,bq]T为MA模型的参数向量，q为MA模型的阶数，B为向后移位算子，B[u(t)]＝u(t-1)，e(t)为高斯白噪声序列且方差为常值
Figure FDA0002410365740000041
根据N个已知的观测数据α(0)，α(1)，…，α(N-1)，归一化的脉动风谱表示为：
Figure FDA0002410365740000051
当采用AR模型时，进一步简化为：
Figure FDA0002410365740000052
式中，△t为采样间隔。
9.根据权利要求8所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，所述AR模型阶数采用线性代数法或信息量准则进行确定，所述线性代数法包括奇异值分解法、Gram-Schmidt正交法，所述信息量准则包括FPE准则、AIC准则，所述AR模型参数采用最小二乘方法、Yule-Walker方法进行求解。
10.根据权利要求9所述的一种基于VMD-ARMA-GARCH模型的非良态风速建模方法，其特征在于，采用奇异值分解法结合总体最小二乘方法进行求解时，包括以下步骤：
C21：用样本的相关矩阵R代替增广矩阵B，计算B的SVD，并储存奇异值σ11≥σ11…≥σhh≥0和矩阵V；
C22：采用归一化比值法，确定增广矩阵B的有效秩p，得到ARMA模型的AR模型阶数估计；
C23：计算矩阵S(p)，
Figure FDA0002410365740000053
Figure FDA0002410365740000054
为矩阵V第j列的加窗段，上标H表示矩阵共轭转置；
C24：计算S(p)的逆矩阵S-(p)，并由ai＝S-(p)(i+1,1)/S-(P)(1,1),i＝1,…,p计算出待求AR模型参数的总体最小二乘估计值。
CN202010174623.6A 2020-03-13 2020-03-13 一种基于vmd-arma-garch模型的非良态风速建模方法 Active CN111310109B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010174623.6A CN111310109B (zh) 2020-03-13 2020-03-13 一种基于vmd-arma-garch模型的非良态风速建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010174623.6A CN111310109B (zh) 2020-03-13 2020-03-13 一种基于vmd-arma-garch模型的非良态风速建模方法

Publications (2)

Publication Number Publication Date
CN111310109A true CN111310109A (zh) 2020-06-19
CN111310109B CN111310109B (zh) 2023-03-21

Family

ID=71160602

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010174623.6A Active CN111310109B (zh) 2020-03-13 2020-03-13 一种基于vmd-arma-garch模型的非良态风速建模方法

Country Status (1)

Country Link
CN (1) CN111310109B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111967203A (zh) * 2020-08-10 2020-11-20 哈尔滨工业大学（深圳） 一种半解析半数值的大气边界层三维台风风场建模方法
CN113688774A (zh) * 2021-09-03 2021-11-23 重庆大学 基于深度学习的高层建筑风致响应预测、训练方法及装置
CN114417750A (zh) * 2022-01-20 2022-04-29 重庆大学 基于主动-被动混合试验技术的三维气动导纳识别方法及系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070168155A1 (en) * 2006-01-13 2007-07-19 Sai Ravela Statistical-deterministic approach to natural disaster prediction
WO2012105973A1 (en) * 2011-02-02 2012-08-09 Michigan Aerospace Corporation Atmospheric measurement system and method
CN103034757A (zh) * 2012-12-02 2013-04-10 中国科学院电工研究所 基于经验模态分解的风电场时频域建模方法
CN104077478A (zh) * 2014-06-26 2014-10-01 华东交通大学 一种下击暴流非平稳脉动风速的数值模拟方法
CN104992008A (zh) * 2015-06-24 2015-10-21 上海大学 基于Hilbert空间多核函数相乘的风速预测方法
US20150302313A1 (en) * 2014-04-22 2015-10-22 State Grid Corporation Of China Method of predicating ultra-short-term wind power based on self-learning composite data source
CN106611243A (zh) * 2016-12-02 2017-05-03 华北电力大学（保定） 一种基于garch模型的风速预测残差修正方法
CN206487208U (zh) * 2017-02-22 2017-09-12 中铁二院工程集团有限责任公司 一种固定拉线式山区测风塔

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070168155A1 (en) * 2006-01-13 2007-07-19 Sai Ravela Statistical-deterministic approach to natural disaster prediction
WO2012105973A1 (en) * 2011-02-02 2012-08-09 Michigan Aerospace Corporation Atmospheric measurement system and method
CN103034757A (zh) * 2012-12-02 2013-04-10 中国科学院电工研究所 基于经验模态分解的风电场时频域建模方法
US20150302313A1 (en) * 2014-04-22 2015-10-22 State Grid Corporation Of China Method of predicating ultra-short-term wind power based on self-learning composite data source
CN104077478A (zh) * 2014-06-26 2014-10-01 华东交通大学 一种下击暴流非平稳脉动风速的数值模拟方法
CN104992008A (zh) * 2015-06-24 2015-10-21 上海大学 基于Hilbert空间多核函数相乘的风速预测方法
CN106611243A (zh) * 2016-12-02 2017-05-03 华北电力大学（保定） 一种基于garch模型的风速预测残差修正方法
CN206487208U (zh) * 2017-02-22 2017-09-12 中铁二院工程集团有限责任公司 一种固定拉线式山区测风塔

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111967203A (zh) * 2020-08-10 2020-11-20 哈尔滨工业大学（深圳） 一种半解析半数值的大气边界层三维台风风场建模方法
CN111967203B (zh) * 2020-08-10 2024-05-24 哈尔滨工业大学（深圳） 一种半解析半数值的大气边界层三维台风风场建模方法
CN113688774A (zh) * 2021-09-03 2021-11-23 重庆大学 基于深度学习的高层建筑风致响应预测、训练方法及装置
CN113688774B (zh) * 2021-09-03 2023-09-26 重庆大学 基于深度学习的高层建筑风致响应预测、训练方法及装置
CN114417750A (zh) * 2022-01-20 2022-04-29 重庆大学 基于主动-被动混合试验技术的三维气动导纳识别方法及系统
CN114417750B (zh) * 2022-01-20 2024-05-24 重庆大学 基于主动-被动混合试验技术的三维气动导纳识别方法及系统

Also Published As

Publication number Publication date
CN111310109B (zh) 2023-03-21

Similar Documents

Publication Publication Date Title
CN111310109B (zh) 一种基于vmd-arma-garch模型的非良态风速建模方法
Wikle et al. A dimension-reduced approach to space-time Kalman filtering
CN107092744B (zh) 基于emd-svr的地表沉降量预测方法
CN107729592B (zh) 基于广义子空间溯踪的时变结构模态参数辨识方法
Smith et al. Adaptive correction of deterministic models to produce probabilistic forecasts
CN107679456B (zh) 一种基于极值-留数分解的海洋平台振动响应消噪方法
Tarinejad et al. Extended FDD-WT method based on correcting the errors due to non-synchronous sensing of sensors
Huang et al. Multi-taper S-transform method for evolutionary spectrum estimation
CN111880159A (zh) 一种基于lstm的雷达序列信号检测方法及系统
Yu et al. Study on the regional prediction model of PM2. 5 concentrations based on multi-source observations
Zhao et al. Nonparametric and parametric methods of spectral analysis
Huang et al. Estimating concurrent climate extremes: A conditional approach
Kang et al. Filtering Partially Observed Multiscale Systems with Heterogeneous Multiscale Methods–Based Reduced Climate Models
CN113609700A (zh) 一种山区桥梁完全非平稳风场的模拟方法
Terejanu et al. Unscented Kalman filter/smoother for a CBRN puff-based dispersion model
Farrelly et al. Determination of uncertainty in environmental noise measurements by bootstrap method
CN112100711A (zh) 一种基于arima和pso-elm的混凝土坝变形组合预报模型构建方法
Zavala et al. Dynamic harmonic regression approach to wind power generation forecasting
Stoynov Structural spectral analysis of electrochemical impedance
CN109117698B (zh) 一种基于最小均方误差准则的噪声背景估计方法
Xianmin A new method with high confidence for validation of computer simulation models of flight systems
Sušelj et al. North Sea near-surface wind climate and its relation to the large-scale circulation patterns
Cui et al. Multipoint vibration response prediction under uncorrelated multiple sources load based on elastic-net regularization in frequency domain
Manoj et al. Reduced-rank sigma-point Kalman filter and its application in ENSO model
CN113823317B (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