龙羊峡重力拱坝模态识别与材料参数反演

侯志英，周天宇，马宁，刘昌伟，向致谦，潘坚文，王进廷

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

龙羊峡重力拱坝模态识别与材料参数反演 PDF

- ORCID：
侯志英 ¹
✉
- ORCID：
周天宇 ²
- ORCID：
马宁 ³
- ORCID：
刘昌伟 ²
- ORCID：
向致谦 ²
- ORCID：
潘坚文 ²
- ORCID：
王进廷 ²
✉

1. 国家电投集团青海黄河电力技术有限责任公司西宁,810016； 2. 清华大学水沙科学与水利水电工程国家重点实验室北京,100084； 3. 国家电投集团大坝管理中心西安,710000

中图分类号： TV642.4⁺4； TH825

最近更新：2024-09-02

DOI：10.16450/j.cnki.issn.1004-6801.2024.04.020

摘要

基于龙羊峡重力拱坝实测振动数据，利用随机子空间法进行了大坝模态参数识别，建立了坝体‑地基‑库水系统的三维流固耦合有限元模型。根据识别出的一阶频率，分别考虑均质地基及分层地基，反演了坝体与基岩材料参数，并通过有限元软件ABAQUS进行了大坝的动力特征分析。结果表明：龙羊峡重力拱坝的一阶频率为2.778~3.196 Hz，反演得到坝体的弹性模量约为50 GPa~52 GPa；对于均值地基，其弹性模量的反演值为30 GPa；对于分层地基，2 490 m高程以下区域的弹性模量为34 GPa；实测振动数据反演得到的坝体弹性模量和地基弹性模量高于其设计值；龙羊峡重力拱坝的前3阶振型分别为正对称、反对称与正对称，其1阶频率和2阶频率相差较大，与重力坝自振频率分布特征相似。

关键词

龙羊峡拱坝; 模态参数识别; 有限元; 动力特征; 流固耦合; 模态分析

引言

随着社会经济的发展，水利工程的安全性需求不断增加，同时经济型需求也不断增长^［

1］。修建的大部分大坝除了防洪发电等任务以外，其自身的抗震性能也至关重要^{［参考文献 2

百度学术}2］。坝体的模态参数包括坝体自振频率、阻尼比和模态振型等，能够为结构系统的振动特性和安全评估提供依据^{［参考文献 3

百度学术}3］，对于坝体健康运行具有重要意义。

通过对实测的振动数据进行模态参数识别是进行大坝动力学分析的重要手段^［

4］。Ellis等^{［参考文献 5

百度学术}5］对美国加利福尼亚州的某大坝进行环境振动试验，利用快速傅里叶变换和最大熵方法的频谱分析和瀑布图分析来识别单拱和全坝的模态特性。Sevim^{［参考文献 6

百度学术}6］在土耳其某拱坝开展环境振动测试，利用增强频域分解法识别坝体结构的模态参数。文献［7］对某双曲拱坝开展环境振动测试，识别出坝体结构的模态参数。基于环境激励数据， Weng等^{［参考文献 8

百度学术}8］使用随机子空间（stochastic subspace identification，简称SSI）算法识别翡翠拱坝的模态参数。李帅等^{［参考文献 9

百度学术}9］分别使用频域分解法和SSI算法识别溪洛渡拱坝的模态参数。由于SSI算法需要调整用户定义参数，参数选择的合理性会显著影响识别结果，而大体积混凝土拱坝涉及到大坝‑库水‑地基的相互作用，故识别其模态参数具有一定困难。除了实测数据处理和模型试验等方法获得模态参数之外，通过有限元模型计算也是大坝模态分析的重要手段^{［参考文献 10

百度学术}10］。张楚汉等^{［参考文献 11

百度学术}11］利用有限元子结构方法分析乐响洪甸拱坝的动水压力，考虑了不同水库模型。邱奕翔等^{［参考文献 12

百度学术}12］建立了拉西瓦拱坝的拱坝‑地基‑库水系统有限元模型，研究了库水模拟对拱坝动力学特征的影响。

基于有限元模型获取的结构振动参数，其数值通常取决于输入的材料弹性模量及密度等。对于大坝等大体积混凝土结构及其基础，弹性模量在实际工作过程中与设计值存在较大差异，且不易直接通过现场测量等方法得到。为了获得真实的坝体混凝土参数，可结合上述2种方法，利用实测的坝体模态参数，建立大坝的三维有限元模型，并进行材料参数反演。程琳等^［

13］提出了一种基于强震观测和多输出支持向量机模型（简称M‑SVM模型）的大坝材料动参数反演方法，验证了基于振动数据反演坝体参数的有效性及合理性。李火坤等^{［参考文献 14

百度学术}14］基于敏感性分析及粒子群算法进行了二滩拱坝的坝体和地基动弹性模量参数反演，其反演结果表明，坝体和地基的弹性模量与其设计值存在较大差异。

笔者利用协方差驱动的随机子空间法进行龙羊峡重力拱坝的模态参数识别。建立龙羊峡重力拱坝坝体‑地基‑库水流固耦合三维有限元模型，根据识别出的模态参数进行了大坝和地基材料参数的反演，基于反演参数分析了龙羊峡重力拱坝的动力学特征。

1 工程概况

龙羊峡水电站位于青海省共和县境内黄河干流上，水库正常蓄水位为2 600 m，相应库容为247亿m³。厂房装有4台单机容量为32万kW的水轮发电机组，总装机容量为128万kW，年发电量为60亿千瓦时。龙羊峡水电站的挡水建筑物由混凝土重力拱坝、左右岸重力墩和混凝土副坝组成。重力拱坝坝顶高程为2 610 m，最大坝高为178 m。最大底宽为80 m，主坝分18个坝段，前沿长为396 m。坝后厂房机组引水道及中、深、底孔泄水道均从坝内通过，左右岸设有重力墩。左右岸挡水副坝均为重力坝，长度分别为375和341.5 m。右岸副坝与重力墩之间设有宽为40 m的两孔开敞式溢洪道，挡水建筑物前沿总长为1 227 m。

龙羊峡拱坝地震台阵共有13个观测点（坝内12个观测点和1个自由场点），共有39个监测方向，测点沿拱冠梁从坝顶到坝基、1/4拱圈处坝肩布置。传感器测量方向为坝体的水平径向。图1为龙羊峡重力拱坝地震台阵。地震台阵自2011年开始采集加速度数据。

图1 龙羊峡重力拱坝地震台阵

Fig.1 Seismologic array of Longyangxia gravity-arch dam

2 模态参数识别

随机子空间法^［

15］是由系统随机响应协方差矩阵构造托普利兹矩阵，对其进行奇异值分解，得到系统的扩展可观和可控矩阵，再对系统状态矩阵进行特征值分解而获得模态参数的方法。

白噪声激励的n自由度系统的离散状态模型可表示为

x_{k + 1} = A_{d} x_{k} + w_{k}

（1）

y_{k} = C_{d} x_{k} + v_{k}

（2）

其中：x_k，y_k，A_d和C_d分别为离散k时刻的状态向量、实测结构响应向量、状态矩阵和输出矩阵；w_k为过程噪声；v_k为测量噪声；w_k和v_k均假定为零均值的白噪声。

利用实测数据y_k构造Hankel矩阵

H = \frac{1}{\sqrt[]{j}} [\frac{\begin{matrix} y_{0} & y_{1} & \dots & y_{j - 1} \\ y_{1} & y_{2} & \dots & y_{j} \\ ⋮ & ⋮ & ⋮ \\ y_{g - 1} & y_{g} & \dots & y_{g + j - 2} \end{matrix}}{\begin{matrix} y_{g} & y_{g + 1} & \dots & y_{g + j - 1} \\ y_{g + 1} & y_{g + 2} & \dots & y_{g + j} \\ ⋮ & ⋮ & ⋮ \\ y_{i - 1} & y_{i} & \dots & y_{i + j - 2} \end{matrix}}] = \frac{1}{\sqrt[]{j}} [\frac{Y_{p}}{Y_{f}}]

（3）

其中：Y_p和Y_f分别为过去输出矩阵和将来输出矩阵。

输出信号的协方差矩阵R_i可表示为

R_{i} = E [y_{k + i} {y_{k}}^{T}] \approx \frac{1}{j} \sum_{k = 0}^{j - 1} y_{k + i} {y_{k}}^{T} = C_{d} {A_{d}}^{i - 1} G

（4）

其中：E［•］为期望算子；G=E［x_k+iy_k^T］；T为转置算子。

根据式（4）构造一个块Toeplitz矩阵

T_{1 | n} = [\begin{matrix} R_{n} & R_{n - 1} & \dots & R_{1} \\ R_{n + 1} & R n & \dots & R_{2} \\ ⋮ & ⋮ & ⋮ \\ R_{2 n - 1} & R_{2 n - 2} & \dots & R_{n} \end{matrix}]

（5）

其中：n为模型阶数。

对式（5）进行奇异值分解，得到状态矩阵 $A_{d}$ 和输出矩阵 $C_{d}$ 。对状态矩阵 $A_{d}$ 进行特征值分解

A_{d} = Ψ Λ Ψ^{- 1}

（6）

其中： $Ψ$ 为奇异向量矩阵； $Λ$ 为包含特征值μ_i的对角矩阵。

基于特征值μ_i，计算得到复特征值为

λ_{i}, {λ_{i}}^{*} = \frac{l n (μ_{i})}{Δ t} = - 2 π ξ_{i} f_{i} \pm 2 j π f_{i} \sqrt[]{1 - {ξ_{i}}^{2}}

（7）

其中：f_i和ξ_i分别为第i阶频率和阻尼比。

结构模态参数的计算公式为

f_{i} = \sqrt[]{({λ_{i}}^{R})^{2} + ({λ_{i}}^{I})^{2}} / 2 π

（8）

ξ_{i} = - λ_{i}^{R} / 2 π f_{i}

（9）

ϕ = C_{d} Ψ

（10）

其中：ϕ为复振型向量；上标I和R分别代表虚部和实部。

在协方差驱动的随机子空间法中，用户定义参数参考文献［

16］方法。生成的稳定图中存在虚假模态，笔者采用文献［17］提出的基于SSI算法和基于密度的噪声应用空间聚类算法来消除虚假模态。该聚类算法可以自动剔除异常模态点，且其用户定义参数具有较强的鲁棒性。

选取2013—2018年的1月和7月各自某一天的10：00~11：00时的加速度数据，并选取坝顶7个测点的径向数据进行随机子空间法的模态识别。图2为2013年7月1日所有测点的实测加速度时域信号。

图2 2013年7月1日所有测点实测加速度时域信号图

Fig.2 Time domain acceleration signals at all measuring points

龙羊峡重力拱坝1阶自振频率如图3所示。由于龙羊峡重力拱坝在2018年发生了较为明显的水位上升，因此其第1阶自振频率出现了较为明显的下降。

图3 龙羊峡重力拱坝1阶自振频率

Fig.3 First order frequency of Longyangxia gravity arch dam

图4为龙羊峡重力拱坝前3阶振型。第1阶和第3阶是正对称，第2阶是反对称。

图4 龙羊峡重力拱坝前3阶振型

Fig.4 First three modes of Longyangxia gravity arch dam

3 计算模型

3.1　有限元模型

建立龙羊峡重力拱坝有限元模型，如图5所示。库水通过声学材料进行模拟，假定库水不可压缩。为综合考虑多种水位的工况，分别建立了水位为138，148及158 m的库水模型。有限元模型共考虑坝址区域8条主要断层，其中：NW组包含F18，F73，F71，F7等断层；NE组包括A2，F120，F57，G4等断层。有限元模型下游边界延伸206.15 m，上游延伸至G4断层以上，左右岸边界各延伸200 m。有限元坝体‑地基部分共包含421 622个节点、2 380 713个单元，其中：坝体部分包含38 386个单元。对于库水部分，水位为138 m模型包含69 084个单元，水位为148 m模型包含73 447个单元，水位为158 m模型包含79 309个单元。水位为148 m的库水模型如图5（b）所示。在有限元模型四周截断边界施加法向约束，模型底部施加全部约束。

图5 龙羊峡重力拱坝有限元模型

Fig.5 FEM of Longyangxia gravity‑arch dam

3.2　材料参数

有限元模型材料参数以设计值作为参数反演的初始状态，坝体混凝土密度取为2 400 kg/m³，动弹性模量取为25 GPa，泊松比为0.2。地基采用无质量地基，密度为0。大坝基岩力学参数如表1所示。

表1 大坝基岩力学参数

Tab.1 Dam bedrock material mechanical parameters

弹性模量/GPa	高程/m
弹性模量/GPa	2 610	2 600	2 590	2 580	2 570	2 550	2 520	2 490	2 463	2 443	2 436	2 432
左岸	8	8	8	8	10	16	16	20	20	20	20	20
右岸	8	8	8	8	10	16	16	20	20	20	20	20

3.3　动水压力模拟方法

坝体‑地基‑库水中，考虑库水与坝体及地基的动力相互作用，其中假设库水为不可压缩的理想流体。根据文献［

17］，以动水压力P为基本变量的流体波动方程为

\frac{\partial^{2} P}{\partial x^{2}} + \frac{\partial^{2} P}{\partial y^{2}} + \frac{\partial^{2} P}{\partial z^{2}} = \frac{1}{C^{2}} \frac{\partial^{2} P}{\partial t^{2}}

（11）

其中：C为水中声波波速。

该方程需满足4个边界条件，分别如下。

库水表面边界条件为

{\dot{u}}_{n} = \dot{P} / ρ g

（12）

其中： ${\dot{u}}_{n}$ 为库水表面法相速度； $ρ$ 为对应的库水密度。

坝体‑库水边界条件为

\partial P / \partial n = - ρ {\ddot{u}}_{d n}

（13）

其中： ${\ddot{u}}_{d n}$ 为坝体外法向加速度。

地基‑库水边界条件为

\partial P / \partial n = - ρ {\ddot{u}}_{d n} + q \dot{P}

（14）

其中：q为水库边界阻抗系数。

库水截断边界辐射条件为

\frac{\partial P}{\partial n} = \frac{c o s θ}{C} \dot{P}

（15）

其中：n为截断边界的内法线方向； $θ$ 为入射角。

在Abaqus中，库水表面边界设定为自由表面边界，坝体‑库水及地基‑库水边界通过设置接触软件Tie进行模拟，库水截断边界通过吸收边界模拟。流固耦合模型中，库水密度取为1 000 kg/m³，体积模量为2×10¹¹ GPa。

4 坝体与基岩材料参数反演

参数反演以龙羊峡拱坝一阶自振频率作为目标函数，采用试算法，选取一定范围的坝体弹性模量与地基弹性模量，计算其坝体的自振频率，选出与实测自振频率最相符的一组材料参数。地基弹性模量采用2种调整方式：①按照初始参数中地基各层弹性模量的值进行等比例放缩；②将地基弹性模量改为统一数值进行反演。记坝体弹性模量为x，2 490 m高程以下地基弹性模量为y，则上述2种地基处理方法在参数反演过程中均只需考虑参数x及y的变化。记实测的坝体自阵频率为f_h，有限元模型计算得到的频率为F_h（x，y），则参数反演目标函数可表示为

δ = \underset{\begin{matrix} a \leq x \leq b \\ c \leq y \leq d \end{matrix}}{m i n} |f_{h} - F_{h} (x, y)|

（16）

其中：h为模型及实测资料对应的上游水深。

对于均质地工况，需要反演的参数为坝体弹性模量及坝基弹性模量。反演过程中，选取坝体弹性模量变化范围为25 GPa~60 GPa，选取坝基弹性模量为20 GPa~30 GPa。对于分层地基工况，需要反演的参数为坝体弹性模量及坝基2 490 m高程以下区域弹性模量（2 490 m高程以上坝基弹性模量按比例放缩）。反演过程中，选取坝体弹性模量变化范围为25 GPa~60 GPa，对于坝基2 490 m高程以下区域，考虑到坝基分层，故选取其弹性模量变化范围为20 GPa~40 GPa。

计算过程中，分别对上游水位为138，148及158 m 3种情况建模，分别对应低水位、中等水位及高水位3种情况。根据识别结果，3种水位分别对应2016年7月、2013年7月及2014年1月的实测数据。

6种情况下龙羊峡重力拱坝的参数反演结果如表2所示。对于龙羊峡重力拱坝，基于坝体‑地基‑库水流固耦合模型反演得到的坝体弹性模量约为50 GPa~52 GPa；对于均值地基，其弹性模量的反演值为30 GPa；对于分层地基，2 490 m高程以下区域的弹性模量为34 GPa。

表2 参数反演结果

Tab.2 Parameter inversion results

工况	自振频率测量值/Hz	自振频率计算值/Hz	坝体弹性模量/GPa	地基弹摸/GPa
均匀地基，水位为138 m	3.193	3.120	50 50 50	30 30 30
均匀地基，水位为148 m	3.044	3.043
均匀地基，水位为158 m	2.853	2.877
分层地基，水位为138 m	3.193	3.119	52 52 52	34 34 34
分层地基，水位为148 m	3.044	3.033
分层地基，水位为158 m	2.853	2.874

水位为148 m时坝体前2阶频率如表3所示。可见，以第1阶频率为目标进行反演得到坝体及地基材料参数，其第2阶频率也与实测数据基本相符，说明根据实测的坝体第1阶自振频率反演出的坝体和地基弹性模量是可靠的。另外，反演的材料参数与其设计值存在较为明显的提升。

表3 水位为148 m时坝体前2阶频率

Tab.3 First and second order frequency of the dam at the water level of 148 m ( Hz )

阶数	自振频率测量值	自振频率计算值
阶数	自振频率测量值	均值地基	分层地基
第1阶	3.044	3.035	3.033
第2阶	3.809	4.277	4.246

5 坝体‑基岩系统模态分析

根据参数反演结果，计算龙羊峡重力拱坝在各工况下的前3阶阵型。图6~8为均匀地基在水位分别为138，148和158 m时的坝体前3阶振型。图9~11为分层地基在水位分别为138，148和158 m时的坝体前3阶振型。对于这6种工况，龙羊峡重力拱坝第1阶振型和第3阶振型为正对称，第2阶振型为反对称。不同工况下振型幅值位置略有不同，实际上坝体的前3阶振型是一致的。

图6 均匀地基在水位为138 m时的坝体前3阶振型

Fig.6 First three modes of the homogeneous foundation model with 138 m water level

图7 均匀地基在水位为148 m时的坝体前3阶振型

Fig.7 First three modes of the homogeneous foundation model with 148 m water level

图8 均匀地基在水位为158 m时的坝体前3阶振型

Fig.8 First three modes of the homogeneous foundation model with 158 m water level

图9 分层地基在水位为138 m时的坝体前3阶振型

Fig.9 First three modes of the layered foundation model with 138 m water level

图10 分层地基在水位为148 m时的坝体前3阶振型

Fig.10 First three modes of the layered foundation model with 148 m water level

图11 分层地基在水位为158 m时的坝体前3阶振型

Fig.11 First three modes of the layered foundation model with 158 m water level

6 结束语

由坝体实测加速度数据识别出的第1阶自振频率为2.778~3.196 Hz，基于实测数据反演得到的坝体及地基弹性模量显著高于其设计值。采用均质地基模型进行参数反演时，坝体弹性模量反演值为50 GPa，地基弹性模量为30 GPa；采用分层地基模型进行参数反演时，坝体弹性模量反演值为52 GPa，2 490 m高程以下地基弹性模量为34 GPa。龙羊峡重力拱坝第1阶振型和第3阶振型为正对称，第2阶振型为反对称，其各阶频率的分布规律与重力坝类似。

参考文献

孙振刚，张岚，段中德.我国水电站工程数量与规模［J］.中国水利，2013（7）：12-13. [百度学术]

SUN Zhengang， ZHANG Lan， DUAN Zhongde. Quantity and scale of hydropower projects in China［J］. China Water Resources， 2013（7）： 12-13. （in Chinese） [百度学术]

陈厚群.汶川地震后对大坝抗震安全的思考［J］. 中国工程科学，2009，11（6）：44-53. [百度学术]

CHEN Houqun. Consideration on seismic safety of dams in China after the Wenchuan earthquake［J］. Strategic Study of CAE， 2009， 11（6）： 44-53. （in Chinese） [百度学术]

冯新. 土木工程中结构识别方法的研究［D］.大连：大连理工大学，2002. [百度学术]

续秀忠，华宏星，陈兆能.基于环境激励的模态参数辨识方法综述［J］.振动与冲击， 2002， 21（3）： 1-5. [百度学术]

XU Xiuzhong， HUA Hongxing， CHEN Zhaoneng. Review of modal identification method based on ambient excitation［J］. Journal of Vibration and Shock， 2002， 21（3）： 1-5. （in Chinese） [百度学术]

ELLIS R， DURON Z， VON G N， et al. Dynamic characterization of multiple arch dams： a finite element approach［C］∥Conference Proceedings of the Society for Experimental Mechanics Series. Lombard： Springer， 2012： 147-163. [百度学术]

SEVIM B. Finite element model calibration of berke arch dam using operational modal testing［J］. Journal of Vibration and Control， 2011， 17（7）： 1065-1079. [百度学术]

GARCÍA-PALACIOS H， SORIA J M， DÍAZ M， et al. Ambient modal testing of a double-arch dam： the experimental campaign and model updating［J］. Journal of Physics： Conference Series， 2016， 744（1）： 012037. [百度学术]

WENG J， LOH C. Structural health monitoring of arch dam from dynamic measurements［C］∥ASCE Earth and Space Conference. Honolulu： ASCE， 2010： 2518-2534. [百度学术]

李帅，潘坚文，罗广衡，等. 溪洛渡拱坝模态参数识别［J］. 水力发电学报， 2020， 39（3）： 86-95. [百度学术]

LI Shuai， PAN Jianwen， LUO Guangheng， et al. Modal parameter identification for Xiluodu arch dam［J］. Journal of Hydroelectric Engineering， 2020， 39（3）： 86-95. （in Chinese） [百度学术]

冯新，周晶，范颖芳.基于模态观测的混凝土坝反演分析［J］. 水利学报， 2004（2）：101-105. [百度学术]

FENG Xin， ZHOU Jing， FAN Yingfang. Inverse analysis of concrete dam with modal measurements［J］. Journal of Hydraulic Engineering， 2004（2）： 101-105. （in Chinese） [百度学术]

张楚汉，阎承大，王光纶.响洪甸拱坝的动水压力响应分析——考虑水库边界吸收作用的影响［J］.水利学报，1991（11）：26-34. [百度学术]

ZHANG Chuhan， YAN Chengda， WANG Guangguan. Hydrodynamic pressure response analysis of Xianghongdian arch dam：taking account of the effect of reservoir boundary absorption［J］. Journal of Hydraulic Engineering， 1991（11）： 26-34. （in Chinese） [百度学术]

邱奕翔，魏楚函，武志刚，等. 库水模拟对拱坝动力特性的影响分析［J］. 水力发电学报， 2020， 39（6）： 109-120. [百度学术]

QIU Yixiang， WEI Chuhan， WU Zhigang， et al. Effects of reservoir water on simulations of arch dam dynamic characteristics［J］. Journal of Hydroelectric Engineering， 2020， 39（6）： 109-120. （in Chinese） [百度学术]

程琳，王冬冬，杨杰，等. 基于强震观测和多输出支持向量机的混凝土坝材料动参数反演［J］. 振动工程学报， 2017， 30（3）： 466-474. [百度学术]

CHENG Lin， WANG Dongdong， YANG Jie， et al. Back-analysis method based on strong-motion record and multiple-output support vector machine used to determine the dynamic material parameters of concrete dams［J］. Journal of Vibration Engineering， 2017， 30（3）： 466-474. （in Chinese） [百度学术]

李火坤，王刚，魏博文，等. 基于敏感性分析与粒子群算法的拱坝原型动弹性模量反演方法［J］. 水利学报， 2020， 51（11）： 1401-1411. [百度学术]

LI Huokun， WANG Gang， WEI Bowen， et al. Inversion of prototype dynamic elastic modulus of arch dam based on sensitivity analysis and particle swarm optimization［J］. Journal of Hydraulic Engineering， 2020， 51（11）： 1401-1411. （in Chinese） [百度学术]

VAN‑OVERSCHEE P， DEMOOR B. Subspace identification for linear systems： theory-implementation-applications［M］. Netherlands： Kluwer Academic Publi⁃shers， 1996： 135-160. [百度学术]

LI S， WANG J T， JIN A Y， et al. Parametric analysis of SSI algorithm in modal identification of high arch dams［J］. Soil Dynamics and Earthquake Engineering， 2020， 129： 105929. [百度学术]

LI S， PAN J W， LUO G H， et al. Automatic modal parameter identification of high arch dams： feasibility verification［J］. Earthquake Engineering and Engineering Vibration， 2020， 19（4）： 953-965. [百度学术]