基于随机传递函数法的泄洪诱发场地振动预测方法与流程

本发明涉及一种基于随机传递函数法的泄洪诱发场地振动预测方法，属于水利水电工程技术领域。

背景技术：

高坝泄洪引起的场地振动反复发生且持续时间长，容易引起地基液化、建筑物基础不均匀沉降、建筑物墙体开裂等危害。对精密仪器的正常工作产生不良影响，引起建筑物产生二次噪声，还容易诱发建筑物产生共振。同时振动会干扰居民的正常生活，并对身心造成一定程度的影响。因此，在水电站实际运行中，对泄洪诱发场地振动做出预测，预判场地振动强度，评估振动带来的影响，是避免或减小上述危害的有效手段。

基于实测数据计算传递函数的场地振动预测方法，由于具有能够预测振动强度和频谱特性、可进行快速在线预测等优点，成为预测场地振动的最主要手段。但是大坝泄流诱发场地振动“动水荷载-坝体-消力池-地基-场地”相互作用和影响下的多因素耦合动力问题，对场地振动进行预测涉及到：1)动水荷载-泄流结构相互作用系统；2)泄流结构-大坝基础相互作用系统；3)地层振动波传播系统，三个子系统相互作用、相互藕合，且，各子系统的振动传递过程中存在各种干扰因素与不确定因素的影响，使得高坝泄洪诱发场地振动问题异常复杂。

以往通过传递函数对场地振动等问题进行预测，多通过建立单一激励与单一响应关系实现，对类似高坝泄洪诱发场地振动这种多激励源激励下的单一响应输出预测研究较少，目前尚没有一种较为合理的泄洪诱发场地振动预测方法，能够得到较为准确的预测结果。

技术实现要素：

鉴于上述原因，本发明的目的在于提供一种基于随机传递函数法的泄洪诱发场地振动预测方法，对每一振源的实测信号进行滤波处理后，组成多振源激励信号作为输入信号，对输出测点的实测信号进行滤波处理后，作为输出响应信号，建立传递函数，基于传递函数的无偏估计计算预测结果信号，并在预测结果信号基础上进行噪声修正，作为最终的预测结果，预测结果较为准确，能够为泄洪诱发场地振动做出预判与评估。

为实现上述目的，本发明采用以下技术方案：

一种基于随机传递函数法的泄洪诱发场地振动预测方法，包括：

对每个单一振源的实测信号进行滤波处理后，组成多振源激励源信号，作为输入信号，

对输出测点的实测信号进行滤波处理后，作为输出响应信号，

基于该输入信号与输出响应信号，建立传递函数，基于传递函数的无偏估计法，计算预测结果信号，

在预测结果信号的基础上叠加噪声序列，作为最终的预测结果。

所述对实测信号进行滤波处理的方法为：

s1：构造与实测信号的信号长度相对应的白噪声信号，对其进行eemd分解，计算噪声指标值ηj；

ηj＝σj/σ1(3)

其中，σj为第j个imf分量的标准差，σ1为第1个imf分量的标准差；

其中，cj(k)为第j个imf分量，为cj(k)的均值，n为信号长度；

s2：对实测信号进行eemd分解，计算其各imf分量的标准差σj’；

s3：计算经eemd分解的实测信号的各imf分量的噪声标准差λj；

噪声标准差λj的计算公式为：

λj＝ηjσ1`(4)

其中，j＝2，3，4，...，n，n为整数；

s4：判别实测信号各imf分量中包含的噪声成分，当λj等于或大于相应的σj’时，将该第j个imf分量直接滤除；当λj小于相应的σj’时，对该第j个imf分量进行小波阈值滤波。

所述小波阈值滤波的计算公式为：

所述噪声序列的方差为：

其中，x为滤波处理前的输出响应信号，y为滤波处理后的输出响应信号，为滤波处理后时间序列的方差，为滤波处理后的噪声序列的方差。

本发明的优点是：

1、对每一振源的实测信号及输出测点的实测信号进行改进的eemd与小波阈值联合的滤波降噪处理，既可有效滤除白噪声，又能准确保留振动信号中的有用成分，提高场地振动预测的准确性；

2、每一振源的实测信号经滤波处理后组成多振源激励信号，作为输入信号，输出测点的实测信号经滤波处理后作为输出响应信号，建立传递函数，基于传递函数的无偏估计法进行场地振动预测，能够准确地反映场地振动传播过程中的能量变化特性，通过实验验证，预测得到的场地振动频谱特性与实测结果一致或接近；

3、在传递函数输出的预测结果信号基础上，进行噪声修正，能够准确反映场地振动传播过程中的振动强度变化特性；

4、本发明的方法适用于多激励源水流脉动荷载联合激励输入系统的场地振动预测，预测结果较为准确，能够为泄洪诱发场地振动的预判与评估提供科学依据。

附图说明

图1是本发明的方法流程示意图。

图2是以多振源激励信号为输入信号计算传递函数，输出的信号振幅示意图，输入信号与输出响应信号未经滤波处理。

图3是以多振源激励信号为输入信号，计算相干传递函数的计算结果示意图，输入信号与输出响应信号未经滤波处理。

图4是以多振源激励信号为输入信号计算传递函数，进行实际工况的振动预测与实测结果的时程结果对比示意图。

图5a是以多振源激励信号为输入信号计算传递函数，进行实际工况的振动预测的傅里叶频谱计算结果示意图；

图5b是图5a所示实际工况的实测结果的傅里叶频谱计算结果示意图。

图6是本发明构造的能量相同、信号长度不同的白噪声组，比值ηj与n值之间的关系示意图。

图7是本发明构造的能量不同、信号长度相同的白噪声组，比值ηj与i值之间的关系示意图。

图8是利用本发明的方法进行实际工况预测，预测信号的振幅示意图，预测结果信号未添加噪声序列。

图9是利用本发明的方法进行实际工况预测，预测信号的频率示意图，预测结果信号未添加噪声序列。

图10是利用本发明的方法进行实际工况预测，预测结果与实测结果的时程结果对比示意图。

图11a是利用本发明的方法进行实际工况预测，预测结果的傅里叶频谱计算结果示意图。

图11b是图11a所示实际工况下，实测结果的傅里叶频谱计算结果示意图。

图12、13是利用本发明的方法，对向家坝水电站t9测点场地振动进行预测的结果示意图。

具体实施方式

以下结合附图和实施例对本发明作进一步详细的描述。

如图1所示，本发明公开的基于随机传递函数法的泄洪诱发场地振动预测方法，包括：利用改进的eemd和小波阈值滤波方法，对泄洪诱发场地的每个单一振源的实测信号进行滤波处理后，组成多振源激励信号，作为输入信号；对场地附近任意测点的实测信号进行滤波处理后，作为输出响应信号，基于该输入信号与输出响应信号建立传递函数，基于传递函数的无偏估计法进行场地振动预测，在传递函数的输出信号基础上添加噪声序列以进行噪声修正，作为最终的预测结果。具体的说：

一、以泄洪场地的多振源激励信号为输入信号，场地附近任意测点的实测信号(振动信号)为输出响应信号，建立传递函数，基于传递函数的无偏估计法进行场地振动预测。

将整个坝区作为传递函数系统的输入端，将所有激励源的垂向加速度振动信号联合作为传递函数的输入信号，选取坝区附近任意位置作为输出测点。由于振动信号的采集同时发生，因此采用时域叠加方法，对所有振源的振动信号进行叠加作为多振源激励信号。于一具体实施例中，选取12次相同泄洪工况下的多振源激励信号作为输入信号，场地特定测点的垂向加速度振动信号作为输出响应信号，根据传递函数的无偏估计方法hn，得到传递函数hn1。

如图2所示，从振动传递过程中的振幅变化过程看，传递函数hn1在0～10.0hz频段幅值较大，在5.0hz左右处出现峰值，在18.5hz、22.0hz以及29.5hz处也存在较大峰值，说明场地振动在这些频率处有放大作用；传递函数hn1在0～40.0hz一直出现小的能量峰值，与多振源激励信号联合输入受到的干扰较大有关。如图3所示，从相干传递函数计算结果看，在0～15.0hz频段，传递函数hn1的相干系数均在0.1～0.6之间，在20.0～40.0hz频段，相干系数整体有所减小，但在一些频率处出现较大值。可以看出，信号受噪声影响较大，因此，难以直接得出某一频段泄流激励荷载传递到场地地面的振动加速度大小程度。

选取某一泄洪工况，对特定测点的垂向振动情况进行预测。如图4所示，从时程结果对比图可以看出，预测结果振幅明显大于原型观测结果，如图5a、5b所示，从频谱计算结果对比图可以看出，预测的振动频率与实测结果相比，在低频、高频处，多了很多能量峰值，实测结果信号仅有0～6.0hz一处优势频带，在8.0hz左右处出现峰值，而预测结果信号有0～6.0hz和8.0～12.0hz两个优势频带。

由以上实验结果数据分析可知，从传递函数的频谱与振幅分布来看，多振源激励信号联合作为传递函数的输入信号时，各单一振源的噪声干扰较大，且各传递系统间的相互影响产生了较多的干扰，非振源引起的振动被放大，预测结果与实测结果相差较大，此种情况下的预测结果并不理想。由于泄洪激励属宽频带随机激励，且随不同泄洪工况的变化而变化，即振动传递的传递函数是随机的，分析随机传递函数的最大难题之一就是噪声问题，因而，若能够解决随机传递函数的噪声问题，有效的滤除噪声等影响因素，将有效提高场地振动的预测结果。

二、利用改进的eemd和小波阈值滤波方法，对每一单一振源的实测信号进行滤波处理后，组成降噪滤波处理后的多振源激励信号，作为输入信号；对输出测点的实测信号进行滤波处理后，作为输出响应信号，基于该输入信号与输出响应信号建立传递函数。

实测水工结构振动响应中往往混有低频干扰和白噪声，通过对实测信号进行eemd分解，低频噪声一般存在于后几阶imf分量中，滤除较为容易。而白噪声将随有用信号一起分解，存在于前几阶imf分量中，为保证信号的完整性，需对前几阶imf分量进行小波阈值滤波，尽量保留有用信号成分。

由于实测信号中混入的噪声是未知的，噪声的标准方差只能是一个估计值，严重影响了小波阈值滤波的精度和可信度。因此，可以利用白噪声eemd分解特性，确定前几阶含噪imf分量中的噪声标准差。

定义σj为第j个imf分量的标准差，则：

式中，cj(k)为信号经eemd分解后第j个imf分量，为cj(k)的均值，n为信号采集的长度。

利用正态分布的随机矩阵randn(m，n)构造白噪声。定义：

si＝2ⁱ×randn(n，1)(2)

其中，i＝0，1，2，3，...，m，2ⁱ表示白噪声强度，n为信号采集的长度，进行水工结构或场地振动原型观测时，n通常取1000～60000。根据式(2)，不同的n值与i值相组合，可构造出不同长度、不同能量的白噪声信号组。

1)构造能量相同，信号长度不同的白噪声组。

以i＝1为例，令i＝1，n＝1000，2000，...，10000，20000，...，60000，按照式(2)构造白噪声信号组。

利用eemd分解其中的每个白噪声信号，计算比值ηj：

ηj＝σj/σ1(3)

其中，σj为第j个imf分量的标准差，σ1为第1个imf分量的标准差。如图6所示，当n＝1000～10000时，比值ηj变化幅度较大，当n＞10000时，比值ηj变化幅度减小。

2)构造能量不同，信号长度相同的白噪声组。

以n＝20000为例，令i＝0～9，n＝20000，按照式(2)构造白噪声信号组。

利用eemd分解其中的每个白噪声信号，利用式(3)计算比值ηj。如图7所示，当信号长度一定时，能量不同的白噪声信号的比值ηj基本恒定。

可见，确定白噪声的信号长度后，不同能量的白噪声信号经eemd分解后的比值ηj基本稳定，或在很小范围内波动。

基于上述实验计算与分析结果，对采集的实测信号(包括每个单一振源的实测信号及特定测点输出的实测信号)进行滤波处理的过程包括如下步骤：

s1：构造与实测信号的信号长度相对应的白噪声信号，对构造的白噪声信号进行eemd分解，通过式(3)计算比值ηj，作为评价imf分量中的噪声水平的噪声指标；

s2：对实测信号进行eemd分解，计算其各imf分量的标准差σj’；

s3：计算经eemd分解的实测信号的各imf分量的噪声标准差λj；

一般情况下，包含白噪声的实测信号经eemd分解后，第一个imf分量全部为白噪声分量，因此，各imf分量的噪声标准差λj可确定为：

λj＝ηjσ1`(4)

其中，j＝2，3，4，...，n，n为整数。

s4：根据实测信号各imf分量的标准差σj’与实测信号各imf分量的噪声标准差λj的关系，判别各imf分量中包含的噪声成分，并进行相应的滤波处理。

具体为，当λj等于或大于相应的σj’时，说明该第j个imf分量全部为噪声分量，可以直接滤除；当λj小于相应的σj’时，说明该第j个imf分量中含有有用信号，应对该第j个imf分量进行小波阈值滤波，小波阈值滤波的原计算公式为：

其中，σ是噪声的标准差，这里取值即为λj，n为信号长度，得到：

利用上述滤波处理方法对采集的每个单一振源的实测信号进行滤波处理后，组成多振源激励信号作为输入信号，利用上述滤波处理方法对采集的输出测点的实测信号进行滤波处理，作为输出响应信号，基于该输入信号与输出响应信号，建立传递函数，基于传递函数的无偏估计法，进行场地振动的预测。由于输入信号与输出响应信号均通过滤波处理去除了噪声及环境影响因素，使得建立的传递函数更加准确，预测结果在振动方式(由频谱角度体现)方面更为准确。

三、基于上述第二条建立传递函数，基于传递函数的无偏估计法，计算得到预测结果信号，在此预测结果信号的基础上叠加噪声序列，作为最终的预测结果。

前述建立的传递函数，其输入信号与输出响应信号均经过滤波处理，滤除了噪声及环境因素影响，预测结果的频谱特性更接近于原型观测的频谱特性，而预测结果的振动强度要小于原型观测的振动强度。因此，可根据信噪比公式构造噪声序列，添加到传递函数输出的预测结果信号中，降低输入、输出端滤波处理对振动幅值预测精度的影响。滤波处理后的输出响应信号的信噪比为：

其中，x为滤波处理前的输出响应信号，y为滤波处理后的输出响应信号，为滤波处理后时间序列的方差，为滤波处理后的噪声序列的方差。因此，噪声序列方差可表示为：

将构造的噪声序列，添加到传递函数输出的预测结果信号中，作为最终的预测结果，使得预测结果于振动方式、振动强度方面均更接近于原型观测结果，预测结果较为准确。

四、方法的验证

对12个相同工况，以孔口、跌坎、导墙、消力底板、尾坎处为输入测点，以特定测点t9为输出测点，分别采集输入测点与输出测点的振动信号，应用前述改进的eemd和小波阈值滤波方法分别进行滤波处理，对各输入测点的实测信号进行滤波处理后，叠加组成多振源激励源，作为输入信号，对输出测点的实测信号进行滤波处理后作为输出响应信号，建立传递函数hn3，基于传递函数的无偏估计法，进行场地振动。

如图8、9所示，传递函数hn3的相干系数在0～8.0hz最高，在3.0hz左右达到峰值，说明能量在该频段传递过程中损失较少，相干系数在10.0hz以后逐渐减小。从振动传递过程中的振幅变化过程看，传递函数hn3在2.8hz处出现明显峰值，主要传递能量集中在1.0～4.0hz区间，在高频处也有少量峰值分布。可见，通过滤波处理后的输入信号与输出响应信号建立的传递函数，振动预测结果更符合实测结果。

在上述传递函数输出的预测结果信号基础上，利用式(8)计算得到噪声序列标准差，对特定测点t9的预测结果信号添加噪声序列，得到最终的时程预测结果(图10所示)及频谱预测结果(图11a、11b所示)，如图所示，依本发明的方法，能够准确的预测出场地振动的主频，同时，在0～10.0hz频段内，预测结果与实测结果的场地振动频谱分布十分相似，而且，添加了噪声序列后的振幅预测结果与实际观测结果吻合较好。

利用本发明的方法，对2013～2015年向家坝水电站t9测点场地振动情况进行预测，如图12、13所示，t9测点的预测结果与实测结果吻合较好：首先，振动强度整体上随消力池过流量增大而增大；其次，振动强度对泄洪方式比较敏感，同一泄流工况不同泄洪方式振动强度有所区别。

以上所述是本发明的较佳实施例及其所运用的技术原理，对于本领域的技术人员来说，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案基础上的等效变换、简单替换等显而易见的改变，均属于本发明保护范围之内。