CN106156404A - 基于emd的非平稳信号作用下瞬态与稳态反应计算方法 - Google Patents

Abstract

本发明提供了一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，包括：1、利用EMD将非平稳信号分解成多个IMF分量；2、对每个IMF分量中每个半波进行拟合简谐；具体：21、将当前IMF分量进行线性插值；将插值后的当前IMF分量每个半波依次拟合为简谐波；每拟合一个半波，将该半波作用在体系中，计算其瞬态反应及稳态反应；22、将当前IMF分量中每个半波求得的瞬态反应及稳态反应存储至数组中；23、对下一个IMF分量进行处理，回到步骤21，直至将每个IMF分量中的每个半波的瞬态反应及稳态反应计算出来，依次存储至数组中；3、分别将每个IMF分量对应的数组进行叠加。本发明计算用时短，结果可靠。

Description

基于EMD的非平稳信号作用下瞬态与稳态反应计算方法

技术领域

本发明涉及一种信号处理技术领域，尤其涉及一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法。

背景技术

针对地震动等非平稳信号作用下系统无法区分瞬态反应和稳态反应的技术难题，提出一种数值计算与理论推导相结合的数值解析方法，较精确、高效地从系统总响应中提取瞬态与稳态成份，为进一步分析各成份对系统的影响奠定基础。在此过程中，可突破直接分析非线性、非平稳信号难以获取工程振动规律的研究瓶颈。

系统在任意动力荷载作用下所产生的响应均可划分为瞬态反应和稳态反应两部分。传统结构动力学理论认为，由于阻尼的影响，由自由振动和伴生自由振动构成的瞬态自由振动随时间很快就消失了，因此通常不予讨论，而只强调稳态振动的作用。关于这一点，现有技术中以简谐荷载为输入，经位移解析解分析认为，当激震频率相对结构自振频率较大时，瞬态自由振动可以引起结构很大的位移反应。这一由简谐振荷载所推导分析的理论成果对于研究工程地震响应分析与隔震结构设计有着明确的工程意义：近场地震动高频加速度产生的速度脉冲将会引起长周期结构产生很大的位移，并出现于瞬态反应阶段，这种潜在破坏需在设计中加以重视；在隔震技术中，荷载频率远大于结构的自振频率，瞬态自由振动在总地震反应中的作用不能被忽略。

简谐荷载作用是结构动力学的一个经典课题，不仅因为这种荷载在工程系统中将遇到(如不平衡转动机器所产生的力)，而且因为简谐荷载下的理论结果可为体系在地震动等其他复杂荷载下的响应研究奠定基础。但是，由于地震动等复杂非平稳信号具有宽频带、非平稳特性，因此无法推导其结构响应的解析解，而只能采用数值计算的方法(如Duhamel积分，wilson-θ，Newmark-β)进行系统响应计算。又由于数值方法只能计算出系统的总体响应，无法区分其中的瞬态与稳态成份。因此要单独提取瞬态响应或稳态响应并对其进行进一步分析在实现技术上具有一定的难度。

为解决非平稳信号作用下系统响应无法区分瞬态与稳态响应的技术难题，现有技术中还提出利用小波包分解的方法来化解这一难题：非平稳信号经足够多层数的分解后，获得若干带宽较窄的小波包分量，并近似认为具有单一频率，进而利用线性调幅的方法将各小波包分量拟合为简谐波，并以简谐地震反应为理论基础推导简谐荷载作用下地震瞬态反应与稳态反应，最后经幅值还原获得非平稳荷载作用下的瞬态与稳态成分。通过线性调幅的方式将信号的小波包分量拟合为简谐波，实现了地震动作用下结构瞬态反应与稳态反应的剥离。但也存在如下缺陷：1)小波包分解需进行足够多层数的分解后，才能得到可近似认为具有单一频率的小波包分量，小波包分量的数量达到一千个之多，需要的计算机内存容量大且计算速度慢；2)小波包分解技术需选择小波基，而小波基的形式灵活多样，不同的小波基具有不同的性质，对信号的分析能力也不同，小波变换的结果很大程度上受所选小波基的影响；3)小波基的有限长会造成信号的能量泄漏，因而较难以对信号作精确的时频分析。

发明内容

本发明要解决的技术问题，在于提供一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，采用EMD分解技术将非平稳信号分解并拟合成简谐信号，地震动等非平稳信号作用下系统区分瞬态反应和稳态反应；同时克服上述小波包分解技术的小波包分量多、计算速度慢、计算精度弱等缺点。

本发明是这样实现的：

一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，包括如下步骤：

步骤1、利用EMD将非平稳信号进行分解，产生复数个具有不同特征尺度的IMF分量，每个所述IMF分量均包括复数个半波；每个所述IMF分量满足两个条件：所述IMF分量的极值点的数目与过零点的数目相等或至多相差一个；所述IMF分量的局部极大值和局部极小值定义的包络均值为零；

步骤2、从第一个所述IMF分量开始，对每个所述IMF分量中每一个半波进行拟合简谐；具体是：

步骤21、将当前IMF分量进行线性插值；将插值后的当前所述IMF分量中的每一个半波依次拟合为简谐波；每拟合一个半波，就将拟合后的该半波简谐荷载作用在体系中，并计算该半波的瞬态反应、稳态反应及全解；

步骤22、将当前所述IMF分量中每一个简谐半波求得的瞬态反应、稳态反应及全解依次存储至当前所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤23、对下一个所述IMF分量进行处理，回到步骤21，直至将每个所述IMF分量中的每一个半波的瞬态反应、稳态反应及全解都计算出来，依次存储至每个所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤3、分别将每个所述IMF分量对应的瞬态反应数组进行叠加、每个所述IMF分量对应的稳态反应数组进行叠加及每个所述IMF分量对应的全解数组进行叠加，计算出体系中非平稳信号的瞬态反应、稳态反应及全解。

进一步地，所述步骤1中利用EMD将非平稳信号进行分解，具体如下：

步骤11、非平稳信号为一原数据序列X(t)，确定该原数据序列X(t)所有的局部极大值点及局部极小值，用三次样条曲线将所有的局部极大值点连接起来，形成上包络线X_max(t)，同时将所有的局部极小值点连接起来，形成下包络线X_min(t)；

步骤12、求出上包络线及下包络线的平均值，记为包络均值m₁₁(t)，将所述原数据序列X(t)去掉该包络均值m₁₁(t)后得到新数据序列h₁₁(t)：

m₁₁(t)＝[X_max(t)+X_min(t)]/2 (1)；

h₁₁(t)＝X(t)-m₁₁(t) (2)；

步骤13、判断h₁₁(t)是否满足IMF分量的两个条件，如果不满足，则将h₁₁(t)作为原数据序列，根据式(1)与式(2)来重复上述处理过程，直到新数据序列：

h_1k(t)＝h_1(k-1)(t)-m_1k(t) (3)；

满足IMF分量的两个条件，从而得到第一个IMF分量c₁(t)：

c₁(t)＝h_1k(t) (4)；

步骤14、从X(t)中分离出c₁(t)，得到剩余序列r₁(t)：

r₁(t)＝X(t)-c₁(t) (5)；

步骤15、将r₁(t)作为一个新的原数据序列，按照以上步骤11至步骤14，依次提取第二、三……直至第n个IMF分量c_n(t)；当残量r_n(t)成为一个单调函数或小于一预定值时，分解结束。

进一步地，所述步骤21中线性插值与半波拟合为简谐波的过程，具体如下：

步骤211、将所述IMF分量进行线性插值，设线性插值后的IMF分量为C_i(t)：

t∈[t(1),t(2),…,t(m)]＝[t₁,t₂,…,t_m] (6)；

C_i(t)∈[x(t₁),x(t₂),…,x(t_m)]＝[x₁,x₂,…,x_m] (7)；

其中，t表示时刻，m表示采样总数，C_i(t)表示插值后的IMF分量；

C_i(t)的零点个数记为k；t_I(j)表示第j个极值点时刻；t_L(j)表示第j个零点时刻；半波的幅值记为U；

步骤212、将插值后的IMF分量中的半波拟合为简谐波；

第j个半波的幅值P_j分为两种情况，当第1个极值点时刻大于第1个零点时刻时，第j个半波的幅值为U(j)；当第1个极值点时刻小于第1个零点时刻时，第j个半波的幅值为U(j+1)；利用公式则表示为：

$P_j=\left\{\begin{array}{cc}U\left(j\right)& t_I\left(1\right)>t_L\left(1\right)\\ U(j+1)& t_I\left(1\right)<t_L\left(1\right)\end{array}\right.,j=1,2,...,k-1---\left(8\right);$

第j个半波的圆频率为：

${\mathrm{\&theta;}}_j=\frac{\mathrm{\&pi;}}{t_L(j+1)-t_L\left(j\right)},j=1,2,...,k-1---\left(9\right);$

第j个半波的拟合简谐函数为：

F_j(t)＝P_jcos(θ_jt-≤φ_j) t∈[t₁，t_m] 1≤j≤k-1 (10)

其中，当j＝1时，相位差

时间区间范围为t∈[t₁,t_L(2))；当2≤j≤k-2时，相位差时间区间范围为t∈[t_L(2),t_L(k-1))；当j＝k-1时，相位差时间区间范围为t∈[t_L(k-1),t_m]。

进一步地，所述步骤21中计算半波的瞬态反应、稳态反应及全解的过程，具体如下：

步骤213、将拟合为简谐波的半波荷载作用在一个阻尼比为ξ、无阻尼体系的自振圆频率为ω的体系中，其动力方程表达为：

${\stackrel{\&CenterDot;\&CenterDot;}{y}}_j\left(t\right)+2\mathrm{\&xi;}\mathrm{\&omega;}{\stackrel{\&CenterDot;}{y}}_j\left(t\right)+{\mathrm{\&omega;}}^2{y}_j\left(t\right)=-F_j\left(t\right)=-P_{j}cos({\mathrm{\&theta;}}_{j}t-{\mathrm{\&phi;}}_j)---\left(11\right);$

式中，及y_j(t)分别为第j个半波荷载作用下的结构加速度、速度和位移时程；

步骤214、经推导，所述动力方程的解析解，由两部分组成：

第j个半波的瞬态反应y_tj(t)：

第j个半波的稳态反应y_sj(t)：

第j个半波的全解y_j(t)：

y_j(t)＝y_tj(t)+y_sj(t) (14)；

式中，y₀和ν₀分别为初始位移和初始速度，为第j个半波的稳态反应y_sj(t)所对应的幅值，为反应比荷载滞后的角度，为有阻尼体系的自振圆频率，且φ_j只限于0～180°的范围；

步骤215、根据式(12)-式(14)，得出第一个半波最后时刻的位移和速度，再将该第一个半波最后时刻的位移和速度作为第二个半波的初始位移和初始速度，依次类推，直至求出最后一个半波的瞬态反应与稳态反应。

进一步地，所述步骤23中每个所述IMF分量中的每一个半波的瞬态反应、稳态反应及全解都计算出来，依次存储至每个所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中，具体如下：

将第一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_t1中、第一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_s1中及第一个所述IMF分量中的每一个半波求得的全解放入全解数组Y₁中，依次类推，直至将最后一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_tn中、最后一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_sn中及最后一个所述IMF分量中的每一个半波求得的全解放入全解数组Y_n中；第i个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i为：

i为整数且1≤i≤n (15)；

式中：k为IMF分量的零点个数，n为非平稳信号经过EMD分解得到IMF分量的个数，y_tj及y_sj分别表示在第j个半波作用下的瞬态反应和稳态反应。

进一步地，所述步骤3中将相应的数组进行叠加，具体如下：

将各个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i分别进行叠加，计算出体系中非平稳信号的瞬态反应D_t、稳态反应D_s及全解D：

$\left\{\begin{array}{c}{D}_t\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{ti}\left(t\right)\\ D_s\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{si}\left(t\right)\\ D\left(t\right)=D_t\left(t\right)+D_s\left(t\right)\end{array}\right.---\left(16\right);$

式中：n为非平稳信号经过EMD分解得到IMF分量的个数，Y_ti及Y_si分别表示在第i个拟简谐IMF分量作用下的瞬态反应和稳态反应。

本发明的优点在于：

本发明与基于小波包变换的瞬态与稳态反应计算方法相比，EMD不需要基函数，具有自适应多分辨率，小波包变换得到的分量可能达到几千个之多，而EMD分解的分量至多只有十几个，计算工作量小；与线性调幅相比，半波长拟简谐不需要拟合的信号具有单一的频率，其适用范围更广，还不需进行调幅-还原步骤，使得操作更为简便，计算精度更高。

本发明提出一种基于EMD分解的半波长拟简谐思想，并结合简谐地震反应理论提取非平稳信号作用下系统的瞬态反应与稳态反应。首先将非线性、非平稳信号通过EMD分解出复数个IMF分量，通过半波长拟合余弦函数的方法，把IMF分量拟合成多段的简谐函数；利用半波长拟简谐的方法计算地震动作用下结构位移反应数值解析解，剥离地震动作用下瞬态反应和稳态反应。这种数值与解析相结合的方法可达到以下几个目地：1)计算结果可靠，能为瞬态反应与稳态反应规律的进一步研究提供较精确的数据；2)计算无需庞大内存空间；3)计算用时短；4)理论简单、算法易实现。

附图说明

下面参照附图结合实施例对本发明作进一步的说明。

图1为本发明一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法的执行流程图。

图2为上包络线与下包络线的示意图。

图3为其中一个IMF分量的示意图。

图4为IMF分量中的第1个极值点时刻大于第1个零点时刻时，半波幅值的示意图。

图5为IMF分量中的第1个极值点时刻小于第1个零点时刻时，半波幅值的示意图。

图6为El Centro地震波的示意图。

图7为对地震波进行EMD分解后的结果示意图。

图8为IMF5与半波长拟简谐比较的结果示意图。

图9为IMF7与半波长拟简谐比较的结果示意图。

图10a为半波长拟简谐后得到的瞬态反应的示意图。

图10b为半波长拟简谐后得到的稳态反应的示意图。

图11为半波长拟简谐后得到的位移数值解析解与Duhamel积分数值解的示意图。

图12a为El Centro地震动作用下瞬态反应的示意图。

图12b为El Centro地震动作用下稳态反应的示意图。

图12c为El Centro地震动作用下全解与Duhamel积分数值解比较后的结果示意图。

具体实施方式

为使得本发明更明显易懂，现以一优选实施例，并配合附图作详细说明如下。

本技术方案可用于地震工程、地球物理探测、结构损伤侦探、机械工程、信息科学、海洋和大气科学、经济学、生态学、医学等领域所有涉及信号处理的科研和工程应用。

术语解释：

1、经验模态分解(Empirical Mode Decomposition，简称EMD)：将信号中不同尺度的波动或趋势逐级分解开来，产生一系列具有不同特征尺度的数据系列。

2、固有模态函数(Intrinsic Mode Function，简称IMF)：由EMD分解产生一系列具有不同特征尺度的数据系列，每个序列称为一个固有模态函数。固有模态函数满足两个条件：(1)信号的极值点的数目与过零点的数目相等或至多相差一个；(2)信号的局部极大值和局部极小值定义的包络均值为零。

3、半波长拟简谐：任意非平稳信号经EMD分解获得IMF分量，将IMF分量中每半个波长(两个零点位置之间)的时间差作为半个周期，两个零点之间的极小值或极大值作为幅值，将该半波长范围内的数值信号拟合为余弦函数。

本发明的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，包括如下步骤：

步骤1、利用EMD将非平稳信号进行分解，产生复数个具有不同特征尺度的IMF分量，每个所述IMF分量均包括复数个半波；每个所述IMF分量满足两个条件：所述IMF分量的极值点的数目与过零点的数目相等或至多相差一个；所述IMF分量的局部极大值和局部极小值定义的包络均值为零；

步骤2、从第一个所述IMF分量开始，对每个所述IMF分量中每一个半波进行拟合简谐；具体是：

步骤21、将当前IMF分量进行线性插值；将插值后的当前所述IMF分量中的每一个半波依次拟合为简谐波；每拟合一个半波，就将拟合后的该半波简谐荷载作用在体系中，并计算该半波的瞬态反应、稳态反应及全解；

步骤22、将当前所述IMF分量中每一个简谐半波求得的瞬态反应、稳态反应及全解依次存储至当前所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤23、对下一个所述IMF分量进行处理，回到步骤21，直至将每个所述IMF分量中的每一个半波的瞬态反应、稳态反应及全解都计算出来，依次存储至每个所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤3、分别将每个所述IMF分量对应的瞬态反应数组进行叠加、每个所述IMF分量对应的稳态反应数组进行叠加及每个所述IMF分量对应的全解数组进行叠加，计算出体系中非平稳信号的瞬态反应、稳态反应及全解。

如图1所示，本发明的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，包括如下步骤：

步骤1、利用EMD将地震动作用下非平稳信号中不同特征尺度的波动或趋势逐级分解开来，产生复数个具有不同特征尺度的IMF分量，每个所述IMF分量均包括复数个半波长；每个所述IMF分量满足两个条件：所述IMF分量的极值点的数目与过零点的数目相等或至多相差一个；所述IMF分量的局部极大值和局部极小值定义的包络均值为零；其中，通过EMD分解出IMF分量的过程具体如下：

步骤11、非平稳信号为一原数据序列X(t)，确定该原数据序列X(t)所有的局部极大值点及局部极小值，用三次样条曲线将所有的局部极大值点连接起来，形成上包络线X_max(t)，同时将所有的局部极小值点连接起来，形成下包络线X_min(t)，如图2所示；

步骤12、求出上包络线及下包络线的平均值，记为包络均值m₁₁(t)，将所述原数据序列X(t)去掉该包络均值m₁₁(t)后得到新数据序列h₁₁(t)：

m₁₁(t)＝[X_max(t)+X_min(t)]/2 (1)；

h₁₁(t)＝X(t)-m₁₁(t) (2)；

步骤13、判断h₁₁(t)是否满足IMF分量的两个条件，如果不满足，则将h₁₁(t)作为原数据序列，根据式(1)与式(2)来重复上述处理过程，直到新数据序列：

h_1k(t)＝h_1(k-1)(t)-m_1k(t) (3)；

满足IMF分量的两个条件，从而得到第一个IMF分量c₁(t)：

c₁(t)＝h_1k(t) (4)；

步骤14、从X(t)中分离出c₁(t)，得到剩余序列r₁(t)：

r₁(t)＝X(t)-c₁(t) (5)；

步骤15、将r₁(t)作为一个新的原数据序列，按照以上步骤11至步骤14，依次提取第二、三……直至第n个IMF分量c_n(t)；当残量r_n(t)成为一个单调函数或小于一预定值时，分解结束；

步骤2、设定分离出的IMF分量的总量为n，i为IMF分量的序号，i为整数且1≤i≤k-1，j为所述IMF分量中半波的序号，j的个数是根据每个IMF分量中半波个数的不同而设定，所述IMF分量中零点的个数为k；初始化i为1；

步骤3、判断i是否小于或等于n，若i≤n，则将第i个IMF分量中每个半波长进行拟简谐，初始化j为1；若i＞n，则进入步骤7；

步骤4、判断j是否小于或等于k-1，若j≤k-1，则将插值后的第i个IMF分量中第j个半波拟合为简谐波；将拟合后的半波荷载作用在体系中，计算获得该半波的瞬态反应、稳态反应及全解，进入步骤5；若j＞k-1，则进入步骤6；

步骤5、逐一增加第i个IMF分量中半波的序号，即j＝j+1，进入步骤4；

步骤6、将第i个IMF分量中每个半波的瞬态反应、稳态反应及全解放入对应的数组中；逐一增加所述IMF分量的序号，即i＝i+1，进入步骤3；

步骤7、将各个IMF分量的瞬态反应、稳态反应及全解分别进行叠加，计算出体系中非平稳信号的瞬态反应D_t、稳态反应D_s及全解D。

其中，步骤3中的线性插值过程具体为：

将所述IMF分量进行线性插值，设线性插值后的IMF分量为C_i(t)：

t∈[t(1),t(2),…,t(m)]＝[t₁,t₂,…,t_m] (6)；

C_i(t)∈[x(t₁),x(t₂),…,x(t_m)]＝[x₁,x₂,…,x_m] (7)；

式中，t表示时刻，m表示采样总数，C_i(t)表示插值后的表示IMF分量；

C_i(t)的零点个数记为k；t_I(j)表示第j个极值点时刻；t_L(j)表示第j个零点时刻；半波的幅值记为U；

其中，步骤4中将半波拟合为简谐波的过程具体为：

如图4和图5所示，将插值后的IMF分量中的半波拟合为简谐波；

第j个半波的幅值P_j分为两种情况，当第1个极值点时刻大于第1个零点时刻时，第j个半波的幅值为U(j)；当第1个极值点时刻小于第1个零点时刻时，第j个半波的幅值为U(j+1)；利用公式则表示为：

$P_j=\left\{\begin{array}{cc}U\left(j\right)& t_I\left(1\right)>t_L\left(1\right)\\ U(j+1)& t_I\left(1\right)<t_L\left(1\right)\end{array}\right.,j=1,2,...,k-1---\left(8\right);$

第j个半波的圆频率为：

${\mathrm{\&theta;}}_j=\frac{\mathrm{\&pi;}}{t_L(j+1)-t_L\left(j\right)},j=1,2,...,k-1---\left(9\right);$

第j个半波的拟合简谐函数为：

F_j(t)＝P_jcos(θ_jt-≤φ_j) t∈[t₁，t_m] 1≤j≤k-1 (10);

其中，当j＝1时，相位差

时间区间范围为t∈[t₁,t_L(2))；当2≤j≤k-2时，相位差时间区间范围为t∈[t_L(2),t_L(k-1))；当j＝k-1时，相位差时间区间范围为t∈[t_L(k-1),t_m]。

其中，所述步骤4中将拟合后的半波荷载作用在体系中，计算获得该半波的瞬态反应、稳态反应及全解的过程，具体如下：

步骤41、将拟合为简谐波的半波荷载作用在一个阻尼比为ξ、无阻尼的自振圆频率为ω的体系中，其动力方程表达为：

${\stackrel{\&CenterDot;\&CenterDot;}{y}}_j\left(t\right)+2\mathrm{\&xi;}\mathrm{\&omega;}{\stackrel{\&CenterDot;}{y}}_j\left(t\right)+{\mathrm{\&omega;}}^2{y}_j\left(t\right)=-F_j\left(t\right)=-P_{j}cos({\mathrm{\&theta;}}_{j}t-{\mathrm{\&phi;}}_j)---\left(11\right);$

式中，及y_j(t)分别为第j个半波荷载作用下的结构加速度、速度和位移时程；

步骤42、经推导，所述动力方程的解析解，由两部分组成：

第j个半波的瞬态反应y_tj(t)：

第j个半波的稳态反应y_sj(t)：

第j个半波的全解y_j(t)：

y_j(t)＝y_tj(t)+y_sj(t) (14)；

式中，y₀和v₀分别为初始位移和初始速度，为第j个半波的稳态反应y_sj(t)所对应的幅值，为反应比荷载滞后的角度，为有阻尼体系的自振圆频率，且φ_j只限于0～１80°

的范围；

步骤43、根据式(12)-式(14)，得出第一个半波最后时刻的位移和速度，再将该第一个半波最后时刻的位移和速度作为第二个半波的初始位移和初始速度，依次类推，直至求出最后一个半波的瞬态反应与稳态反应。

其中，所述步骤6中将第i个IMF分量中每个半波的瞬态反应、稳态反应及全解放入对应的数组中，具体如下：

将第一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_t1中、第一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_s1中及第一个所述IMF分量中的每一个半波求得的全解放入全解数组Y₁中，依次类推，直至将最后一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_tn中、最后一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_sn中及最后一个所述IMF分量中的每一个半波求得的全解放入全解数组Y_n中；第i个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i为：

i为整数且1≤i≤n (15)；

式中：k为IMF分量的零点个数，n为非平稳信号经过EMD分解得到IMF分量的个数，y_tj及y_sj分别表示在第j个半波作用下的瞬态反应和稳态反应。

其中，所述步骤7中将相应的数组进行叠加，具体如下：

将各个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i分别进行叠加，计算出体系中非平稳信号的瞬态反应D_t、稳态反应D_s及全解D：

$\left\{\begin{array}{c}{D}_t\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{ti}\left(t\right)\\ D_s\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{si}\left(t\right)\\ D\left(t\right)=D_t\left(t\right)+D_s\left(t\right)\end{array}\right.---\left(16\right);$

式中：n为非平稳信号经过EMD分解得到IMF分量的个数，Y_ti及Y_si分别表示在第i个拟简谐IMF分量作用下的瞬态反应和稳态反应。

本发明的一较佳实施例为：(以结构在实际地震动作用下瞬态反应与稳态反应计算中的应用为例进行分析说明。)

如图6所示，选用El Centro南北向地震动，加速度时间间隔为0.02s，持时30s，输入到自振周期为3s，阻尼比为0.05的单自由度系统中，通过本发明中的计算方法提取该体系的瞬态反应和稳态反应。各步骤及有效性分析如下：

(1)对该地震动进行EMD分解，形成9个IMF分量(IMF1、IMF2、IMF3、IMF4、IMF5、IMF6、IMF7、IMF8及IMF9)和1个残余量(RES)，(如图7所示)；

(2)为保证高频分量中半个波长内有足够的分析数据；将IMF分量进行线性插值。如图3所示，IMF分量中每个半波长曲线与正余弦函数曲线具有一定的相似性，利用两者的这种局部相似性，将插值后的IMF分量利用公式(6)-(10)进行半波长拟简谐；

(3)比较任意两个IMF分量(IMF5高频分量和IMF7低频分量)分别与半波长拟简谐的结果，如图8和图9所示；图8和图9表明，经半波长拟简谐的结果与IMF5、IMF7分量几乎重合，说明半波长拟简谐的方法是可行的；

(4)将图8中的IMF5拟简谐分量输入到自振周期为3s、阻尼比为0.0的单自由度系统中，经公式(12)-(14)计算获得半波长拟简谐位移数值解析解。其中，瞬态反应与稳态反应如图10a与图10b所示，半波长拟简谐位移数值解析解与Duhamel积分数值解的比较如图11所示。图11显示，二者具有很高的吻合度，说明本发明提出的方法是可行的且准确的；

(5)通过公式(16)叠加所有分量拟简谐作用的结果，获得完整Elcentro地震动作用下体系的瞬态反应、稳态反应，并将总位移数值解析解与Duhamel积分数值解进行比较；如图12a、图12b与图12c所示，本发明提出的方法与Duhamel积分数值解也具有较高的吻合度，说明基于EMD拟简谐计算地震动的瞬态反应与稳态反应是合理的。

本发明的优点如下：

本发明与基于小波包变换的瞬态与稳态反应计算方法相比，EMD不需要基函数，具有自适应多分辨率，小波包变换得到的分量可能达到几千个之多，而EMD分解的分量至多只有十几个，计算工作量小；与线性调幅相比，半波长拟简谐不需要拟合的信号具有单一的频率，其适用范围更广，还不需进行调幅-还原步骤，使得操作更为简便，计算精度更高。

本发明提出一种基于EMD分解的半波长拟简谐思想，并结合简谐地震反应理论提取非平稳信号作用下系统的瞬态反应与稳态反应。首先将非线性、非平稳信号通过EMD分解出复数个IMF分量，通过半波长拟合余弦函数的方法，把IMF分量拟合成多段的简谐函数；利用半波长拟简谐的方法计算地震动作用下结构位移反应数值解析解，剥离地震动作用下瞬态反应和稳态反应。这种数值与解析相结合的方法可运用于其他非平稳信号作用下的瞬态反应与稳态反应计算，同时可达到以下几个目地：1)计算结果可靠，能为瞬态反应与稳态反应规律的进一步研究提供较精确的数据；2)计算无需庞大内存空间；3)计算用时短；4)理论简单、算法易实现。

虽然以上描述了本发明的具体实施方式，但是熟悉本技术领域的技术人员应当理解，我们所描述的具体的实施例只是说明性的，而不是用于对本发明的范围的限定，熟悉本领域的技术人员在依照本发明的精神所作的等效的修饰以及变化，都应当涵盖在本发明的权利要求所保护的范围内。

Claims (6)

1.一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：包括如下步骤：

步骤1、利用EMD将非平稳信号进行分解，产生复数个具有不同特征尺度的IMF分量，每个所述IMF分量均包括复数个半波；每个所述IMF分量满足两个条件：所述IMF分量的极值点的数目与过零点的数目相等或至多相差一个；所述IMF分量的局部极大值和局部极小值定义的包络均值为零；

步骤2、从第一个所述IMF分量开始，对每个所述IMF分量中每一个半波进行拟合简谐；具体是：

步骤21、将当前IMF分量进行线性插值；将插值后的当前所述IMF分量中的每一个半波依次拟合为简谐波；每拟合一个半波，就将拟合后的该半波简谐荷载作用在体系中，并计算该半波的瞬态反应、稳态反应及全解；

步骤22、将当前所述IMF分量中每一个拟合简谐半波求得的瞬态反应、稳态反应及全解依次存储至当前所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤23、对下一个所述IMF分量进行处理，回到步骤21，直至将每个所述IMF分量中的每一个半波的瞬态反应、稳态反应及全解都计算出来，依次存储至每个所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中；

步骤3、分别将每个所述IMF分量对应的瞬态反应数组进行叠加、每个所述IMF分量对应的稳态反应数组进行叠加及每个所述IMF分量对应的全解数组进行叠加，计算出体系中非平稳信号的瞬态反应、稳态反应及全解。

2.根据权利要求1所述的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：所述步骤1中利用EMD将非平稳信号进行分解，具体如下：

步骤11、非平稳信号为一原数据序列X(t)，确定该原数据序列X(t)所有的局部极大值点及局部极小值，用三次样条曲线将所有的局部极大值点连接起来，形成上包络线X_max(t)，同时将所有的局部极小值点连接起来，形成下包络线X_min(t)；

步骤12、求出上包络线及下包络线的平均值，记为包络均值m₁₁(t)，将所述原数据序列X(t)去掉该包络均值m₁₁(t)后得到新数据序列h₁₁(t)：

m₁₁(t)＝[X_max(t)+X_min(t)]/2 (1)；

h₁₁(t)＝X(t)-m₁₁(t) (2)；

步骤13、判断h₁₁(t)是否满足IMF分量的两个条件，如果不满足，则将h₁₁(t)作为原数据序列，根据式(1)与式(2)来重复上述处理过程，直到新数据序列：

h_1k(t)＝h_1(k-1)(t)-m_1k(t) (3)；

满足IMF分量的两个条件，从而得到第一个IMF分量c₁(t)：

c₁(t)＝h_1k(t) (4)；

步骤14、从X(t)中分离出c₁(t)，得到剩余序列r₁(t)：

r₁(t)＝X(t)-c₁(t) (5)；

步骤15、将r₁(t)作为一个新的原数据序列，按照以上步骤11至步骤14，依次提取第二、三……直至第n个IMF分量c_n(t)；当残量r_n(t)成为一个单调函数或小于一预定值时，分解结束。

3.根据权利要求2所述的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：所述步骤21中线性插值与半波拟合为简谐波的过程，具体如下：

步骤211、将所述IMF分量进行线性插值，设线性插值后的IMF分量为C_i(t)：

t∈[t(1),t(2),…,t(m)]＝[t₁,t₂,…,t_m] (6)；

C_i(t)∈[x(t₁),x(t₂),…,x(t_m)]＝[x₁,x₂,…,x_m] (7)；

其中，t表示时刻，m表示采样总数，C_i(t)表示插值后的IMF分量；

C_i(t)的零点个数记为k；t_I(j)表示第j个极值点时刻；t_L(j)表示第j个零点时刻；半波的幅值记为U；

步骤212、将插值后的IMF分量中的半波拟合为简谐波；

第j个半波的幅值P_j分为两种情况，当第1个极值点时刻大于第1个零点时刻时，第j个半波的幅值为U(j)；当第1个极值点时刻小于第1个零点时刻时，第j个半波的幅值为U(j+1)；利用公式则表示为：

$P_j=\left\{\begin{array}{cc}U\left(j\right)& t_I\left(1\right)>t_L\left(1\right)\\ U(j+1)& t_I\left(1\right)<t_L\left(1\right)\end{array}\right.,j=1,2,...,k-1---\left(8\right);$

第j个半波的圆频率为：

${\mathrm{\&theta;}}_j=\frac{\mathrm{\&pi;}}{t_L(j+1)-t_L\left(j\right)},j=1,2,...,k-1---\left(9\right);$

第j个半波的拟合简谐函数为：

F_j(t)＝P_jcos(θ_jt-φ_j)t∈[t₁，t_m]1≤j≤k-1 (10)；

其中，当j＝１时,相位差时间区间范围为t∈[t₁,t_L(2))；当2≤j≤k-2时，相位差时间区间范围为t∈[t_L(2),t_L(k-1);当j＝k-1时，相位差时间区间范围为t∈[t_L(k-1),t_m]。

4.根据权利要求3所述的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：所述步骤21中计算半波的瞬态反应、稳态反应及全解的过程，具体如下：

步骤213、将拟合为简谐波的半波荷载作用在一个阻尼比为ξ、无阻尼的自振圆频率为ω的体系中，其动力方程表达为：

${\stackrel{\&CenterDot;\&CenterDot;}{y}}_j\left(t\right)+2\mathrm{\&xi;}\mathrm{\&omega;}{\stackrel{\&CenterDot;}{y}}_j\left(t\right)+{\mathrm{\&omega;}}^2{y}_j\left(t\right)=-F_j\left(t\right)=-P_{j}cos({\mathrm{\&theta;}}_{j}t-{\mathrm{\&phi;}}_j)---\left(11\right);$

式中，及y_j(t)分别为第j个半波荷载作用下的结构加速度、速度和位移时程；

步骤214、经推导，所述动力方程的解析解，由两部分组成：

第j个半波的瞬态反应y_tj(t)：

第j个半波的稳态反应y_sj(t)：

第j个半波的全解y_j(t)：

y_j(t)＝y_tj(t)+y_sj(t) (14)；

式中，y₀和v₀分别为初始位移和初始速度，为第j个半波的稳态反应y_sj(t)所对应的幅值，为反应比荷载滞后的角度，为有阻尼体系的自振圆频率，且φ_j只限于0～180°的范围；

步骤215、根据式(12)-式(14)，得出第一个半波最后时刻的位移和速度，再将该第一个半波最后时刻的位移和速度作为第二个半波的初始位移和初始速度，依次类推，直至求出最后一个半波的瞬态反应与稳态反应。

5.根据权利要求4所述的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：所述步骤23中每个所述IMF分量中的每一个半波的瞬态反应、稳态反应及全解都计算出来，依次存储至每个所述IMF分量对应的瞬态反应数组、稳态反应数组及全解数组中，具体如下：

将第一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_t1中、第一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_s1中及第一个所述IMF分量中的每一个半波求得的全解放入全解数组Y₁中，依次类推，直至将最后一个所述IMF分量中的每一个半波求得的瞬态反应放入瞬态反应数组Y_tn中、最后一个所述IMF分量中的每一个半波求得的稳态反应放入稳态反应数组Y_sn中及最后一个所述IMF分量中的每一个半波求得的全解放入全解数组Y_n中；第i个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i为：

式中：k为IMF分量的零点个数，n为非平稳信号经过EMD分解得到IMF分量的个数，yt_j及y_sj分别表示在第j个半波作用下的瞬态反应和稳态反应。

6.根据权利要求5所述的一种基于EMD的非平稳信号作用下瞬态与稳态反应计算方法，其特征在于：所述步骤3中将相应的数组进行叠加，具体如下：

将各个IMF分量的瞬态反应Y_ti、稳态反应Y_si及全解Y_i分别进行叠加，计算出体系中非平稳信号的瞬态反应D_t、稳态反应D_s及全解D：

$\left\{\begin{array}{c}{D}_t\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{ti}\left(t\right)\\ D_s\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{Y}_{si}\left(t\right)\\ D\left(t\right)=D_t\left(t\right)+D_s\left(t\right)\end{array}\right.---\left(16\right);$

式中：n为非平稳信号经过EMD分解得到IMF分量的个数，Y_ti及Y_si分别表示在第i个拟简谐IMF分量作用下的瞬态反应和稳态反应。