CN102520071B

CN102520071B - Transmission tower modal parameter identification method based on improved subspace algorithm

Info

Publication number: CN102520071B
Application number: CN201110428858.4A
Authority: CN
Inventors: 刘晓锋; 卢修连; 黄磊; 孙和泰
Original assignee: NANJING ELECTRIC POWER PLANT OF DATANG GROUP; State Grid Corp of China SGCC; State Grid Jiangsu Electric Power Co Ltd; Jiangsu Fangtian Power Technology Co Ltd
Current assignee: NANJING ELECTRIC POWER PLANT OF DATANG GROUP; State Grid Corp of China SGCC; State Grid Jiangsu Electric Power Co Ltd; Jiangsu Fangtian Power Technology Co Ltd
Priority date: 2011-12-20
Filing date: 2011-12-20
Publication date: 2014-10-15
Anticipated expiration: 2031-12-20
Also published as: CN102520071A

Abstract

The invention discloses a transmission tower modal parameter identification method based on an improved subspace algorithm. For a large-span transmission tower structure, under the excitation of the field environment, the modal parameter identification of the structure is completed solely by using response data. The identification process is based on the principle of random subspace identification, and realizes the integration of multi-group group response data into synchronous impulse response data. The integrated data avoids the fitting error caused by group identification and realizes the synchronization of data as a whole. Identification, and finally use the modal parameter identification results under different system orders to construct a stability diagram, which solves the order determination problem of modal parameter identification under environmental excitation, avoids the modal omission and repetition phenomenon of traditional methods such as peak method, and gives More accurate and stable identification results, and lay the foundation for subsequent structural damage detection and life estimation.

Description

Modal Parameter Identification Method of Transmission Tower Based on Improved Subspace Algorithm

技术领域 technical field

本发明涉及一种大型结构试验模态分析技术，尤其是一种输电塔模态参数识别方法。The invention relates to a large-scale structure test modal analysis technology, in particular to a transmission tower modal parameter identification method.

背景技术 Background technique

大型结构试验模态分析技术是新世纪开展起来的结构健康安全评估的有效方法之一。针对输电线路铁塔原型动力特性与动力响应现场实测是近年才开展起来的。由于试验人员在输电塔原型高压电场附近高空作业，测量工作条件十分困难；微风时振动信号微弱，对测量传感器及仪器的信噪比要求高；整塔模态试验为了获得较可靠的振形，需要较多的测点，对传感器及数采仪要求通道数很多。因此，目前有关输电塔原型的振动模态现场实测据甚少，整塔的试验模态分析工作几乎是空白，个别研究工作只做了一阶弯曲模态。Large-scale structural test modal analysis technology is one of the effective methods for structural health and safety assessment developed in the new century. The on-site measurement of the dynamic characteristics and dynamic response of the transmission line tower prototype has only been carried out in recent years. Because the testers work at high altitude near the high-voltage electric field of the transmission tower prototype, the measurement working conditions are very difficult; the vibration signal is weak in the breeze, and the signal-to-noise ratio of the measurement sensor and instrument is high; the whole tower modal test In order to obtain a more reliable vibration shape, More measuring points are required, and a large number of channels are required for sensors and data acquisition instruments. Therefore, at present, there is very little field measurement data on the vibration mode of the prototype of the transmission tower, and the experimental modal analysis of the whole tower is almost blank, and only the first-order bending mode has been done in individual research work.

对桥梁、建筑物、输电塔等大型结构进行模态试验，无论是正弦、随机或者脉冲方式的人为激励都是不可能或不允许的。但是任何大型结构物都存在一定的振动环境，例如风、水流冲击、大地脉动、移动的车辆引起的振动等，在这些自然环境的激励下，结构物都会产生微弱的振动。在强风作用，输电铁塔的环境振动还可能是很大的。虽然我们对这些激励特性无法精确定量，但并非一无所知，可合理的假定这些激励是近似的平稳随机信号，其频谱是具有一定带宽的连续谱，在带宽内基本覆盖了对结构物感兴趣的频带，从而在结构物的自然环境激励下的振动信号中包含了这些模态。It is impossible or not allowed to conduct modal tests on large structures such as bridges, buildings, and transmission towers, whether it is sinusoidal, random, or pulsed artificial excitation. However, any large structure has a certain vibration environment, such as wind, water impact, ground pulsation, vibration caused by moving vehicles, etc. Under the excitation of these natural environments, the structure will produce weak vibrations. Under the effect of strong wind, the environmental vibration of the transmission tower may be very large. Although we cannot precisely quantify the characteristics of these excitations, we are not ignorant of them. It is reasonable to assume that these excitations are approximately stationary random signals, and their spectrum is a continuous spectrum with a certain bandwidth. The frequency band of interest, so that these modes are included in the vibration signal excited by the structure's natural environment.

对结构物振动信号中的振型、阻尼进行识别，可以为后续的有限元建模、损伤和安全寿命估计等提供可靠的模态参数识别结果。目前流行且在工程上比较有效地获得的自然环境激励下的振型的方法是将全部测点在环境激励下的振动响应和某一固定的参考点的振动响应分别作双通道频谱分析。首先在自功率谱图上识别出共振频率f_i，再将各测点与参考点在共振频率上的幅值谱之比Φ(f_i)作为该点的振型的相对值，将它们的互功率谱的实部在此频率上的正负作为该点振型的相位。Identifying the mode shape and damping in the vibration signal of the structure can provide reliable modal parameter identification results for subsequent finite element modeling, damage and safety life estimation. The currently popular and more effective method to obtain the mode shape under natural environment excitation is to perform dual-channel spectrum analysis on the vibration response of all measuring points under environmental excitation and the vibration response of a fixed reference point. First identify the resonant frequency f _i on the self-power spectrum diagram, then take the ratio Φ(f _i ) of the amplitude spectrum at the resonant frequency of each measuring point and the reference point as the relative value of the mode shape at this point, and take their The positive and negative values of the real part of the cross power spectrum at this frequency are taken as the phase of the point mode shape.

$| | Φ Φ (({f f}_{i i})) | | = = | | \frac{B B (({f f}_{i i}))}{A A (({f f}_{i i}))} | | = = {| | \frac{B B (({f f}_{i i})) \overset{&OverBar; &OverBar;}{B B} (({f f}_{i i}))}{A A (({f f}_{i i})) \overset{&OverBar; &OverBar;}{A A} (({f f}_{i i}))} | |}^{\frac{11}{22}} = = {| | \frac{{G G}_{bb bb} (({f f}_{i i}))}{{G G}_{aa aa} (({f f}_{i i}))} | |}^{\frac{11}{22}} - - - - - - ((11))$

sgn(Φ(f_i))＝sgn(Real(G_ba(f_i)) (2)sgn(Φ(f _i ))＝sgn(Real(G _ba (f _i )) (2)

其中A(f_i)为参考点信号的傅氏变换，B(f_i)为测量点信号的傅氏变换，G_aa(f_i)、G_bb(f_i)分别为参考点信号、测量点信号的自功率谱，G_ba(f_i)为测量点信号与参考点信号的互功率谱。Among them, A(f _i ) is the Fourier transform of the reference point signal, B(f _i ) is the Fourier transform of the measurement point signal, G _aa (f _i ), G _bb (f _i ) are the reference point signal, the measurement point The self-power spectrum of the signal, G _ba (f _i ) is the cross-power spectrum of the measurement point signal and the reference point signal.

环境激励的模态识别方法单独利用响应数据进行模态参数识别，不需人为激励，从而大大节省了试验设备的成本。并且避免了人为激励，如发射火箭或切断预拉力钢丝，造成突发冲击等引起的结构损伤。而且方法在结构物真实支承状态下进行，便于与理论计算进行比较。The modal identification method of environmental excitation uses the response data alone to identify the modal parameters without artificial excitation, thus greatly saving the cost of test equipment. And it avoids the structural damage caused by artificial excitation, such as launching a rocket or cutting off the pre-tension steel wire, causing sudden impact and so on. Moreover, the method is carried out under the real support state of the structure, which is convenient for comparison with theoretical calculations.

近年来关于环境激励下工程结构的模态参数识别方法主要分为时域分析和频域分析。In recent years, the modal parameter identification methods of engineering structures under environmental excitation are mainly divided into time domain analysis and frequency domain analysis.

频域识别方法主要有峰值发、频域分解法。这两种方法具有快速、方便和经济的优点，但是前者峰值的拾取往往是主观的，得到的是结构工作挠曲形状而不是振型，仅限于实模态和比例阻尼结构，不能识别复杂模态和阻尼比。后者是峰值拾取法的延伸，克服了峰值拾取法的不足，但不能识别密集模态。Frequency domain identification methods mainly include peak detection and frequency domain decomposition methods. These two methods have the advantages of fast, convenient and economical, but the peak picking of the former is often subjective, and the working deflection shape of the structure is obtained instead of the mode shape, which is limited to the real mode and proportional damping structure, and cannot identify complex modes. state and damping ratio. The latter is an extension of the peak-picking method, which overcomes the shortcomings of the peak-picking method, but cannot identify dense modes.

目前，大型结构环境激励模态参数的识别一般采用频域方法，各阶模态阻尼根据全部响应点信号的集总平均谱，人为进行模态定阶。At present, the frequency domain method is generally used to identify the excitation modal parameters of large-scale structural environments, and the modal order of each order of modal damping is artificially determined according to the aggregated average spectrum of all response point signals.

传统的频域识别算法，如：峰值法、NExT算法等可以较为准确的识别固有频率值，但是在环境激励的条件下很难有效确定模态阶数，很容易受到周期信号干扰，而且振型的识别精度较差。Traditional frequency domain identification algorithms, such as: peak value method, NExT algorithm, etc., can identify the natural frequency value more accurately, but it is difficult to effectively determine the mode order under the condition of environmental excitation, and it is easy to be interfered by periodic signals, and the mode shape The recognition accuracy is poor.

随机子空间识别方法是1995年以来国内外模态分析的一个热点。它充分利用了矩阵QR分解，奇异值分解SVD，最小二乘法等强大的数学工具，使得该算法非常有效的进行环境激励下参数识别，是目前最先进的结构环境振动模态参数识别方法。The random subspace identification method has been a hot spot in modal analysis at home and abroad since 1995. It makes full use of powerful mathematical tools such as matrix QR decomposition, singular value decomposition SVD, least square method, etc., which makes the algorithm very effective for parameter identification under environmental excitation. It is currently the most advanced method for identifying structural environment vibration modal parameters.

由于随机子空间识别方法是一种时域算法，当结构所有测试点数据同步采集、测量时对于识别精度的提高有很好的效果。但是由于受到测试条件的限制，大型输电塔框架尺寸较大，为了表征系统的固有特性，需要很多测量自由度。而现场难于多测点同时测量，一般采用，固定参考点，移动测量点的试验方法。当前通常采用的方法是选择好固定的参考点位置后，分批次逐步完成整个结构的测试。但分批次测量后存在数据整合的问题，难以避免整合过程带来的拟合误差。Since the random subspace identification method is a time-domain algorithm, when the data of all test points of the structure are collected and measured synchronously, it has a good effect on improving the identification accuracy. However, due to the limitation of test conditions, large transmission tower frame size is large, in order to characterize the inherent characteristics of the system, many measurement degrees of freedom are required. However, it is difficult to measure multiple measuring points at the same time in the field. Generally, the test method of fixing the reference point and moving the measuring point is adopted. The current method usually used is to gradually complete the testing of the entire structure in batches after selecting a fixed reference point position. However, there is a problem of data integration after batch measurement, and it is difficult to avoid the fitting error caused by the integration process.

发明内容Contents of the invention

本发明的目的就是为了解决克服现有基于环境激励下频域识别算法的局限性以及大型结构无法同步测量所有测试点的困难，对现有随机子空间模态参数识别算法进行改进，提高模态参数的识别精度，改善人为模态定阶的缺陷。The purpose of the present invention is to solve the limitations of the existing frequency-domain identification algorithm based on environmental excitation and the difficulty that all test points cannot be measured synchronously for large structures, improve the existing random subspace modal parameter identification algorithm, and improve the modal The identification accuracy of parameters improves the defects of artificial mode order determination.

本发明的技术方案是提供一种基于改进子空间算法的输电塔模态参数识别方法，其特征在于：其包括以下步骤：The technical solution of the present invention is to provide a transmission tower modal parameter identification method based on the improved subspace algorithm, which is characterized in that: it includes the following steps:

1)在输电塔上设置一个参考点和多个测量点，分步采集环境激励下参考点和测量点的振动响应信号，计算其自功率谱和互功率谱；1) Set a reference point and multiple measurement points on the transmission tower, collect the vibration response signals of the reference point and measurement points under environmental excitation step by step, and calculate their self-power spectrum and cross-power spectrum;

2)由该互功率谱得到互相关序列组成hankel矩阵；2) Obtaining a cross-correlation sequence from the cross-power spectrum to form a Hankel matrix;

3)计算该hankel矩阵的投影矩阵；3) Calculate the projection matrix of the Hankel matrix;

4)对该投影矩阵加权处理；4) weighting the projection matrix;

5)对步骤4)中加权处理后的投影矩阵进行奇异值分解；5) carrying out singular value decomposition to the weighted projection matrix in step 4);

6)根据随机子空间理论计算输电塔系统的离散状态空间方程的卡尔曼滤波状态向量；6) Calculate the Kalman filter state vector of the discrete state space equation of the transmission tower system according to the stochastic subspace theory;

7)根据输电塔系统离散状态空间方程求解系统状态矩阵和输出矩阵；7) Solve the system state matrix and output matrix according to the discrete state space equation of the transmission tower system;

8)根据系统状态矩阵和输出矩阵求解系统模态参数；8) Solve the system modal parameters according to the system state matrix and output matrix;

9)重复步骤7)和8)，获得不同阶次的输电塔系统的模态参数，建立系统模态参数的稳定图。9) Repeat steps 7) and 8) to obtain the modal parameters of the transmission tower system of different orders, and establish a stability diagram of the system modal parameters.

优选的，所述步骤2)中组成的hankel矩阵为：Preferably, the hankel matrix formed in the step 2) is:

该Hankel矩阵是由测量点的响应数据y_k组成的2mi×j矩阵，y_i表示i时刻所有测量点的响应，下标p表示过去，下标f表示将来，Y_p＝Y_0|i-1为Hankel矩阵第一列中的下标起始为0，终点为i-1的元素对应的所有的行和列组成的Hankel矩阵的块，表示过去行空间，Y_f＝Y_i-1|2i-1为Hankel矩阵的后i行，表示将来行空间。The Hankel matrix is a 2mi×j matrix composed of response data y _k of measurement points, y _i represents the response of all measurement points at time i, subscript p represents the past, subscript f represents the future, Y _p =Y _{0|i- 1} is the block of the Hankel matrix composed of all rows and columns corresponding to the subscript in the first column of the Hankel matrix whose subscript starts from 0 and ends at i-1, indicating the past row space, Y _f =Y _{i-1| 2i-1} is the last i row of the Hankel matrix, representing the future row space.

优选的，所述步骤3)中求解hankel矩阵的投影矩阵的步骤为：Preferably, the step of solving the projection matrix of the Hankel matrix in the step 3) is:

3.1)求振动响应信号的将来行空间到过去行空间的投影O_i，求解公式为：3.1) Find the projection O _i from the future row space to the past row space of the vibration response signal, and the solution formula is:

${O o}_{i i} = = {Y Y}_{f f} / / {Y Y}_{p p} = = {Y Y}_{f f} {Y Y}_{p p}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{+ +} {Y Y}_{p p},,$

式中，(·)⁺表示矩阵求广义逆，(·)^T表示矩阵求转置；In the formula, ( ) ⁺ means seeking generalized inversion of matrix, ( ) ^T means seeking transpose of matrix;

3.2)根据随机子空间理论，所识别输电塔系统的离散状态空间方程为：3.2) According to the stochastic subspace theory, the discrete state space equation of the identified transmission tower system is:

x_k+1＝Ax_k+ρ_w x _k+1 ＝Ax _k +ρ _w

y_k＝Cx_k+ρ_v y _k ＝Cx _k +ρ _v

式中A，C分别表示n×n阶系统状态矩阵和m×n阶输出矩阵，该离散状态空间模型的阶数为n，ρ_w，ρ_v分别为过程噪音和测量噪音；x_k为k时刻的系统状态向量，y_k为系统的k时刻的振动响应向量；系统模态参数由系统状态矩阵A的特征值和特征向量，以及系统输出矩阵C表示。In the formula, A and C represent the n×n order system state matrix and m×n order output matrix respectively, the order of the discrete state space model is n, ρ _w , ρ _v are process noise and measurement noise respectively; x _k is k The system state vector at time, y _k is the vibration response vector of the system at time k; the system modal parameters are represented by the eigenvalues and eigenvectors of the system state matrix A, and the system output matrix C.

3.3)根据随机子空间识别理论，投影矩阵O_i可以分解为可控可观矩阵Γ_i和卡尔曼滤波状态向量的乘积，即：3.3) According to the random subspace identification theory, the projection matrix O _i can be decomposed into a controllable observable matrix Γ _i and a Kalman filter state vector The product of , that is:

${O o}_{i i} = = {Γ Γ}_{i i} {\overset{^^}{X x}}_{i i} = = [\begin{matrix} C C \\ CA CA \\ {CA CA}^{22} \\ \cdot \cdot \cdot \cdot \cdot \cdot \\ {CA CA}^{i i - - 11} \end{matrix}] [\begin{matrix} {\overset{^^}{x x}}_{i i} & {\overset{^^}{x x}}_{i i + + 11} & {\overset{^^}{x x}}_{i i + + 22} & {\overset{^^}{x x}}_{i i + + j j - - 11} \end{matrix}]$

式中，为卡尔曼滤波状态向量，是系统状态向量x_k的最优估计值；Γ_i为由系统状态矩阵A和输出矩阵C构成的可控可观矩阵。In the formula, is the Kalman filter state vector, which is the optimal estimated value of the system state vector x _k ; Γ _i is a controllable and observable matrix composed of the system state matrix A and the output matrix C.

优选的，步骤4)中采用主分量PC方法对投影矩阵进行加权处理，其方法是：首先，设加权矩阵为：Preferably, in step 4), adopt principal component PC method to carry out weighting process to projection matrix, and its method is: at first, set weighting matrix as:

W₁＝IW ₁ =I

${W W}_{22} = = {Y Y}_{P P}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{11 / / 22} {Y Y}_{p p}$

式中，I为单位矩阵；In the formula, I is the identity matrix;

代入投影矩阵O_i，获得加权后的投影矩阵为：Substituting into the projection matrix O _i , the weighted projection matrix obtained is:

${O o}_{i i}^{' '} = = {W W}_{11} {O o}_{i i} {W W}_{22} = = {Y Y}_{f f} {Y Y}_{P P}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{- - 11 / / 22} {Y Y}_{p p} . .$

优选的，步骤5)对步骤4)中加权处理后的投影矩阵进行奇异值分解的步骤为：Preferably, step 5) the step of carrying out singular value decomposition to the weighted projection matrix in step 4) is:

5.1)对投影矩阵进行奇异分解，确定投影矩阵的系统阶数：5.1) Perform singular decomposition on the projection matrix to determine the system order of the projection matrix:

${O o}_{i i}^{' '} = = {W W}_{11} {O o}_{i i} {W W}_{22} = = UΣV UΣV = = [\begin{matrix} {U u}_{r r} & {U u}_{s the s} \end{matrix}] [\begin{matrix} {S S}_{r r} \\ {S S}_{s the s} = = 00 \end{matrix}] [\begin{matrix} {V V}_{r r} \\ {V V}_{s the s} \end{matrix}] {{U u}_{r r}^{T T} S S}_{r r} {V V}_{r r}^{T T},,$

式中，U_r，S_r，V_r为主分量的左奇异值向量、奇异值和右奇异值向量；U_s，S_s，V_s为噪声分量的左奇异值向量、奇异值和有奇异值向量；U，V为正交矩阵；S为由从大到小排列的奇异值组成的对角矩阵；In the formula, U _r , S _r , V _r are the left singular value vector, singular value and right singular value vector of the main component; U _s , S _s , V _s are the left singular value vector, singular value and singular value of the noise component Value vector; U, V are orthogonal matrices; S is a diagonal matrix composed of singular values arranged from large to small;

根据投影矩阵的奇异值分解得到以下结果：According to the singular value decomposition of the projection matrix, the following results are obtained:

1)加权后的投影矩阵O′_i等于可控可观矩阵Γ_i和卡尔曼滤波状态向量的乘积： 1) The weighted projection matrix O′ _i is equal to the controllable observable matrix Γ _i and the Kalman filter state vector The product of:

2)投影矩阵的系统阶数为公式中奇异值矩阵中不为零的奇异值数；2) The system order of the projection matrix is the number of singular values that are not zero in the singular value matrix in the formula;

3)可控可观矩阵为： 3) The controllable and observable matrix is:

4)i时刻的卡尔曼滤波状态向量为：4) The Kalman filter state vector at time i is:

5.2)根据权利要求2所述的Hankel矩阵的第二种分块方法，投影矩阵用下式表达：5.2) according to the second kind of blocking method of Hankel matrix described in claim 2, projection matrix is expressed with following formula:

${O o}_{i i - - 11}^{' '} = = {Y Y}_{f f}^{- -} / / {Y Y}_{p p}^{+ +} = = {Γ Γ}_{i i - - 11} \cdot &Center Dot; {\overset{^^}{X x}}_{i i + + 11};;$

则i+1时刻的卡尔曼滤波状态向量为： Then the Kalman filter state vector at time i+1 is:

Γ_i-1为Γ_i剔除最后一块行所得的矩阵。Γ _i-1 is the matrix obtained by excluding the last row of Γ _i .

优选的，由可知，卡尔曼滤波状态向量仅用振动响应数据即可得到；preferred by It can be seen that the Kalman filter state vector It can be obtained only by vibration response data;

最后将求得的卡尔曼滤波状态向量代入输电塔系统离散状态空间方程得： $(\begin{matrix} {\hat{X}}_{i + 1} \\ Y_{i | i} \end{matrix}) = (\begin{matrix} A \\ C \end{matrix}) ({\hat{X}}_{i}) + (\begin{matrix} ρ_{w} \\ ρ_{v} \end{matrix}),$ 式中，ρ_w，ρ_v为卡尔曼滤波最佳估计的残差；Y_i|i表示Hankel矩阵第一列中的下标起始为i，终点为i的元素对应的所有的行和列组成的Hankel矩阵的块，表示只有一个块行的Hankel矩阵。Finally, substitute the obtained Kalman filter state vector into the discrete state space equation of the transmission tower system to get: $(\begin{matrix} {\hat{x}}_{i + 1} \\ Y_{i | i} \end{matrix}) = (\begin{matrix} A \\ C \end{matrix}) ({\hat{x}}_{i}) + (\begin{matrix} ρ_{w} \\ ρ_{v} \end{matrix}),$ In the formula, ρ _w , ρ _v are the residuals of the best estimate of the Kalman filter; Y _i|i represents all the rows and columns corresponding to the elements whose subscript starts at i and ends at i in the first column of the Hankel matrix The blocks that make up the Hankel matrix represent a Hankel matrix with only one block row.

优选的，由于卡尔曼滤波状态向量和振动响应数据已知，ρ_ω，ρ_v和相互独立，利用数值计算领域广泛应用的最小二乘法求解输电塔系统离散状态空间方程，得到线性方程：由此方程求解出系统状态矩阵A和输出矩阵C。Preferably, since the Kalman filter state vector and vibration response data are known, ρ _ω , ρ _v and Independently of each other, the discrete state space equation of the transmission tower system is solved by using the least square method widely used in the field of numerical calculation, and the linear equation is obtained: From this equation, the system state matrix A and output matrix C are solved.

优选的，所述步骤6)还包括以下步骤：Preferably, said step 6) also includes the following steps:

6.1)对系统状态矩阵A进行特征值分解：6.1) Decompose the eigenvalue of the system state matrix A:

A＝ψΛψ^-1 A＝ψΛψ ^-1

式中：Λ＝diag(λ_i)，是系统状态矩阵的复特征值的对角阵，ψ是特征向量矩阵；In the formula: Λ=diag(λ _i ), is the diagonal matrix of the complex eigenvalues of the system state matrix, and ψ is the eigenvector matrix;

6.2)根据系统状态矩阵的特征值求系统的固有频率和阻尼比：6.2) Calculate the natural frequency and damping ratio of the system according to the eigenvalues of the system state matrix:

其中，系统状态矩阵的复特征值λ_i，是λ_i的共轭复数。Among them, the complex eigenvalue λ _i of the system state matrix, is the complex conjugate number of λ _i .

根据随机子空间理论，对于单位采样间隔时间的输电塔系统，系统第i阶固有频率 According to the stochastic subspace theory, for the transmission tower system with unit sampling interval time, the i-th order natural frequency of the system is

系统第i阶阻尼比 The i-th order damping ratio of the system

8.3)由输出矩阵C和系统状态矩阵的调整向量得到系统第i阶振型为： 8.3) From the output matrix C and the adjustment vector of the system state matrix, the i-th order mode shape of the system is obtained as:

优选的，步骤9)中建立稳定图并识别系统稳定、正确模态参数的方法为：Preferably, in step 9), the method for establishing a stable diagram and identifying system stability and correct modal parameters is:

设所识别输电塔系统的离散状态空间方程阶次为n，得到多个不同阶次的离散状态空间方程，对每个方程进行模态参数识别，将得到的所有模态参数绘制在二维坐标图中，坐标图的横坐标为频率值，纵坐标为离散方程的阶数，便可得到稳定图；Assuming that the order of the discrete state space equation of the identified transmission tower system is n, multiple discrete state space equations of different orders are obtained, and the modal parameters of each equation are identified, and all the obtained modal parameters are drawn on the two-dimensional coordinates In the figure, the abscissa of the coordinate diagram is the frequency value, and the ordinate is the order of the discrete equation, and the stability diagram can be obtained;

稳定图把不同阶次方程的模态参数画在同一幅图上，在相应于系统某阶固有频率的轴上，高一阶状态空间方程识别的模态参数同低一阶状态空间方程识别的模态参数相比较，如果固有频率、阻尼比和模态振型的差异小于预设的限定值，则这个点就称为稳定点，稳定点组成的轴为稳定轴，相应的模态参数即为系统的模态参数；限定值取固有频率＜1％，阻尼比＜5％，模态振型＜2％，即稳定轴要求为：The stability diagram draws the modal parameters of different order equations on the same graph. On the axis corresponding to a certain order natural frequency of the system, the modal parameters identified by the higher-order state-space equations are the same as those identified by the lower-order state-space equations. Compared with the parameters, if the difference of natural frequency, damping ratio and mode shape is less than the preset limit value, then this point is called stable point, the axis composed of stable points is the stable axis, and the corresponding modal parameters are the system modal parameters; the limit values are natural frequency < 1%, damping ratio < 5%, mode shape < 2%, that is, the stable axis requirements are:

$\frac{| | {ω ω}^{((n no))} - - {ω ω}^{((n no + + 11))} | |}{{f f}^{((n no))}} \times \times 100100 % % < < 11 % %$

$\frac{| | {ξ ξ}^{((n no))} - - {ξ ξ}^{((n no + + 11))} | |}{{ξ ξ}^{((n no))}} \times \times 100100 % % < < 55 % %$

(1-MAC(n，n+1))×100％＜2％(1-MAC(n,n+1))×100%＜2%

式中n表示输电塔系统离散状态空间方程的阶次，阶次为n的状态空间方程计算得到系统的模态参数为：固有频率ω⁽ⁿ⁾，阻尼比ξ⁽ⁿ⁾，模态振型MAC表示模态保证准则：In the formula, n represents the order of the discrete state space equation of the transmission tower system, and the modal parameters of the system obtained by calculating the state space equation with the order n are: natural frequency ω ⁽ⁿ⁾ , damping ratio ξ ⁽ⁿ⁾ , mode shape MAC represents the Modality Assurance Criteria:

优选的，所述步骤2)中设参考点为j，响应点数为i，两点之间的互相关函数的表达式为：Preferably, in said step 2), set the reference point as j, the number of response points as i, and the expression of the cross-correlation function between the two points is:

${R R}_{ij ij} = = {Σ Σ}_{r r = = 11}^{N N} {A A}_{ijr ijr} exp exp ((- - {ξ ξ}_{r r} {ω ω}_{nr nr} t t)) sin sin (({ω ω}_{dr dr} t t + + θ θ));;$

A_cmr为幅值，ξ_r为阻尼比，ω_nr为固有频率，ω_dr共振频率，r为模态阶数，N为总的模态阶数；A _cmr is the amplitude, ξ _r is the damping ratio, ω _nr is the natural frequency, ω _dr is the resonance frequency, r is the modal order, and N is the total modal order;

用该互相关函数代替振动响应信号，对步骤1)中计算得到的互功率谱进行反傅里叶变换，得到互相关函数：R_ij＝ifft(F_ij)，式中F_ij为互功率谱，由响应点i和参考点j的自功率谱计算得到；Use this cross-correlation function to replace the vibration response signal, perform inverse Fourier transform on the cross-power spectrum calculated in step 1), and obtain the cross-correlation function: R _ij =ifft(F _ij ), where F _ij is the cross-power spectrum , calculated from the autopower spectrum of response point i and reference point j;

建立所述Hankel矩阵：Build the Hankel matrix:

式中，r_i为互相关函数矩阵中的第i列，与i时刻的响应信号相对应，R_p＝R_0|i-1为Hankel矩阵的前i行，表示过去行空间，R_f＝R_i-1|2i-1为Hankel矩阵的后i行，表示将来行空间。In the formula, r _i is the i-th column in the cross-correlation function matrix, corresponding to the response signal at time i, R _p =R _0|i-1 is the first i row of the Hankel matrix, representing the past row space, R _f = R _i-1|2i-1 is the last i row of the Hankel matrix, representing the future row space.

本发明的基于改进子空间算法的输电塔模态参数识别方法对大跨越输电塔结构，在野外环境激励下，单独利用响应数据完成结构的模态参数识别。识别过程在随机子空间识别原理的基础之上，实现了将多组分步响应数据整合为同步的脉冲响应数据，整合后的数据避免了分组识别带来的拟合误差，实现了数据的同步整体识别，最后利用不同阶数下的识别结果构建稳定图，解决了环境激励下模态参数识别的定阶问题，避免了峰值法等传统方法的模态遗漏和重复现象，给出更加准确稳定的识别结果，并为后续的结构损伤检测和寿命估计奠定了基础。The modal parameter identification method of the transmission tower based on the improved subspace algorithm of the present invention, for the large-span transmission tower structure, under the excitation of the field environment, independently uses the response data to complete the modal parameter identification of the structure. The identification process is based on the principle of random subspace identification, and realizes the integration of multi-group step response data into synchronous impulse response data. The integrated data avoids the fitting error caused by group identification and realizes data synchronization. Overall identification, and finally use the identification results under different orders to construct a stability diagram, which solves the order determination problem of modal parameter identification under environmental excitation, avoids the modal omission and repetition phenomenon of traditional methods such as peak method, and gives a more accurate and stable The identification results of the structure lay the foundation for the subsequent structural damage detection and life estimation.

附图说明 Description of drawings

图1是本发明的基于改进子空间算法的输电塔模态参数识别方法的流程图；Fig. 1 is the flow chart of the transmission tower modal parameter identification method based on the improved subspace algorithm of the present invention;

图2是稳定图生成原理图。Figure 2 is a schematic diagram of the stable map generation.

具体实施方式 Detailed ways

下面对本发明的具体实施方式作进一步详细的描述。Specific embodiments of the present invention will be further described in detail below.

如图1和图2所示，本发明的一种基于改进子空间算法的输电塔模态参数识别方法，其包括以下步骤：As shown in Figure 1 and Figure 2, a kind of transmission tower modal parameter identification method based on improved subspace algorithm of the present invention, it comprises the following steps:

1、在输电塔上设置参考点和多个测量点，采集输电塔在环境激励下的振动响应信号，并计算其互动功率谱。1. Set reference points and multiple measurement points on the transmission tower, collect vibration response signals of the transmission tower under environmental excitation, and calculate its interactive power spectrum.

2、应用随机子空间算法理论求解所有测点的响应的投影矩阵2. Apply the random subspace algorithm theory to solve the projection matrix of the responses of all measuring points

1)首先是建立Hankel矩阵，将测点振动响应数据y(k)组成2m×j的Hankel矩阵(Hankel矩阵是交叉对角线上具有元素相同的矩阵)。把Hanke l矩阵的行空间分成“过去”行空间和“将来”行空间：1) Firstly, the Hankel matrix is established, and the vibration response data y(k) of the measuring point is composed into a 2m×j Hankel matrix (the Hankel matrix is a matrix with the same elements on the cross diagonal). Divide the row space of the Hanke l matrix into "past" and "future" row spaces:

y_i表示i时刻所有测量点的响应，下标p表示过去，下标f表示将来，Y_p＝Y_0|i-1为Hankel矩阵第一列中的下标起始为0，终点为i-1的元素对应的所有的行和列组成的Hankel矩阵的块，表示过去行空间，Y_f＝Y_i-1|2i-1为Hankel矩阵的后i行，表示将来行空间。y _i represents the response of all measurement points at time i, the subscript p represents the past, the subscript f represents the future, Y _p =Y _0|i-1 is the subscript in the first column of the Hankel matrix, starting from 0 and ending at i The block of the Hankel matrix composed of all rows and columns corresponding to the element of -1 represents the past row space, and Y _f =Y _i-1|2i-1 is the last i row of the Hankel matrix, representing the future row space.

Hankel矩阵的另外一种划分方法是，将“将来”行空间的第一块行移到“过去”行空间中，即“过去”的块行数为i+1，“将来”的块行数为i-1。Another division method of the Hankel matrix is to move the first row of the "future" row space to the "past" row space, that is, the number of block rows in the "past" is i+1, and the number of block rows in the "future" for i-1.

$H h = = {Y Y}_{00 | | 22 i i - - 11} = = = = [[\frac{{Y Y}_{00 | | i i}}{{Y Y}_{i i + + 11 | | 22 i i - - 11}}]] = = [[\frac{{Y Y}_{p p}^{+ +}}{{Y Y}_{f f}^{- -}}]] - - - - - - ((22))$

2)求响应的将来行空间到过去行空间的投影2) Find the projection of the corresponding future row space to the past row space

${O o}_{i i} = = {Y Y}_{f f} / / {Y Y}_{p p} = = {Y Y}_{f f} {Y Y}_{p p}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{+ +} {Y Y}_{p p} - - - - - - ((33))$

其中(·)⁺表示求广义逆，(·)^T表示求转置。Among them (·) ⁺ means seeking generalized inverse, (·) ^T means seeking transpose.

3)根据随机子空间理论，所识别输电塔系统的离散状态空间方程表示为：3) According to the stochastic subspace theory, the discrete state space equation of the identified transmission tower system is expressed as:

x_k+1＝Ax_k+ρ_w (4)x _k+1 ＝Ax _k +ρ _w (4)

y_k＝Cx_k+ρ_v y _k ＝Cx _k +ρ _v

式中，A，C分别表示n×n阶系统状态矩阵和m×n阶和输出矩阵，离散状态空间方程的阶次为n，ρ_w，ρ_v分别为过程噪音和测量噪音。x_k为k时刻的系统状态向量，y_k为系统的k时刻的振动响应向量。系统的模态参数由系统状态矩阵A的特征值和特征向量，以及系统输出矩阵C表示。In the formula, A and C represent the n×n order system state matrix and the m×n order sum output matrix respectively, the order of the discrete state space equation is n, ρ _w , ρ _v are process noise and measurement noise respectively. x _k is the system state vector at time k, and y _k is the vibration response vector of the system at time k. The modal parameters of the system are represented by the eigenvalues and eigenvectors of the system state matrix A and the system output matrix C.

根据随机子空间识别理论，投影矩阵O_i可以分解为可观矩阵Γ_i和卡尔曼滤波状态向量的乘积。According to the stochastic subspace identification theory, the projection matrix O _i can be decomposed into an observable matrix Γ _i and a Kalman filter state vector product of .

${O o}_{i i} = = {Γ Γ}_{i i} {X x}_{i i} = = [\begin{matrix} C C \\ CA CA \\ {CA CA}^{22} \\ \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; \\ {CA CA}^{i i - - 11} \end{matrix}] [\begin{matrix} {\overset{^^}{x x}}_{i i} & {\overset{^^}{x x}}_{i i + + 11} & {\overset{^^}{x x}}_{i i + + 22} & {\overset{^^}{x x}}_{i i + + j j - - 11} \end{matrix}] - - - - - - ((55))$

式中，为卡尔曼滤波滤波状态向量，是系统状态向量x_k的最优估计值；Γ_i为由系统状态矩阵A和输出矩阵C构成的可控可观矩阵。In the formula, is the Kalman filter filter state vector, which is the optimal estimated value of the system state vector x _k ; Γ _i is a controllable and observable matrix composed of the system state matrix A and the output matrix C.

3、投影矩阵的模态参数识别3. Modal parameter identification of projection matrix

随机子空间模态参数识别方法通过对投影矩阵采用不同的加权处理方式，将产生不同的识别算法。采用主分量PC方法对投影矩阵进行加权处理：The random subspace modal parameter identification method will produce different identification algorithms by adopting different weighting methods for the projection matrix. The projection matrix is weighted using the principal component PC method:

设加权矩阵为：Let the weighting matrix be:

W₁＝IW ₁ =I

(6)(6)

式中，I为单位矩阵；代入投影矩阵O_i，获得加权后的投影矩阵为：In the formula, I is the unit matrix; substituting into the projection matrix O _i , the weighted projection matrix obtained is:

${O o}_{i i}^{' '} = = {W W}_{11} {O o}_{i i} {W W}_{22} = = {Y Y}_{f f} {Y Y}_{P P}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{- - 11 / / 22} {Y Y}_{p p} . . - - - - - - ((77))$

对加权处理后的投影矩阵进行奇异值分解，确定投影矩阵的系统阶数：Singular value decomposition is performed on the weighted projection matrix to determine the system order of the projection matrix:

${O o}_{i i}^{' '} = = {W W}_{11} {O o}_{i i} {W W}_{22} = = UΣV UΣV = = [\begin{matrix} {U u}_{r r} & {U u}_{s the s} \end{matrix}] [\begin{matrix} {S S}_{r r} \\ {S S}_{s the s} = = 00 \end{matrix}] [\begin{matrix} {V V}_{r r} \\ {V V}_{s the s} \end{matrix}] {{U u}_{r r}^{T T} S S}_{r r} {V V}_{r r}^{T T},, - - - - - - ((88))$

A.加权后的投影矩阵O′_i等于可控可观矩阵Γ_i和卡尔曼滤波状态向量的乘积： A. The weighted projection matrix O′ _i is equal to the controllable observable matrix Γ _i and the Kalman filter state vector The product of:

B.投影矩阵的系统阶数为公式中奇异值矩阵中不为零的奇异值数；B. The system order of the projection matrix is the number of singular values that are not zero in the singular value matrix in the formula;

C.可控可观矩阵为：C. The controllable observable matrix is:

${Γ Γ}_{i i} = = {U u}_{r r} {S S}_{r r}^{11 / / 22};;$

D.状态向量的卡尔曼滤波估计为：D. The Kalman filter estimate of the state vector is:

根据Hankel矩阵的第二种分块方法，投影矩阵用下式表达：According to the second block method of Hankel matrix, the projection matrix is expressed by the following formula:

则下一时刻的卡尔曼滤波状态向量为：Then the Kalman filter state vector at the next moment is:

计算投影矩阵的卡尔曼状态向量：Compute the Kalman state vector of the projection matrix:

由可知，卡尔曼滤波状态向量仅用振动响应数据即可得到；Depend on It can be seen that the Kalman filter state vector It can be obtained only by vibration response data;

最后将求得的卡尔曼滤波状态向量代入系统离散状态空间方程：Finally, the obtained Kalman filter state vector is substituted into the discrete state space equation of the system:

$(\begin{matrix} {\overset{^^}{X x}}_{i i + + 11} \\ {Y Y}_{i i | | i i} \end{matrix}) = = (\begin{matrix} A A \\ C C \end{matrix}) (({\overset{^^}{X x}}_{i i})) + + (\begin{matrix} {ρ ρ}_{w w} \\ {ρ ρ}_{v v} \end{matrix}) - - - - - - ((1111))$

式中，ρ_w，ρ_v为卡尔曼滤波最佳估计的残差；Y_i|i是只有一个块行的Hankel矩阵。In the formula, ρ _w , ρ _v are the residuals of the best estimate of the Kalman filter; Y _i|i is a Hankel matrix with only one block row.

由于卡尔曼状态向量和振动响应数据已知，ρ_ω，ρ_v和相互独立，因此根据最小二乘求解状态方程，得到线性方程：Since the Kalman state vector and vibration response data are known, ρ _ω , ρ _v and are independent of each other, so solving the equation of state according to least squares yields the linear equation:

由此方程可求解出系统状态矩阵A和输出矩阵CFrom this equation, the system state matrix A and output matrix C can be solved

对系统状态矩阵A进行特征值分解：Eigenvalue decomposition of the system state matrix A:

A＝ψΛψ^-1 (13)A＝ψΛψ ^-1 (13)

式中：Λ＝diag(λ_i)，是离散空间系统状态矩阵的复特征值的对角阵，Ψ是特征向量矩阵；In the formula: Λ=diag(λ _i ), is the diagonal matrix of the complex eigenvalues of the state matrix of the discrete space system, and Ψ is the eigenvector matrix;

根据随机子空间理论输电塔系统第i阶固有频率：According to the stochastic subspace theory, the i-th order natural frequency of the transmission tower system is:

${ω ω}_{i i} = = \sqrt{ln ln (({λ λ}_{i i})) ln ln ((\overset{&OverBar; &OverBar;}{{λ λ}_{i i}}))} - - - - - - ((1616))$

系统第i阶阻尼比：The i-th order damping ratio of the system:

${ξ ξ}_{i i} = = \frac{ln ln (({λ λ}_{i i} \overset{&OverBar; &OverBar;}{{λ λ}_{i i}}))}{22} - - - - - - ((1717))$

由离散空间输出矩阵C和离散空间系统状态矩阵的特征向量得到系统第i阶振型为：From the output matrix C of the discrete space and the eigenvector of the state matrix of the discrete space system, the mode shape of the i-th order of the system is obtained as:

4、建立稳态图4. Establish a steady state diagram

在子空间识别结构模态参数的方法中，确定系统的阶次是该方法的关键。奇异值分解确定系统阶次，得到的结果不是很理想。稳定图方法是一种比较新颖的确定系统阶次的方法，同时可以从诸多模态中鉴别真假模态。可用于有较强噪声情形下的模态识别，尤其适合环境激励下的模态定阶。那些满足条件的稳定极点理论上被认为是系统的真实极点。In the method of subspace identification of structural modal parameters, determining the order of the system is the key to the method. Singular value decomposition determines the order of the system, and the results obtained are not very ideal. The stability diagram method is a relatively novel method to determine the order of the system, and can distinguish true and false modes from many modes at the same time. It can be used for modal identification in the case of strong noise, especially for modal order determination under environmental excitation. Those stable poles satisfying the conditions are theoretically considered to be the real poles of the system.

设所识别系统的离散状态空间方程阶次为n，得到多个不同阶次的离散状态空间方程，对每个模型进行模态参数识别，将得到的所有模态参数绘制在二维坐标图中，坐标图的横坐标为频率值，纵坐标为离散状态方程的阶次，便可得到稳定图。稳定图把不同阶次状态方程的模态参数画在同一幅图上，在相应于系统某阶固有频率的轴上，高一阶状态方程识别的模态参数同低一阶状态方程识别的模态参数相比较，如果固有频率、阻尼比和模态振型的差异小于预设的限定值，则这个点就称为稳定点，稳定点组成的轴为稳定轴，相应的模态参数即为系统的模态参数。限定值一般取固有频率＜1％，阻尼比＜5％，模态振型＜2％，即稳定轴要求为：Assuming that the order of the discrete state space equation of the identified system is n, multiple discrete state space equations of different orders are obtained, and the modal parameters of each model are identified, and all the obtained modal parameters are drawn in a two-dimensional coordinate diagram , the abscissa of the coordinate diagram is the frequency value, and the ordinate is the order of the discrete state equation, and the stability diagram can be obtained. The stability diagram draws the modal parameters of different order state equations on the same graph. On the axis corresponding to a certain order natural frequency of the system, the modal parameters identified by the higher order state equation are the same as the modal parameters identified by the lower order state equation. In comparison, if the difference between the natural frequency, damping ratio and mode shape is less than the preset limit value, this point is called a stable point, the axis formed by the stable point is the stable axis, and the corresponding modal parameters are the system Modal parameters. The limit value is generally taken as natural frequency < 1%, damping ratio < 5%, mode shape < 2%, that is, the stable axis requirement is:

$\frac{| | {ξ ξ}^{((n no))} - - {ξ ξ}^{((n no + + 11))} | |}{{ξ ξ}^{((n no))}} \times \times 100100 % % < < 55 % % - - - - - - ((1919))$

(1-MAC(n，n+1))×100％＜2％(1-MAC(n,n+1))×100%＜2%

5、构造Hankel矩阵的方法还可以进行改进：5. The method of constructing Hankel matrix can also be improved:

对于时域算法，各测点的振动数据同步采集测量，对于识别精度的提高效果很好，但由于受到测试条件和手段的限制，对于大型输电塔结构很难同步测量所有测试点振动数据。选择好固定的参考点位置后，分批次逐步完成整个结构的测试是当前采用的主要方法。但分批次测量后存在数据整合的问题，难以避免整合过程带来的拟合误差。为提高振型的拟合识别精度，对现有算法进行改良。For the time-domain algorithm, the vibration data of each measuring point is collected and measured synchronously, which is very effective in improving the recognition accuracy. However, due to the limitation of test conditions and methods, it is difficult to measure the vibration data of all test points synchronously for large transmission tower structures. After selecting a fixed reference point position, the main method currently used is to gradually complete the test of the entire structure in batches. However, there is a problem of data integration after batch measurement, and it is difficult to avoid the fitting error caused by the integration process. In order to improve the fitting recognition accuracy of the mode shape, the existing algorithm is improved.

从采集的数据入手，设参考点为j，响应点数为i，两点之间的互相关函数的表达式为：Starting from the collected data, set the reference point as j, the number of response points as i, and the expression of the cross-correlation function between the two points is:

${R R}_{ij ij} = = {Σ Σ}_{r r = = 11}^{N N} {A A}_{ijr ijr} exp exp ((- - {ξ ξ}_{r r} {ω ω}_{nr nr} t t)) sin sin (({ω ω}_{dr dr} t t + + θ θ)) - - - - - - ((21 twenty one))$

它同脉冲响应函数具有相似的表达式，二者都能表示为衰减正弦函数的和，并且每个衰减正弦都有一个自然频率和阻尼比同结构的各阶模态相对应。由此可见，上式右边包含了系统的全部模态参数，这就确定了响应信号的互相关函数与模态参数之间的关系。It has a similar expression to the impulse response function, both of which can be expressed as the sum of attenuated sine functions, and each attenuated sine has a natural frequency and damping ratio corresponding to each mode of the structure. It can be seen that the right side of the above formula contains all the modal parameters of the system, which determines the relationship between the cross-correlation function of the response signal and the modal parameters.

为了使得不同时间采集得到的各组实验数据同步，对互功率谱进行反傅里叶变换，得到互相关函数：In order to synchronize the experimental data collected at different times, the cross-power spectrum is inversely Fourier transformed to obtain the cross-correlation function:

R_ij＝ifft(F_ij) (22)R _ij ＝ifft(F _ij ) (22)

F_ij为互功率谱，由响应点i和参考点j的自功率谱计算得到。F _ij is the cross power spectrum, which is calculated from the autopower spectrum of the response point i and the reference point j.

用互相关函数代替响应信号，重新建立Hankel矩阵：Re-establish the Hankel matrix by replacing the response signal with a cross-correlation function:

式中，r_i为互相关函数矩阵中的第i列，与i时刻的响应信号相对应。In the formula, _ri is the i-th column in the cross-correlation function matrix, which corresponds to the response signal at time i.

利用新的Hankel矩阵，进行分解和计算，这样就得到了分步测量，整体识别的输电塔模态参数随机子空间识别算法。The new Hankel matrix is used to decompose and calculate, so that the step-by-step measurement and overall identification of the transmission tower modal parameter random subspace identification algorithm is obtained.

以上实施例仅为本发明其中的一种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。The above embodiment is only one implementation mode of the present invention, and its description is relatively specific and detailed, but it should not be construed as limiting the patent scope of the present invention. It should be pointed out that those skilled in the art can make several modifications and improvements without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent for the present invention should be based on the appended claims.

Claims

1. The transmission tower modal parameter identification method based on the improved subspace algorithm, is characterized in that: it comprises the following steps:

1) Set a reference point and multiple measurement points on the transmission tower, collect the vibration response signals of the reference point and measurement points under environmental excitation step by step, and calculate their self-power spectrum and cross-power spectrum;

2) Obtaining a cross-correlation sequence from the cross-power spectrum to form a Hankel matrix;

3) Calculate the projection matrix of the Hankel matrix;

4) weighting the projection matrix;

5) carrying out singular value decomposition to the weighted projection matrix in step 4);

6) Calculate the Kalman filter state vector of the discrete state space equation of the transmission tower system according to the stochastic subspace theory;

7) Solve the system state matrix and output matrix according to the discrete state space equation of the transmission tower system;

8) Solve the system modal parameters according to the system state matrix and output matrix;

9) repeat steps 7) and 8), obtain the modal parameters of the transmission tower system of different orders, and establish the stability diagram of the system modal parameters;

The hankel matrix formed in the step 2) is:

H h = = {Y Y}_{00 | | 22 i i - - 11} = = \frac{11}{\sqrt{j j}} \overset{\overset{j j}{&LeftRightArrow; &LeftRightArrow;}}{[[\frac{\begin{matrix} {y the y}_{00} & {y the y}_{11} & {y the y}_{22} & . . . . . . & {y the y}_{j j - - 11} \\ {y the y}_{11} & {y the y}_{22} & {y the y}_{33} & . . . . . . & {y the y}_{j j} \\ . . . . . . & . . . . . . & . . . . . . & . . . . . . & . . . . . . \\ {y the y}_{i i - - 11} & {y the y}_{i i} & {y the y}_{i i + + 11} & . . . . . . & {y the y}_{j j + + i i - - 22} \end{matrix}}{\begin{matrix} {y the y}_{i i} & {y the y}_{i i + + 11} & {y the y}_{i i + + 11} & . . . . . . & {y the y}_{j j + + i i - - 11} \\ {y the y}_{i i + + 11} & {y the y}_{i i + + 11} & {y the y}_{i i + + 22} & . . . . . . & {y the y}_{j j + + i i} \\ . . . . . . & . . . . . . & . . . . . . & . . . . . . & . . . . . . \\ {y the y}_{22 i i - - 11} & {y the y}_{22 i i} & {y the y}_{22 i i} & . . . . . . & {y the y}_{j j + + 22 i i - - 22} \end{matrix}}]]} = = [\begin{matrix} \frac{{Y Y}_{00 | | i i - - 11}}{{Y Y}_{i i - - 11 | | 22 i i - - 11}} \end{matrix}] = = [\begin{matrix} \frac{{Y Y}_{p p}}{{Y Y}_{f f}} \end{matrix}]

The Hankel matrix is a 2mi×j matrix composed of response data y _k of measurement points, y _i represents the response of all measurement points at time i, subscript p represents the past, subscript f represents the future, Y _p =Y _{0|i- 1} is the block of the Hankel matrix composed of all rows and columns corresponding to the subscript in the first column of the Hankel matrix whose subscript starts from 0 and ends at i-1, indicating the past row space, Y _f =Y _{i-1| 2i-1} is the back i row of Hankel matrix, represents the future row space; Set reference point as j in the described step 2), the number of response points is i, and the expression of the cross-correlation function between two points is:

{R R}_{ij ij} = = {Σ Σ}_{r r = = 11}^{N N} {A A}_{ijr ijr} exp exp ((- - {ξ ξ}_{r r} {ω ω}_{nr nr} t t)) sin sin (({ω ω}_{dr dr} t t + + θ θ));;

A _cmr is the amplitude, ξ _r is the damping ratio, ω _nr is the natural frequency, ω _dr is the resonance frequency, r is the modal order, and N is the total modal order;

Use this cross-correlation function to replace the vibration response signal, perform inverse Fourier transform on the cross-power spectrum calculated in step 1), and obtain the cross-correlation function: R _ij =ifft(F _ij ), where F _ij is the cross-power spectrum , calculated from the autopower spectrum of response point i and reference point j;

Build the Hankel matrix:

H h = = {R R}_{00 | | 22 i i - - 11} = = \frac{11}{\sqrt{j j}} \overset{\overset{j j}{&LeftRightArrow; &LeftRightArrow;}}{[[\frac{\begin{matrix} {r r}_{00} & {r r}_{11} & {r r}_{22} & . . . . . . & {r r}_{j j - - 11} \\ {r r}_{11} & {r r}_{22} & {r r}_{33} & . . . . . . & {r r}_{j j} \\ . . . . . . & . . . . . . & . . . . . . & . . . . . . & . . . . . . \\ {r r}_{i i - - 11} & {r r}_{i i} & {r r}_{i i + + 11} & . . . . . . & {r r}_{j j + + i i - - 22} \end{matrix}}{\begin{matrix} {r r}_{i i} & {r r}_{i i + + 11} & {r r}_{i i + + 11} & . . . . . . & {r r}_{j j + + i i - - 11} \\ {r r}_{i i + + 11} & {r r}_{i i + + 11} & {r r}_{i i + + 22} & . . . . . . & {r r}_{j j + + i i} \\ . . . . . . & . . . . . . & . . . . . . & . . . . . . & . . . . . . \\ {r r}_{22 i i - - 11} & {r r}_{22 i i} & {r r}_{22 i i} & . . . . . . & {r r}_{j j + + 22 i i - - 22} \end{matrix}}]]} = = [\begin{matrix} \frac{{R R}_{00 | | i i - - 11}}{{R R}_{i i - - 11 | | 22 i i - - 11}} \end{matrix}] = = [\begin{matrix} \frac{{R R}_{p p}}{{R R}_{f f}} \end{matrix}]

In the formula, r _i is the i-th column in the cross-correlation function matrix, corresponding to the response signal at time i, R _p =R _0|i-1 is the first i row of the Hankel matrix, representing the past row space, R _f = R _i-1|2i-1 is the last i row of the Hankel matrix, representing the future row space.

2. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 1, is characterized in that: the step of solving the projection matrix of Hankel matrix in described step 3) is:

3.1) Find the projection O _i from the future row space to the past row space of the vibration response signal, and the solution formula is:

{Q Q}_{i i} = = {Y Y}_{f f} / / {Y Y}_{p p} = = {Y Y}_{f f} {Y Y}_{p p}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{+ +} {Y Y}_{p p},,

In the formula, ( ) ⁺ means seeking generalized inversion of matrix, ( ) ^T means seeking transpose of matrix;

3.2) According to the stochastic subspace theory, the discrete state space equation of the identified transmission tower system is:

x _k+1 ＝Ax _k +ρ _w

y _k ＝Cx _k +ρ _v

In the formula, A and C represent the n×n order system state matrix and m×n order output matrix respectively, the order of the discrete state space equation is n, ρ _w , ρ _v are process noise and measurement noise respectively; x _k is k The system state vector at time, y _k is the vibration response vector of the system at time k; the system modal parameters are represented by the eigenvalues and eigenvectors of the system state matrix A, and the system output matrix C;

3.3) According to the random subspace identification theory, the projection matrix O _i can be decomposed into a controllable observable matrix Γ _i and a Kalman filter state vector The product of , that is:

{O o}_{i i} = = {Γ Γ}_{i i} {\overset{^^}{X x}}_{i i} [\begin{matrix} C C \\ CA CA \\ {CA CA}^{22} \\ . . . . . . \\ {CA CA}^{i i - - 11} \end{matrix}] [\begin{matrix} {\overset{^^}{x x}}_{11} & {\overset{^^}{x x}}_{i i + + 11} & {\overset{^^}{x x}}_{i i + + 22} & {\overset{^^}{x x}}_{i i + + j j - - 11} \end{matrix}]

In the formula, is the Kalman filter state vector, which is the optimal estimated value of the system state vector x _k ; Γ _i is a controllable and observable matrix composed of the system state matrix A and the output matrix C.

3. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 2, is characterized in that: adopt principal component PC method to carry out weighting process to projection matrix in step 4), its method is: at first, set The weighting matrix is:

W ₁ =I

{W W}_{22} = = {Y Y}_{P P}^{T T} {(({Y Y}_{p p} {Y Y}_{P P}^{T T}))}^{11 / / 22} {Y Y}_{p p}

In the formula, I is the identity matrix;

Substituting into the projection matrix O _i , the weighted projection matrix obtained is:

o_{i}^{'} = W_{1} o_{i} W_{2} = Y_{f} Y_{P}^{T} {(Y_{p} Y_{P}^{T})}^{- 1 / 2} Y_{p} .

4. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 3, is characterized in that: step 5) carries out the step of singular value decomposition to the projection matrix after weighting process in step 4) is:

5.1) Perform singular decomposition on the projection matrix to determine the system order of the projection matrix:

{O o}_{i i}^{' '} {= = W W}_{11} {O o}_{i i} {W W}_{22} = = UΣV UΣV = = [\begin{matrix} {U u}_{r r} & {U u}_{s the s} \end{matrix}] [\begin{matrix} {S S}_{r r} \\ {S S}_{s the s} = = 00 \end{matrix}] [\begin{matrix} {V V}_{r r} \\ {V V}_{s the s} \end{matrix}] = = {U u}_{r r}^{T T} {S S}_{r r} {V V}_{r r}^{T T},,

In the formula, U _r , S _r , V _r are the left singular value vector, singular value and right singular value vector of the main component; U _s , S _s , V _s are the left singular value vector, singular value and singular value of the noise component Value vector; U, V are orthogonal matrices; S is a diagonal matrix composed of singular values arranged from large to small;

According to the singular value decomposition of the projection matrix, the following results are obtained:

1) The weighted projection matrix O′ _i is equal to the controllable observable matrix Γ _i and the Kalman filter state vector The product of:

2) The system order of the projection matrix is the number of singular values that are not zero in the singular value matrix in the formula;

3) The controllable and observable matrix is:

4) The state vector of the Kalman filter system at time i is:

{\overset{^^}{X x}}_{i i} = = {S S}_{r r}^{11 / / 22} {V V}_{11}^{T T} = = {Γ Γ}_{i i}^{+ +} {O o}_{i i}^{' '};;

5.2) according to the second kind of blocking method of Hankel matrix described in claim 2, projection matrix is expressed with following formula:

{O o}_{i i - - 11}^{' '} = = {Y Y}_{f f}^{- -} / / {Y Y}_{p p}^{+ +} = = {Γ Γ}_{i i - - 11} \cdot \cdot {\overset{^^}{X x}}_{i i + + 11};;

Then the Kalman filter state vector at moment i+1 is:

Γ _i-1 is the matrix obtained by removing the last row element of Γ _i .

5. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 4, is characterized in that: by It can be seen that the Kalman filter state vector It can be obtained only by vibration response data;

Finally, substitute the obtained Kalman filter state vector into the discrete state space equation of the transmission tower system to get:

(\begin{matrix} {\hat{x}}_{i + 1} \\ Y_{i | i} \end{matrix}) = (\begin{matrix} A \\ C \end{matrix}) ({\hat{x}}_{i}) + (\begin{matrix} ρ_{w} \\ ρ_{v} \end{matrix}),

In the formula, ρ _w , ρ _v are the residuals of the best estimate of the Kalman filter; Y _i|i represents all the rows and columns corresponding to the elements whose subscript starts at i and ends at i in the first column of the Hankel matrix The blocks that make up the Hankel matrix represent a Hankel matrix with only one block row.

6. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 5, is characterized in that: because Kalman filtering state vector and vibration response data are known, ρ _ω , ρ _v and Independently of each other, the discrete state space equation of the transmission tower system is solved by using the least square method widely used in the field of numerical calculation, and the linear equation is obtained: From this equation, the system state matrix A and output matrix C are solved.

7. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 6, is characterized in that: described step 6) also comprises the following steps:

6.1) Decompose the eigenvalue of the system state matrix A:

A＝ψΛψ- ¹

In the formula: Λ=diag(λ _i ), is the diagonal matrix of the complex eigenvalues of the system state matrix, and ψ is the eigenvector matrix;

6.2) Calculate the natural frequency and damping ratio of the system according to the eigenvalues of the system state matrix:

Among them, the complex eigenvalue λ _i of the system state matrix, is the conjugate complex number of λ _i ;

According to the stochastic subspace theory, for the transmission tower system with unit sampling interval time, the i-th order natural frequency of the system is

ω_{i} = \sqrt{\ln (λ_{i}) \ln (\overset{&OverBar;}{λ_{i}})},

The i-th order damping ratio of the system

8.3) From the output matrix C and the eigenvector of the system state matrix, the i-th order mode shape of the system is obtained as:

8. the transmission tower modal parameter identification method based on improved subspace algorithm according to claim 7, is characterized in that: step 9) in setting up stable figure and identifying system stability, the method for correct modal parameter is:

Assuming that the order of the discrete state space equation of the identified transmission tower system is n, multiple discrete state space equations of different orders are obtained, and the modal parameters of each equation are identified, and all the obtained modal parameters are drawn on the two-dimensional coordinates In the figure, the abscissa of the coordinate diagram is the frequency value, and the ordinate is the order of the discrete state space equation of the transmission tower system, and the stability diagram can be obtained;

The stability diagram draws the modal parameters of different order equations on the same graph. On the axis corresponding to a certain order natural frequency of the system, the modal parameters identified by the higher-order state-space equations are the same as those identified by the lower-order state-space equations. Compared with parameters, if the difference of natural frequency, damping ratio and mode shape is less than the preset limit value, then this point is called stable point, the axis composed of stable points is stable axis, and the corresponding modal parameters are the system modal parameters; the limit values are natural frequency < 1%, damping ratio < 5%, mode shape < 2%, that is, the stable axis requirements are:

\frac{| | {ω ω}^{((n no))} - - {ω ω}^{((n no + + 11))} | |}{{f f}^{((n no))}} \times \times 100100 % % < < 11 % %

\frac{| | {ξ ξ}^{((n no))} - - {ξ ξ}^{((n no + + 11))} | |}{{ξ ξ}^{((n no))}} \times \times 100100 % % < < 55 % %

(1-MAC(n,n+1))×100%＜2%

In the formula, n represents the order of the discrete state space equation of the transmission tower system, and the modal parameters of the system calculated by the state space equation of order n are: natural frequency ω ⁽ⁿ⁾ , damping ratio ξ ⁽ⁿ⁾ , modal vibration type MAC represents the Modality Assurance Criterion: