Underdetermined operating modal parameter identification method and system based on tensor decomposition recursion

Abstract

The present invention discloses a method and system for identifying underdetermined working modal parameters based on tensor decomposition recursion, which relates to the technical field of working modal parameter identification. The method comprises: obtaining a linear time-varying structural vibration response signal of a sensor measuring point under a set environmental excitation, and extracting the principal component score of the response signal through PCA; determining an appropriate lag period by calculating the autocorrelation of the principal component score and combining the peak value in the autocorrelation function, and constructing it and the response signal into a third-order tensor; under time-varying conditions, establishing a weighted observation matrix through recursive decomposition and dynamic forgetting factor selection mechanism, recursively updating the decomposed results, thereby obtaining more accurate modal vibration shapes and natural frequencies. The present invention significantly improves the recognition accuracy, stability and robustness under complex time-varying working modes, and is suitable for large-scale civil engineering structural health monitoring, such as bridges, buildings and wind power towers, etc., to help timely discover structural damage.

Description

Technical Field

[0001] The present invention relates to the technical field of operating modal parameter identification, and in particular to an underdetermined operating modal parameter identification method and system based on tensor decomposition recursion. Background Art

[0002] Operational modal parameter identification (OMPI) is a blind identification method that identifies system modal parameters by relying solely on measured structural responses. Compared to traditional experimental modal parameter identification, OMPI offers the following advantages: It only requires obtaining time-domain vibration output data, eliminating the need for input load measurement; and because the measured data is derived directly from the structure's actual operating state, the identification results are more consistent with realistic boundary conditions and operating conditions. Therefore, OMPI is a powerful tool, particularly suitable for analyzing the dynamic characteristics of large, complex structures. Because it uses data from actual operating conditions for modal parameter identification, the identification results are more accurate and closer to reality. More importantly, OMPI also enables online modal analysis of structures, allowing for real-time monitoring of their dynamic performance. Due to its practicality and ability to be performed online, OMPI has been widely used in analyzing large-scale engineering structures, such as aircraft, automobiles, and civil engineering.

[0003] With the emergence of more large and high-speed engineering structures, the dynamic challenges of linear systems involving time-varying mass, stiffness, and damping are increasing. Therefore, accurately predicting the vibration response of structures has become increasingly critical.

[0004] Although the theory of operating modal parameter identification has matured, it is primarily based on the assumption that the system is linearly time-varying, which can lead to limitations in certain situations. The identification of operating modal parameters for time-varying systems is still in its early stages of development. In the real world, many engineering structures exhibit time-varying characteristics due to the time-varying structural parameters (such as mass, stiffness, and damping) of their operational parameters. For example, in a model of a train crossing a bridge, the mass of the train itself causes the dynamic characteristics of the train-bridge coupling system to continuously change as the train moves. During the launch of a rocket or missile, changes in system mass caused by fuel consumption also affect its dynamic characteristics. Furthermore, time-varying structural characteristics are increasingly common in fields such as energy and shipbuilding. For example, the rotation of a wind turbine's rotor blades during operation can cause changes in the dynamic characteristics of the entire structural system. Current industrial development has placed new demands on structural dynamics analysis, and therefore, the study of parameter identification for time-varying structures has become increasingly urgent. Summary of the Invention

[0005] In response to the problems of the prior art, the purpose of the present invention is to provide a method and system for identifying underdetermined working modal parameters based on tensor decomposition recursion. The obtained signal matrix is ​​optimized and constructed into a third-order tensor through the PCA (autocorrelation optimized principal component) method, and the constructed tensor is simultaneously diagonalized and recursively decomposed. The decomposition results are analyzed to obtain the corresponding modal vibration shapes and natural frequencies; without destroying the original connection of the signal, the modal vibration shapes and natural frequencies under time-varying conditions are better identified, thereby improving the recognition accuracy and stability.

[0006] The present invention adopts the following technical solutions:

[0007] On the one hand, a method for identifying underdetermined operating modal parameters based on tensor decomposition recursion includes:

[0008] The vibration response signal acquisition step is to acquire the vibration response signal X(t) of the time-varying structure collected at the selected sensor measurement point at time t under the set environmental excitation; wherein, n m is the number of sensors, N(t) is the number of time samples;

[0009] The principal component score calculation step is to process the obtained vibration response signal using principal component analysis (PCA) to calculate the principal component score;

[0010] The lag period determination step is to calculate the autocorrelation function of the principal component scores and determine the lag period τ by detecting the peak in the autocorrelation function. k ;

[0011] The third-order tensor construction step is based on the modal vibration matrix Φ(t), modal response matrix Q(t) and lag period τ obtained by decomposing the vibration response signal X(t) at time t k Constructed as a dimension of n m ×N(t)×τ k The third-order tensor of H (t);

[0012] The third-order tensor update step, the third-order tensor at time t H Add new data slices in the second dimension of (t) to construct a third-order tensor at time t+1 H (t+1);

[0013] The weighted observation matrix construction step is based on the third-order tensor at time t+1 H (t+1), a dynamic forgetting factor is introduced to construct the weighted observation matrix H(t+1) at time t+1, and the observation matrix is ​​windowed and weighted;

[0014] The factor matrix recursive update step is based on the two factorizations of the weighted observation matrix H(t) at time t, the diagonalization of the connection matrix, the tracking of singular values, and the common blocks between the modal response matrices Q(t) and Q(t+1). The modal vibration shape matrix Φ(t+1), the modal response matrix Q(t+1), and the signal to modal contribution matrix Θ(t+1) at time t+1 are obtained.

[0015] In the working modal parameter identification step, the modal vibration matrix Φ(t+1), the modal response matrix Q(t+1) and the signal to modal contribution matrix Θ(t+1) are decomposed to identify the time-varying working modal parameters.

[0016] Preferably, the principal component score calculation step specifically includes:

[0017] The collected vibration response signal X(t) is normalized to obtain the normalized matrix B(t), as follows:

[0018] B(t)=(X(t)-μ) / σ;

[0019] Among them, μ is the mean of the original data, σ is the standard deviation;

[0020] The covariance matrix Z(t) is obtained by calculating the standardized matrix B(t), and the covariance matrix Z(t) is diagonalized to obtain

[0021] The matrix V(t) composed of eigenvectors is as follows:

[0022] V(t)Λ(t)V T (t) = Z(t);

[0023] Where T is the transpose of the matrix, Λ(t) is a diagonal matrix,

[0024] The principal component score matrix S(t) is obtained based on the normalized matrix B(t) and the matrix V(t) composed of eigenvectors, as follows:

[0025] S(t)=B(t)V(t).

[0026] Preferably, the hysteresis period determination step specifically includes:

[0027] Calculate the autocorrelation function R(l) of the principal component score as follows:

[0028]

[0029] Among them, i is the index of a selected principal component, that is, the i-th column of the principal component score matrix S(t), S t，irepresents the value of the i-th principal component score at time t; S t+l，i represents the value of the i-th principal component score at time t+1; represents the mean of the i-th principal component score, that is, N is the length of the data, l is the number of lags;

[0030] Using the autocorrelation function R(l), the peak detection is used to determine the periodic characteristics and determine the appropriate lag period τ k :

[0031] locs={k|R(l)>R(l-1)and R(l)>R(l+1)};

[0032]

[0033] Where k represents the peak index; len(R(l)) represents the length of the autocorrelation function R(l), that is, the number of autocorrelation function values; locs represents the peak array; and len(locs) represents the number of peaks.

[0034] Preferably, the modal shape matrix Φ(t), modal response matrix Q(t) and lag period τ obtained by decomposing the vibration response signal X(t) at time t k Constructed as a dimension of n m ×N(t)×τ k The third-order tensor of H (t), is expressed as follows:

[0035]

[0036] Where R is the rank of CP decomposition; Φ r (t) is a tensor H (t) is the vector in the first dimension, i.e., the rth modal vibration shape vector matrix of the modal vibration shape matrix Φ(t); Q r (t) is a tensor H (t) is the vector in the second dimension, i.e., the rth modal response vector matrix of the modal response matrix Q(t); Θ r (t) is a tensor H (t) is a vector in the third dimension, i.e., the rth vector matrix of the signal contribution matrix Θ(t); operator Represents the vector outer product.

[0037] Preferably, the third-order tensor at time t is H Add new data slices in the second dimension of (t) to construct a third-order tensor at time t+1 H (t+1), expressed as follows:

[0038] H(t+1)=[ H (t), h(t+1)];

[0039] Wherein, h(t+1) represents the new data slice.

[0040] Preferably, the weighted observation matrix construction step specifically includes:

[0041] Calculate the data change rate change_rate of the vibration response signal input at the current time point, and dynamically adjust the forgetting factor λ(t+1) at time t+1 based on the data change rate and the forgetting factor λ(t) at time t, as follows:

[0042]

[0043] Among them, λ base represents the initial value of the forgetting factor;

[0044] Use the dynamically adjusted forgetting factor λ(t+1) to weight the historical data and construct a new weighted observation matrix H(t+1) as follows:

[0045]

[0046] Among them, H (1) (t+1) is the tensor H (t+1) The matrix after expansion along the first dimension, dimension is n m ×(N(t+1)×τ k ) to extract modal information from the data using matrix decomposition techniques; is a weighted matrix consisting of the current forgetting factor, as follows:

[0047]

[0048] Wherein, diag(λ(t+1),λ(t),...,1) is a diagonal matrix constructed using the forgetting factor λ(t+1);

[0049] The weighted observation matrix H(t+1) is decomposed to obtain the initial modal vibration shape matrix Φ(t+1), the modal response matrix Q(t+1), and the signal to modal contribution matrix Θ(t+1). Let A(t+1) = Φ(t+1)⊙Θ(t+1) to obtain the loading matrix A(t+1) at time t+1; the operator ⊙ represents the Hadamard product.

[0050] Preferably, the factor matrix recursive updating step specifically includes:

[0051] The two factorizations of the weighted observation matrix H(t) at time t are as follows:

[0052]

[0053] The first decomposition is to transform the formula H in (1) (t) is replaced by the result obtained by its CP decomposition, ; The second decomposition is the economic-scale SVD decomposition of H(t), where and U(t) is a combination of singular vectors containing the tensor in the first and third dimensions; ∑(t) contains the singular values ​​in the decomposition of U(t); L(t) contains the singular vectors containing the tensor in the second dimension;

[0054] There exists a non-singular connectivity matrix So that:

[0055]

[0056] in, is the auxiliary matrix; the connection between the two equations W(t) is found by solving the simultaneous diagonalization of a set of matrices of size R×R;

[0057] Combining the above two formulas, we can calculate E(t), L(t), W(t) and W -1 (t); W -1 (t) is the inverse matrix of W(t);

[0058] Perform singular value decomposition on the loading matrix A(t+1), as follows:

[0059] [Θ r (t+1),σ r , Φ r (t+1)]=svd(A r (t+1));

[0060] Among them, A r (t+1) represents the rth vector of A(t+1); Θ r (t+1) is the matrix A r The left singular vector matrix of (t+1) is A. r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the row space; σ r It represents the diagonal matrix of singular values, including A r The singular values ​​of (t+1) reflect the strength of the matrix in each orthogonal direction; Φ r (t+1) is the matrix A r The right singular vector matrix of (t+1) is A. r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the column space; the updated factor matrix Θr (t+1) and Φ r (t+1) is used to adjust the tensor decomposition results in the recursive process;

[0061] Based on the singular value decomposition results, we obtain the updated L(t+1) of L(t) at time t+1 and the updated E(t+1) of E(t) at time t+1 as follows:

[0062] Assign the left singular vector matrix to U(t+1) to update U(t), and assign the decomposed singular value to ∑(t+1) to update ∑(t); based on Get the updated E(t+1) of E(t) at time t+1;

[0063] Assign part of the right singular vector matrix to L(t+1) to update L(t);

[0064] At time t+1, the formula The conversion is as follows:

[0065]

[0066] Using the time-shift structure of Q(t+1), time t and time t+1 are connected, where E(t+1) represents the auxiliary matrix at time t+1; W(t+1) represents the non-singular connection matrix at time t+1; W -1 (t+1) is the inverse matrix of W(t+1);

[0067] Since the time-varying structures Q(t) and Q(t+1) have a common block, we can write The dimensions of L(t) and L(t+1) are N(t)×R and N(t+1)×R respectively, so the dimensions of Q(t) and Q(t+1) are also the same; in this case, the common block is Q(t) itself;

[0068] Define the following matrix based on the specified window:

[0069]

[0070] Indicates that L(t) at time t does not change;

[0071] It means to intercept the first N(t) rows from the L(t+1) matrix to form

[0072] It means extracting the N(t+1)th row of L(t+1), forming

[0073] By public block get:

[0074]

[0075] Where λ(t) 1 / 2 represents the square root of the forgetting factor, which is used to weight historical data;

[0076] get:

[0077]

[0078] in, represents the pseudo-inverse of a matrix; is an identity matrix, so

[0079] Update Φ(t) and Θ(t): Based on the formula Φ(t+1)⊙Θ(t+1)=E(t+1)W(t+1), update Φ(t) and Θ(t) to obtain Φ(t+1) and Θ(t+1);

[0080] Update of Q(t): Based on the formula Update qT(t+1), which is eventually appended to Q(t) to construct Q(t+1);

[0081] Determine whether the change of Q(t+1) relative to Q(t) is less than a preset second threshold. If so, further determine whether the factor matrices of all ranks have been updated. If so, output the updated Φ(t+1), Q(t+1) and Θ(t+1).

[0082] In another aspect, a system for identifying under-timed variable operating modal parameters includes:

[0083] A motor system and a support platform for providing mechanical vibration and structural support;

[0084] Measurement and analysis system, including force and motion sensors, used to acquire the vibration response signals of the time-varying structure;

[0085] An excitation system, used to apply control excitations on the time-varying structure to simulate the real working environment;

[0086] The processing unit analyzes the collected vibration response signal based on the underdetermined working modal parameter identification method based on tensor decomposition recursion and identifies the time-varying working modal parameters.

[0087] Compared with the prior art, the present invention has the following beneficial effects:

[0088] The present invention obtains the linear time-varying structural vibration response signal of the sensor measuring point and constructs it into a tensor. The adaptive tensor decomposition technology can more accurately capture the modal parameters in the time-varying structure, thereby significantly improving the recognition accuracy; by combining the simultaneous diagonalization tracking technology and the dynamic selection of forgetting factor technology, it can maintain high stability and robustness when facing complex working modes; the method of the present invention is particularly suitable for use in the health monitoring of large-scale civil engineering structures (such as bridges, buildings, wind power towers, etc.), and similar neighborhoods where the modal parameters of the structure (such as frequency, damping ratio and modal shape) change with time and external environment help monitor the health status of the structure and promptly discover potential structural damage. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 Flowchart of a method for identifying underdetermined operating modal parameters based on tensor decomposition recursion according to an embodiment of the present invention;

[0090] Figure 2 A diagram of a time-varying three-degree-of-freedom spring oscillator model provided by an embodiment of the present invention;

[0091] FIG3( a ) is a diagram of a Gaussian white noise excitation simulation signal according to an embodiment of the present invention;

[0092] FIG3( b ) is a simulation signal diagram of a displacement response signal 1 ( m1 ) provided by an embodiment of the present invention;

[0093] FIG3( c ) is a simulation signal diagram of a displacement response signal 2 ( m2 ) provided by an embodiment of the present invention;

[0094] FIG3( d ) is a simulation signal diagram of a displacement response signal 3 ( m3 ) provided by an embodiment of the present invention;

[0095] Figure 4 Comparison between the experimental method provided in the embodiment of the present invention for identifying the theoretical natural frequency of a three-degree-of-freedom time-varying structure and the identification of a natural frequency variation curve;

[0096] Figure 5 The experimental method provided in the embodiment of the present invention identifies the MAC value change curve of a three-degree-of-freedom time-varying system;

[0097] FIG6( a ) is a comparison of the natural vibration mode of a three-degree-of-freedom time-varying system at 38.125 s using the experimental method provided by an embodiment of the present invention and the instantaneous vibration mode identified by the proposed method;

[0098] FIG6( b ) is a comparison of the natural vibration mode of the three-degree-of-freedom time-varying system at 450 s using the experimental method provided by an embodiment of the present invention and the instantaneous vibration mode identified by the proposed method;

[0099] FIG6( c ) is a comparison of the natural vibration mode of the three-degree-of-freedom time-varying system in 900 s using the experimental method provided by an embodiment of the present invention and the instantaneous vibration mode identified by the proposed method;

[0100] FIG6( d ) is a comparison of the natural vibration mode of the three-degree-of-freedom time-varying system at 1500 s using the experimental method provided by an embodiment of the present invention and the instantaneous vibration mode identified by the proposed method;

[0101] Figure 7 A schematic structural diagram of an under-timed variable operating modal parameter identification system provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0102] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention for which protection is sought, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0103] like Figure 1 As shown, the present invention provides an underdetermined operating modal parameter identification method based on tensor decomposition recursion, which includes the following steps.

[0104] S1, a vibration response signal acquisition step, obtaining the vibration response signal X(t) of the time-varying structure collected at the selected sensor measurement point at time t under the set environmental excitation; wherein, n m is the number of sensors, and N(t) is the number of time samples.

[0105] Specifically, multiple sensors are arranged on the time-varying structure to be monitored, and these sensors are used to collect the vibration response signals of the structure in real time. Assume that the structure has n m These sensors record vibration signals at a fixed sampling rate to form the vibration response matrix X(t) at time t.

[0106] S2, a principal component score calculation step, uses principal component analysis (PCA) to process the acquired vibration response signal and calculate the principal component score.

[0107] Specifically, the collected vibration response signal X(t) is normalized to obtain a normalized matrix B(t), as follows:

[0108] B(t)=(X(t)-μ) / σ;

[0109] Among them, μ is the mean of the original data, σ is the standard deviation;

[0110] The covariance matrix Z(t) is obtained by calculating the standardized matrix B(t), and the covariance matrix Z(t) is diagonalized to obtain

[0111] The matrix V(t) composed of eigenvectors is as follows:

[0112] V(t)Λ(t)V T (t) = Z(t);

[0113] Where T is the transpose of the matrix, Λ(t) is a diagonal matrix,

[0114] The principal component score matrix S(t) is obtained based on the normalized matrix B(t) and the matrix V(t) composed of eigenvectors, as follows:

[0115] S(t)=B(t)V(t).

[0116] S3, lag period determination step, calculates the autocorrelation function of the principal component score and determines the lag period τ by detecting the peak in the autocorrelation function k .

[0117] Specifically, the autocorrelation function R(l) of the principal component score is calculated as follows:

[0118]

[0119] Among them, i is the index of a selected principal component, that is, the i-th column of the principal component score matrix S(t), S t，i represents the value of the i-th principal component score at time t; S t+l，i represents the value of the i-th principal component score at time t+1; represents the mean of the i-th principal component score, that is, N is the length of the data, l is the number of lags;

[0120] Using the autocorrelation function R(l), the peak detection is used to determine the periodic characteristics and determine the appropriate lag period τ k :

[0121] locs={k|R(l)>R(l-1)and R(l)>R(l+1)};

[0122]

[0123] Where k represents the peak index; len(R(l)) represents the length of the autocorrelation function R(l), that is, the number of autocorrelation function values; locs represents the peak array; and len(locs) represents the number of peaks.

[0124] S4, the third-order tensor construction step, based on the modal vibration matrix Φ(t), modal response matrix Q(t) and lag period τ obtained by decomposing the vibration response signal X(t) at time t k Constructed as a dimension of n m ×N(t)×τ k The third-order tensor of H (t).

[0125] Specifically, the third-order tensor H (t), is expressed as follows:

[0126]

[0127] Where R is the rank of CP decomposition; Φ r (t) is a tensor H (t) is the vector in the first dimension, i.e., the rth modal vibration shape vector matrix of the modal vibration shape matrix Φ(t); Q r (t) is a tensor H (t) is the vector in the second dimension, i.e., the rth modal response vector matrix of the modal response matrix Q(t); Θ r (t) is a tensor H (t) is a vector in the third dimension, that is, the rth vector matrix of the signal contribution matrix Θ(t).

[0128] S5, the third-order tensor update step, the third-order tensor at time t H Add new data slices in the second dimension of (t) to construct a third-order tensor at time t+1 H (t+1).

[0129] Specifically, the third-order tensor at time t+1 is H (t+1), expressed as follows:

[0130] H (t+1)=[ H (t), h(t+1)];

[0131] Wherein, h(t+1) represents the new data slice.

[0132] S6, weighted observation matrix construction step, based on the third-order tensor at time t+1 H (t+1), a dynamic forgetting factor is introduced to construct the weighted observation matrix H(t+1) at time t+1, and the observation matrix is ​​windowed and weighted.

[0133] Specifically, at the initial stage of the algorithm, the initial value of the forgetting factor λ(t) is set to λ base , and use the initial data sample number M to build the initial model.

[0134] Calculate the data change rate change_rate of the vibration response signal input at the current time point, and adjust the forgetting factor λ(t+1) at time t+1 according to the change rate as follows:

[0135]

[0136] Use the dynamically adjusted forgetting factor λ(t+1) to weight the historical data and construct a new weighted observation matrix H(t+1):

[0137]

[0138] Among them, H (1) (t+1) is the tensor H (t+1) The matrix after expansion along the first dimension has a dimension of n m ×(N(t+1)×τ k ), the purpose of this processing is to use matrix decomposition technology to extract modal information from the data. is a weighted matrix consisting of the current forgetting factor:

[0139]

[0140] Wherein, diag(λ(t+1),λ(t),...,1) is a diagonal matrix constructed using the forgetting factor λ(t+1).

[0141] The weighted observation matrix H(t+1) is decomposed to obtain the initial factor matrices Φ(t+1), Q(t+1) and Θ(t+1).

[0142] The factor matrix Q(t+1) is recursively updated, and the initial estimate q of Q(t+1) at time t+1 is estimated by the least squares method. T (t+1):

[0143]

[0144] Among them, A(t+1)=Φ(t+1)⊙Θ(t+1) is the loading matrix at the new moment, is the pseudo-inverse matrix of the loading matrix A(t+1).

[0145] S7, the factor matrix recursive update step, is based on the two factorizations of the weighted observation matrix H(t) at time t, the diagonalization of the connection matrix, the singular value tracking, and the common blocks between the modal response matrices Q(t) and Q(t+1), and recursively updates to obtain the modal vibration shape matrix Φ(t+1), the modal response matrix Q(t+1), and the signal to modal contribution matrix Θ(t+1) at time t+1.

[0146] Specifically, consider the two factorizations of the weighted observation matrix H(t) at time t, as follows:

[0147]

[0148] The first decomposition is to transform the formula H in (1) (t) is replaced by the result obtained by its CP decomposition, ; The second decomposition is the economic-scale SVD decomposition of H(t), where and U(t) is a combination of singular vectors of the tensor in the first and third dimensions; ∑(t) is a combination of singular values ​​in the decomposition of U(t); L(t) is a combination of singular vectors of the tensor in the second dimension.

[0149] There exists a non-singular connectivity matrix So that:

[0150]

[0151] in, is the auxiliary matrix; the connection between the two equations W(t) is found by solving the simultaneous diagonalization of a set of matrices of size R×R.

[0152] Combining the above two formulas, we can calculate E(t), L(t), W(t) and W -1 (t); W -1 (t) is the inverse matrix of W(t).

[0153] Next, SVD tracking is performed. First, both the left and right subspaces (the left subspace is the left singular vectors of the modal shape matrix, and the right subspace is the right singular vectors of the modal response matrix) must be tracked. Furthermore, the orthogonal (orthogonalized) right subspace basis vectors must be calculated. These basis vectors are then used to derive the recursive update of W(t), H(t) with increasing dimensions, using a dual iteration technique to update E(t) and L(t). L(t) represents a matrix obtained from the simultaneous diagonalization recursion, specifically a portion of the right singular vector matrix.

[0154] During specific tracking, the loading matrix A(t+1) is subjected to singular value decomposition as follows:

[0155] [Θ r (t+1),σ r , Φ r (t+1)]=svd(A r (t+1));

[0156] Among them, A r (t+1) represents the rth vector of A(t+1); Θ r (t+1) is the matrix A r The matrix composed of the left singular vectors of (t+1), whose column vector is A r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the row space; σ r It represents the diagonal matrix of singular values, including A r The singular values ​​of (t+1) reflect the strength of the matrix in each orthogonal direction; Φ r (t+1) is the matrix A r The matrix composed of the right singular vectors of (t+1), whose column vector is A r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the column space; the updated factor matrix Θ r (t+1) and Φ r (t+1) is used to adjust the tensor decomposition results in the recursive process; Θ r (t+1) and Φ r (t+1) represents the modal components used for updating in each recursive iteration. During the recursive process, they are iteratively updated by the least squares method to ensure that the model can gradually approach the true modal characteristics. Based on the singular value decomposition results, the updated L(t) at time t+1 L(t+1) and the updated E(t) at time t+1 E(t+1) are obtained as follows:

[0157] Assign the left singular vector matrix to U(t+1) to update U(t), and assign the decomposed singular value to ∑(t+1) to update ∑(t); based on Get the updated E(t+1) of E(t) at time t+1;

[0158] Assign the portion of the right singular vector matrix to L(t+1) to update L(t).

[0159] At time t+1, the formula The conversion is as follows:

[0160]

[0161] Using the time-shift structure of Q(t+1), time t and time t+1 are connected, where E(t+1) represents the auxiliary matrix at time t+1; W(t+1) represents the non-singular connection matrix at time t+1; W -1 (t+1) is the inverse matrix of W(t+1).

[0162] Since the time-varying structures Q(t) and Q(t+1) have a common block, we can write The dimensions of L(t) and L(t+1) are N(t)×R and N(t+1)×R respectively, so the dimensions of Q(t) and Q(t+1) are also the same; in this case, the common block is Q(t) itself.

[0163] Define the following matrix based on the specified window:

[0164]

[0165] Indicates that L(t) at time t does not change;

[0166] It means to intercept the first N(t) rows from the L(t+1) matrix to form

[0167] It means extracting the N(t+1)th row of L(t+1), forming

[0168] By public block get:

[0169]

[0170] Where λ(t) 1 / 2 represents the square root of the forgetting factor, which is used to weight historical data. Specifically, λ(t) 1 / 2 It is based on the attenuation of the forgetting factor λ(t), which makes the weighted decrease smoother.

[0171] get:

[0172]

[0173] in, represents the pseudo-inverse of a matrix; is an identity matrix, so

[0174] Update of Φ(t) and Θ(t): Based on the formula Φ(t+1)⊙Θ(t+1)=E(t+1)W(t+1), update Φ(t) and Θ(t) to obtain Φ(t+1) and Θ(t+1).

[0175] Update of Q(t): Based on the formula Update qT(t+1) and eventually append to Q(t) to construct Q(t+1).

[0176] Determine whether the change of Q(t+1) relative to Q(t) is less than a preset second threshold. If so, further determine whether the factor matrices of all ranks have been updated. If so, output the updated Φ(t+1), Q(t+1) and Θ(t+1).

[0177] S8, the working modal parameter identification step, decomposes the modal vibration matrix Φ(t+1), the modal response matrix Q(t+1) and the signal to modal contribution matrix Θ(t+1) to identify time-varying working modal parameters, such as modal vibration shapes and natural frequencies.

[0178] Through the above embodiments, the under-time-varying working modal parameter identification method based on tensor decomposition with forgotten aging factors and simultaneous diagonalization recursion proposed in the present invention can adaptively and recursively update the Φ(t) modal array matrix, Q(t) modal response matrix and Θ(t) and other weight matrices, thereby effectively processing the identification of time-varying structure working modal parameters.

[0179] Furthermore, in another preferred embodiment of the present invention, a linear time-varying three-degree-of-freedom spring oscillator system is used to verify the modal parameter identification method proposed in the present invention.

[0180] like Figure 2 The linear time-varying three-degree-of-freedom spring oscillator system shown.

[0181] Assume that the sampling frequency f is 40 Hz, the sampling interval Δt is 0.025 s, and the sampling time is t=2000 s.

[0182] Assume that the initial conditions of the system's three degrees of freedom displacement are all zero, apply Gaussian white noise excitation to block 1, and set the mass to:

[0183]

[0184] m2=m3=1kg, stiffness kl=1000N / m, k2=1000N / m, k3=1000N / m; damping c1=0.01Ns / m,

[0185] c2=0.01Ns / m, c3=0.01Ns / m.

[0186] Then, the three-degree-of-freedom spring oscillator was tested. The response signal was simulated using MATLAB / Newmark-β. The sampling frequency of the signal was 40 Hz, the sampling time was t = 2000 s, and the sampling interval was 0.025 s. The Gaussian white noise excitation applied to the object m1 and the simulated response signal are shown in Figure 1. Figure 3(a) 、 3(b) , 3(c) and 3(d).

[0187] like Figure 4 As shown, this embodiment compares and evaluates the working modal parameter identification methods of linear time-varying structures through the Modal Assurance Criterion (MAC) and the relative error of the modal natural frequency.

[0188] The formula for the modal confidence criterion is as follows:

[0189]

[0190] in, and They represent the estimated modal shape vector and the true modal shape vector of the i-th order, respectively, and the MAC value varies between 0 and 1. The larger the MAC value (closer to 1), the better the mode shape.

[0191] The formula for the natural frequency relative error evaluation method is as follows:

[0192]

[0193] Among them, ω i With ω i-theory denote the estimated modal natural frequency and theoretical natural frequency of the i-th order, ζ i The closer the value is to 0, the more accurate the identified natural frequency is.

[0194] The MAC value identification results of the three-degree-of-freedom spring oscillator obtained in this experiment at the instantaneous moment of 900s are shown in Table 1. The natural frequency error is shown in Table 2. The MAC value change curve of the three-degree-of-freedom time-varying system identified by the experimental method provided by the embodiment of the present invention can be found in Figure 5 shown.

[0195] Table 1. Confidence criteria of each modal at 900 s;

[0196] Order MAC 1 1.0000 2 <![CDATA[0 . 9972]]> 3 <![CDATA[0 . 9856]]>

[0197] Table 2. Natural frequency errors of various orders at 900 s;

[0198] Order Theoretical natural frequency Identifying natural frequencies Natural frequency error 1 <![CDATA[2 . 2399]]> <![CDATA[2 . 2600]]> <![CDATA[0 . 8973%]]> 2 <![CDATA[6 . 2760]]> <![CDATA[6 . 8520]]> 9.1778％ 3 9.0690 10.0400 10.7068％

[0199] Figure 6(a) 、 6(b), 6(c) and 6(d) are comparisons of the natural vibration modes of the three-degree-of-freedom time-varying system at different times of 38.125s, 450s, 900s and 1500s using the experimental method provided in an embodiment of the present invention and the instantaneous vibration modes identified by the proposed method.

[0200] From Table 1, Table 2, Figure 6(a) 、 6(b) , 6(c) and 6(d), it can be seen that the method proposed in the present invention can well realize the identification of modal parameters under under-time variation conditions.

[0201] Table 3 lists the global symbols, meanings and specifications in this invention.

[0202] Table 3. Global symbol table;

[0203]

[0204]

[0205] like Figure 7 As shown, an embodiment of the present invention further provides an under-timed variable working modal parameter identification system, comprising: a motor system and a support platform (DC micro-motor motion controller) for providing mechanical vibration and structural support;

[0206] Measurement and analysis system, including force and motion sensors, to collect vibration response signals of time-varying structures;

[0207] Excitation system (vibrator), used to apply controlled excitation to simulate the real working environment;

[0208] The processing unit (computer) analyzes the collected signals and identifies the time-varying modal parameters based on the underdetermined working modal parameter identification method based on tensor decomposition recursion.

[0209] Furthermore, the present invention also provides an underdetermined working modal parameter identification device, which includes a memory and a processor, wherein the memory stores a computer program, and the computer program can be executed by the processor to implement the underdetermined working modal parameter identification method based on tensor decomposition recursion as described above.

[0210] Furthermore, the present invention also provides a computer-readable storage medium, which includes computer-readable instructions stored on the computer-readable storage medium, and when the computer-readable instructions are executed by the processor of the device where the computer-readable storage medium is located, the underdetermined working modal parameter identification method based on tensor decomposition recursion as described above is implemented.

[0211] In the several embodiments provided in the embodiments of the present invention, it should be understood that the disclosed devices and methods can also be implemented in other ways. The device and method embodiments described above are merely illustrative. For example, the flowcharts in the accompanying drawings show the possible architectures, functions, and operations of the devices, methods, and computer program products according to multiple embodiments of the present invention. In this regard, each box in the flowchart or block diagram can represent a module, program segment, or part of a code, which contains one or more executable instructions for implementing the specified logical functions. It should also be noted that in some alternative implementations, the functions marked in the boxes can also occur in an order different from that marked in the drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the block diagram and / or flowchart, as well as the combination of boxes in the block diagram and / or flowchart, can be implemented using a dedicated hardware-based system that performs the specified functions or actions, or can be implemented using a combination of dedicated hardware and computer instructions.

[0212] In addition, the functional modules in the various embodiments of the present invention may be integrated together to form an independent part, or each module may exist independently, or two or more modules may be integrated to form an independent part.

[0213] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, electronic device, or network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes various media that can store program code, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk. It should be noted that, in this article, the terms "include", "comprising" or any other variant thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device that includes a series of elements includes not only those elements, but also other elements that are not explicitly listed, or also includes elements inherent to such a process, method, article or device. Without further constraints, an element defined by the phrase "comprises a..." does not preclude the existence of additional identical elements in the process, method, article or apparatus that includes the element.

[0214] The terms used in the embodiments of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The singular forms "a", "an", "the" and "the" used in the embodiments of the present invention and the appended claims are also intended to include plural forms unless the context clearly indicates otherwise.

[0215] It should be understood that the term "and / or" as used herein is merely a description of the relationship between associated objects, indicating that three possible relationships exist. For example, "A and / or B" can represent: A exists alone, A and B exist simultaneously, or B exists alone. Furthermore, the character " / " in this document generally indicates that the associated objects are in an "or" relationship.

[0216] The word "if," as used herein, may be interpreted as "at the time of" or "when" or "in response to determining" or "in response to detecting," depending on the context. Similarly, the phrases "if it is determined" or "if (stated condition or event) is detected" may be interpreted as "when it is determined" or "in response to the determination" or "when detecting (stated condition or event)" or "in response to detecting (stated condition or event)," depending on the context.

[0217] The "first" and "second" mentioned in the embodiments are merely used to distinguish similar objects and do not represent a specific ordering of the objects. It is understood that the specific order or precedence of "first" and "second" can be interchanged where appropriate. It should be understood that the objects distinguished by "first" and "second" can be interchanged where appropriate, so that the embodiments described herein can be implemented in an order other than that illustrated or described herein.

[0218] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for identifying underdetermined operating modal parameters based on tensor decomposition recursion, characterized in that: include: The vibration response signal acquisition step is to acquire the vibration response signal X(t) of the time-varying structure collected at the selected sensor measurement point at time t under the set environmental excitation; wherein, n m is the number of sensors, N(t) is the number of time samples; The principal component score calculation step is to process the acquired vibration response signal using principal component analysis (PCA) to calculate the principal component score; The lag determination step is to calculate the autocorrelation function of the principal component scores and determine the lag τ by detecting the peak in the autocorrelation function. k ; The third-order tensor construction step is based on the modal vibration matrix Φ(t), modal response matrix Q(t) and lag period τ obtained by decomposing the vibration response signal X(t) at time t k Constructed as a dimension of n m ×N(t)×τ k The third-order tensor H(t); The third-order tensor update step adds a new data slice to the second dimension of the third-order tensor H(t) at time t to construct the third-order tensor H(t+1) at time t+1; The weighted observation matrix construction step is based on the third-order tensor H(t+1) at time t+1, introduces a dynamic forgetting factor to construct the weighted observation matrix H(t+1) at time t+1, and performs windowing and weighting processing on the observation matrix; The factor matrix recursive update step is based on the two factorizations of the weighted observation matrix H(t) at time t, the diagonalization of the connection matrix, the singular value tracking, and the common blocks between the modal response matrices Q(t) and Q(t+1), and recursively updates the modal vibration matrix Φ(t+1), the modal response matrix Q(t+1), and the signal to modal contribution matrix Θ(t+1) at time t+1; In the working modal parameter identification step, the modal vibration matrix Φ(t+1), the modal response matrix Q(t+1) and the signal to modal contribution matrix Θ(t+1) are decomposed to identify the time-varying working modal parameters.

2. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 1 is characterized in that: The principal component score calculation step specifically includes: The collected vibration response signal X(t) is standardized to obtain the standardized matrix B(t), as follows: B(t)=(X(t)-μ) / σ; Among them, μ is the mean of the original data, σ is the standard deviation; The covariance matrix Z(t) is obtained by calculating the standardized matrix B(t), and the covariance matrix Z(t) is diagonalized to obtain the matrix V(t) composed of eigenvectors, as follows: V(t)Λ(t)V T (t)=Z(t); Where T is the transpose of the matrix, Λ(t) is a diagonal matrix, The principal component score matrix S(t) is obtained based on the standardized matrix B(t) and the matrix V(t) composed of eigenvectors, as follows: S(t)=B(t)V(t).

3. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 1 is characterized in that: The hysteresis period determination step specifically includes: Calculate the autocorrelation function R(l) of the principal component score as follows: Among them, i is the index of a selected principal component, that is, it represents the i-th column of the principal component score matrix S(t), S t,i represents the value of the i-th principal component score at time t; S t+l,i Represents the value of the i-th principal component score at time t+l; represents the mean of the i-th principal component score, that is, N is the length of the data, l is the number of lags; Using the autocorrelation function R(l), the appropriate lag period τ is determined by peak detection to determine the periodic characteristics k : locs={k|R(l)>R(l-1)and R(l)>R(l+1)}; Where k represents the peak index; len(R(l)) represents the length of the autocorrelation function R(l), that is, the number of autocorrelation function values; locs represents the peak array; and len(locs) represents the number of peaks.

4. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 1 is characterized in that: The modal shape matrix Q(t), modal response matrix Q(t) and lag period τ obtained by decomposing the vibration response signal X(t) at time t k Constructed as a dimension of n m ×N(t)×τ k The third-order tensor of H (t), expressed as follows: Where R is the rank of CP decomposition; Φ r (t) is a tensor H (t) is the vector in the first dimension, i.e., the rth modal vibration shape vector matrix of the modal vibration shape matrix Φ(t); Q r (t) is a tensor H (t) is the vector in the second dimension, i.e., the rth modal response vector matrix of the modal response matrix Q(t); Θ r (t) is a tensor H (t) is a vector in the third dimension, i.e., the rth vector matrix of the signal contribution matrix Θ(t); operator Represents the vector outer product.

5. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 4 is characterized in that: The third-order tensor at time t H Add new data slices in the second dimension of (t) to construct a third-order tensor at time t+1 H (t+1), expressed as follows: H (t+1)=[ H (t),h(t+1)]; Among them, h(t+1) represents the new data slice.

6. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 5 is characterized in that: The weighted observation matrix construction step specifically includes: Calculate the data change rate change_rate of the vibration response signal input at the current time point, and dynamically adjust the forgetting factor λ(t+1) at time t+1 according to the data change rate and the forgetting factor λ(t) at time t, as follows: Among them, λ base represents the initial value of the forgetting factor; Use the dynamically adjusted forgetting factor λ(t+1) to weight the historical data and construct a new weighted observation matrix H(t+1), as follows: Among them, H (1) (t+1) is the tensor H (t+1) The matrix after expansion along the first dimension, dimension is n m ×(N(t+1)×τ k ) to extract modal information from the data using matrix decomposition techniques; is a weighted matrix consisting of the current forgetting factor, as follows: Where diag(λ(t+1),λ(t),…,1) is the diagonal matrix constructed using the forgetting factor λ(t+1); The weighted observation matrix H(t+1) is decomposed to obtain the initial modal vibration matrix Φ(t+1), the modal response matrix Q(t+1) and the signal to modal contribution matrix Θ(t+1). Let A(t+1)=Φ(t+1)⊙Θ(t+1) to obtain the loading matrix A(t+1) at time t+1; the operator ⊙ represents the Hadamard product.

7. The underdetermined operating modal parameter identification method based on tensor decomposition recursion according to claim 6 is characterized in that: The factor matrix recursive updating step specifically includes: The two factorizations of the weighted observation matrix H(t) at time t are as follows: The first decomposition is achieved by replacing the formula H (1) (t) is replaced by the result obtained by its CP decomposition, The second decomposition is the economic-scale SVD decomposition of H(t), where and U(t) is a combination of singular vectors of the tensor in the first and third dimensions; Σ(t) is a singular value in the decomposition of U(t); L(t) is a singular vector of the tensor in the second dimension; There exists a non-singular connectivity matrix So that: in, is the auxiliary matrix; the connection between the two equations W(t) is found by solving the simultaneous diagonalization of a set of matrices of size R×R; Combining the above two formulas, we can calculate E(t), L(t), W(t) and W -1 (t); W -1 (t) is the inverse matrix of W(t); Perform singular value decomposition on the loading matrix A(t+1), as follows: [I r (t+1),σ r ,F r (t+1)]=svd(A r (t1)); Among them, A r (t+1) represents the rth vector of A(t+1); Θ r (t+1) is the matrix A r The left singular vector matrix of (t+1) is A. r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the row space; σ r represents the diagonal matrix of singular values, containing A r The singular values ​​of (t+1) reflect the strength of the matrix in each orthogonal direction; Φ r (t+1) is the matrix A r The right singular vector matrix of (t+1) is A. r The eigenvector of (t+1) represents the orthogonal basis vector of the data in the column space; the updated factor matrix Θ r (t+1) and Φ r (t+1) is used to adjust the tensor decomposition results in the recursive process; Based on the singular value decomposition results, the updated L(t+1) of L(t) at time t+1 and the updated E(t+1) of E(t) at time t+1 are obtained as follows: Assign the left singular vector matrix to U(t+1) to update U(t), and assign the decomposed singular values ​​to Σ(t+1) to update Σ(t); based on Get the updated E(t+1) of E(t) at time t+1; Assign part of the right singular vector matrix to L(t+1) to update L(t); At time t+1, the formula The conversion is as follows: Using the time-shift structure of Q(t+1), time t and time t+1 are connected, where E(t+1) represents the auxiliary matrix at time t+1; W(t+1) represents the non-singular connection matrix at time t+1; W -1 (t+1) is the inverse matrix of W(t+1); Since the time-varying structures Q(t) and Q(t+1) have a common block, we write The dimensions of L(t) and L(t+1) are N(t)×R and N(t+1)×R respectively, and so are the dimensions of Q(t) and Q(t+1); in this case, the common block is Q(t) itself; Define the following matrix based on the set window: It means that L(t) at time t does not change; It means to intercept the first N(t) rows from the L(t+1) matrix to form It means extracting the N(t+1)th row of L(t+1), forming By common block get: Where λ(t) 1 / 2 represents the square root of the forgetting factor, which is used to weight historical data; get: in, represents the pseudo-inverse of a matrix; is an identity matrix, so Update of Φ(t) and Θ(t): Based on the formula Φ(t+1)⊙Θ(t+1)=E(t+1)W(t+1), update Φ(t) and Θ(t) to obtain Φ(t+1) and Θ(t+1); Update of Q(t): Based on the formula Update q T (t+1), which is finally appended to Q(t) to construct Q(t+1); Determine whether the change of Q(t+1) relative to Q(t) is less than a preset second threshold. If so, further determine whether the factor matrices of all ranks have been updated. If so, output the updated Φ(t+1), Q(t+1) and Θ(t+1).

8. An under-timed variable working modal parameter identification system, characterized in that: include: A motor system and a support platform for providing mechanical vibration and structural support; Measurement and analysis systems, including force and motion sensors, for collecting vibration response signals of time-varying structures; An excitation system, used to apply control excitations on the time-varying structure to simulate the real working environment; The processing unit analyzes the collected vibration response signal and identifies the time-varying working modal parameters based on the method according to any one of claims 1 to 7.

Images