A mechanical structure excitation identification method based on a map estimator

Abstract

The present application relates to the field of structural health monitoring, in particular to a mechanical structure excitation identification method based on a MAP estimator, which is used to solve the inverse problem of mechanical structure external excitation identification. The method comprises the following steps: S1: constructing a transfer matrix H between the excitation input position and the response measuring point; S2: measuring the monitoring mechanical structure vibration response signal y caused by the dynamic excitation; S3: establishing a MAP estimator based on the Bayesian model and the structure vibration response signal y; S4: solving the MAP estimator to obtain the mechanical structure dynamic excitation time history x. The method has the advantages of high identification accuracy, high calculation efficiency and strong stability.

Description

Technical Field

[0001] This invention relates to the field of structural health monitoring, and more specifically to a method for identifying mechanical structure excitations based on a MAP estimator. Background Technology

[0002] Dynamic load identification plays a crucial role in multiple fields, including dynamic optimization design, reliability analysis, load transfer path analysis, active vibration control, mechanical fault diagnosis, and structural health monitoring. However, research on these issues can be somewhat unpredictable without explicitly given load values.

[0003] Mechanical equipment such as aircraft engines and gas turbines may experience malfunctions during operation. Minor malfunctions can lead to economic losses, while severe malfunctions can result in catastrophic consequences. Therefore, online health monitoring and prediction of large mechanical equipment are extremely necessary, and dynamic load identification is an indispensable supporting component of this system. Wind turbine blades used in wind power generation are subjected to impacts from external objects such as sandstorms, birds, hail, and maintenance tools during operation and maintenance. These impacts can damage the turbine blades, affecting their integrity and load-bearing capacity. Therefore, accurate identification and location of impact loads are crucial. Real-time monitoring of atmospheric pulsating pressure acting on aircraft wings, as well as impacts from foreign objects such as birds and maintenance tools on wings or engine blades, is vital for ensuring the operational safety of aircraft. Identifying impact loads between train wheelsets and rails is crucial for monitoring train operation status and planning maintenance. Vibration transmission path analysis of vehicle transmission systems requires identifying the dynamic loads on engine support components to provide a basis for vehicle vibration isolation. In active vibration isolation control of marine power equipment, dynamic load identification is very important for evaluating vibration effects.

[0004] Therefore, many scientific challenges remain to be solved in the field of load identification. In dynamic environments, direct measurement methods suffer from limitations such as installation space constraints and inaccessibility of measurement locations. There is an urgent need for a method that can accurately invert and calculate loads by combining measured responses with a system mathematical model. Summary of the Invention

[0005] The purpose of this invention is to provide a mechanical structure excitation identification method based on MAP estimator, which can use easily obtainable responses such as displacement, velocity, acceleration or strain to invert the input dynamic excitation and obtain high-precision inversion of the dynamic excitation time history under different excitation forms.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] A method for identifying mechanical structure excitations based on a MAP estimator includes the following steps:

[0008] S1: Construct the transfer matrix H between the excitation input location and the response measurement point;

[0009] S2: Measure the vibration response signal y of the monitored mechanical structure caused by dynamic excitation;

[0010] S3: Establish a MAP estimator based on the Bayesian model and the structural vibration response signal y;

[0011] S4: Solve the MAP estimator to obtain the time history of dynamic excitation of the mechanical structure.

[0012] Further, step S1 includes:

[0013] S1.1: Obtain the frequency response function H(ω) between the excitation input position and the response measurement position;

[0014] S1.2: Perform an inverse fast Fourier transform on the frequency response function H(ω) obtained in step 1.1 to obtain the corresponding unit impulse response function h(t), and discretize the unit impulse response function h(t) by combining the sampling interval to obtain the transfer matrix H.

[0015] Furthermore, in step S2, the vibration response signal y of the monitored mechanical structure caused by dynamic excitation is measured by a vibration sensor.

[0016] Furthermore, the MAP estimator established based on the Bayesian model and the structural vibration response signal y includes:

[0017]

[0018] Where x is the dynamic excitation; p(xy) is the posterior probability of x given the vibration response signal y; p(y|x) is the likelihood function; and p(x) is the prior probability. MAP The maximum posterior probability of the input dynamic stimulus x.

[0019] Furthermore, step S4 includes solving the MAP estimator using the near-point gradient descent algorithm.

[0020] Furthermore, solving for the MAP estimator using the near-point gradient descent algorithm includes:

[0021] S4.1: Initialization: Let the initial value of the excitation iteration be x. (0) =0, intermediate variable z (0) =0; regularization parameter λ>0, generalized Gaussian distribution topology parameter p>0; (λ,p) is a hyperparameter pair, which can be selected by generalized cross-validation; iteration step size α=1 / maxeig(H T H), where maxeig(H) T H) represents taking matrix H TThe maximum eigenvalue of H; the iteration termination threshold ε = 10 -5 The number of iterations k=1, and the objective function is...

[0022] S4.2: Update intermediate variable z (i) :

[0023] z (k) =x (k) -αH T (Hx (k) -y);

[0024] S4.3: Update the current iteration value x (i+1) :

[0025] x (k+1) =prox(z) (k) ;λα);

[0026] Where prox(·) represents the nearest neighbor operator, specifically expressed as:

[0027]

[0028] Where b* is the solution to the first-order optimal condition, expressed as the following equation:

[0029]

[0030] In the formula Let z be the vector of the kth iteration. (k) The i-th element; b is a variable.

[0031] Furthermore, in the Bayesian model, the noise vector is a Gaussian white noise vector, i.e. The dynamic excitation follows a generalized Gaussian distribution, i.e., x ~ GG(0,σ). x I), and all elements of the Gaussian white noise vector follow an independent and identically distributed pattern:

[0032]

[0033]

[0034] Where σ1,c x ,σ x >0 and is a constant; Γ(a1,b1) and Γ(a2,b2) represent the Gamma distribution, which are the prior distributions introduced by the parameters of the Bayesian model; p is the morphological parameter of the generalized Gaussian distribution;

[0035] Substituting p(y|x) and p(x) into the MAP estimator, the MAP estimator can be further expressed as:

[0036]

[0037] Further, the solution of b > 0 is obtained through an iterative algorithm, and the specific steps are as follows:

[0038] Input:

[0039] Initialization: j = 0, t = abs(y i );

[0040] When 0 < p ≤ 1,

[0041] If Otherwise, b * is the solution generated by fixed-point iteration:

[0042]

[0043] When p > 1, b* is the solution generated by Newton - Raphson iteration:

[0044]

[0045] Output b*;

[0046] S4.4: Determine the iteration termination condition according to the following formula:

[0047]

[0048] where J is the total number of iterations; T p is the threshold during the threshold iteration; y i is the i-th element of the measurement response vector y; ‖·‖ ∞ represents the infinity norm; j is the current iteration number;

[0049] If the current solution x (k+1) satisfies the above iteration termination condition, terminate the iterative calculation and output the current solution x (k+1) as the excitation recognition result; otherwise, return to step S4.2 to continue the iterative calculation until the above iteration termination condition is satisfied, or the maximum number of iterations k max .

[0050] Compared with the prior art, the beneficial effects of the present invention are:

[0051] 1. This paper presents a mechanical structure excitation identification method based on MAP estimator, which can invert the input dynamic excitation using easily obtainable responses such as displacement, velocity, acceleration, or strain, and obtain high-precision inversion of the dynamic excitation time history under different excitation forms. 2. A novel excitation identification method is derived from the perspective of Bayesian maximum a posteriori estimation. A novel algorithm is developed to solve the excitation identification model, enabling reconstruction of the time history of different types of excitations. 3. An excitation sparse regularization identification model is established using a non-convex regularization term with convexity preservation, which can obtain sparser identification results and improve excitation identification accuracy. The computational efficiency of the algorithm is improved through a non-convex optimization acceleration strategy. Attached Figure Description

[0052] Figure 1 This is an overall flowchart of a mechanical structure excitation identification method based on MAP estimator;

[0053] Figure 2 This is a flowchart illustrating the specific steps of a mechanical structure excitation identification method based on a MAP estimator.

[0054] Figure 3 This is a schematic diagram of a typical cantilever beam structure excitation recognition device in an embodiment of the present invention;

[0055] Figure 4 This is the acceleration response signal of the response measurement point in the embodiment of the present invention;

[0056] Figure 5 This is a comparison between the excitation identification result of the cantilever beam structure and the actual excitation signal in an embodiment of the present invention. Detailed Implementation

[0057] See Figure 1 A mechanical structure excitation identification method based on MAP estimator includes the following steps:

[0058] S1: Construct the transfer matrix H between the excitation input location and the response measurement point;

[0059] S2: Measure the vibration response signal y of the monitored mechanical structure caused by dynamic excitation;

[0060] S3: Establish a MAP estimator based on the Bayesian model and the structural vibration response signal y;

[0061] S4: Solve the MAP estimator to obtain the time history of dynamic excitation of the mechanical structure.

[0062] Further, step S1 includes:

[0063] S1.1: Obtain the frequency response function H(ω) between the excitation input position and the response measurement position; it can be obtained by performing modal tests on the monitored object structure, such as hammer impact tests / frequency sweep tests, or by obtaining it through modal analysis using a finite element model.

[0064] S1.2: Perform an inverse fast Fourier transform on the frequency response function H(ω) obtained in step 1.1 to obtain the corresponding unit impulse response function h(t), and discretize the unit impulse response function h(t) by combining the sampling interval to obtain the transfer matrix H.

[0065] In step S2, the vibration response signal y of the monitored mechanical structure caused by dynamic excitation is measured using a vibration sensor. The vibration sensor can be an accelerometer, strain gauge, etc.

[0066] The MAP estimator based on the Bayesian model and the structural vibration response signal y includes:

[0067]

[0068] Where x is the dynamic excitation; p(xy) is the posterior probability of x given the vibration response signal y; p(y|x) is the likelihood function; and p(x) is the prior probability. MAP The maximum posterior probability of the input dynamic stimulus x.

[0069] Step S4 includes solving the MAP estimator using the near-point gradient descent algorithm.

[0070] Solving the MAP estimator using the near-point gradient descent algorithm includes:

[0071] S4.1: Initialization: Let the initial value of the excitation iteration be x. (0) =0, intermediate variable z (0) =0; regularization parameter λ>0, generalized Gaussian distribution topology parameter p>0; (λ,p) is a hyperparameter pair, which can be selected by generalized cross-validation; iteration step size α=1 / maxeig(H T H), where maxeig(H) T H) represents taking matrix H T The maximum eigenvalue of H; the iteration termination threshold ε = 10 -5 The number of iterations k=1, and the objective function is...

[0072] S4.2: Update intermediate variable z (i) :

[0073] z (k) =x (k) -αH T (Hx (k) -y);

[0074] S4.3: Update the current iteration value x (i+1) :

[0075] x (k+1) = prox(z (k) ; λα);

[0076] where prox(•) represents the proximal operator, and the specific expression is:

[0077]

[0078] where b* is the solution of the first-order optimality condition, expressed as the following equation:

[0079]

[0080] In the formula is the i-th element of the k-th iteration vector z (k) ; b is a variable.

[0081] In the Bayesian model, the noise vector is a Gaussian white noise vector, that is The dynamic excitation follows a generalized Gaussian distribution, that is x ~ GG(0,σ x I), and the elements of the Gaussian white noise vector are all independently and identically distributed:

[0082]

[0083]

[0084] where, σ1,c x ,σ x > 0 and are constants; Γ(a1,b1), Γ(a2,b2) represent the Gamma distribution, which is the prior distribution introduced by the parameters of the Bayesian model; p is the shape parameter of the generalized Gaussian distribution;

[0085] Substituting p(y|x) and p(x) into the MAP estimator, the MAP estimator can be further expressed as:

[0086]

[0087] The solution of b > 0 is obtained through an iterative algorithm, and the specific steps are as follows:

[0088] Input: J = 2;

[0089] Initialization: j = 0, t = abs(y i );

[0090] When 0 < p ≤ 1,

[0091] if Otherwise, b * The solution generated by the fixed-point iteration:

[0092]

[0093] When p>1, b* is the solution generated by the Newton-Raphson iteration:

[0094]

[0095] Output b*;

[0096] S4.4: Determine the iteration termination condition according to the following formula:

[0097]

[0098] Where J is the total number of iterations; T p The threshold in the threshold iteration process; y i To measure the i-th element of the response vector y; ||·|| ∞ Represents the infinity norm; j is the current iteration number;

[0099] In this embodiment, a typical mechanical cantilever beam structure model is used as the research object to identify the excitation, such as... Figure 2 As shown, one end is fixed and constrained, while the other end is free. The composite material plate has dimensions of 500×50×6mm. Taking this cantilever beam model as the research object, the simulation verification process for identifying the excitation is as follows:

[0100] Established on composite material plates, such as Figure 2 In the coordinate system shown, R1 is the fixed point of the accelerometer, with coordinates (100, 25) mm. A simulated excitation is applied at coordinate E1, with coordinates (500, 25) mm. The frequency response function between the excitation point and the measurement point is obtained through harmonic response analysis.

[0101] The sampling frequency for obtaining the frequency response function is 2048Hz, the sampling time is 0.25s, and the data length is 512. The condition number of the transfer matrix between the excitation and response points approaches each other. The condition number is a commonly used indicator to measure the ill-conditioning of a matrix; the larger the condition number, the more ill-conditioned the matrix. Therefore, it can be concluded that the excitation recognition inverse problem in this case is severely ill-conditioned.

[0102] A simulated excitation was applied at E1, and the acceleration response signal of the structure was acquired at R1. Simultaneously, the excitation time history and strain signal time history were recorded. The sampling frequency was 2048 Hz, the sampling time was 0.25 s, and the corresponding sampling length was 512. To simulate a real environment, the signal-to-noise ratio of the response signal was set to 30 dB. Figure 3As shown, the actual time history of the real excitation and the corresponding noisy acceleration signal time history are shown in the simulation process. The excitation time history is used as a comparison object for the identification results of the mechanical structure excitation identification method based on MAP estimator.

[0103] To quantitatively evaluate the performance of the method-based approach in stimulus identification, the following relative error is defined:

[0104]

[0105] Where x and These represent the actual time history of the stimulus and the recognition result of the stimulus, respectively, where ||·||2 represents the L2 norm of the vector. For example... Figure 4 The figure shows the excitation identification results for composite material plate structures. The mechanical structure excitation identification method based on the MAP estimator has a relative error of 4.82%, accurately identifying the excitation time history.

[0106] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for identifying mechanical structure excitations based on a MAP estimator, characterized in that, Includes the following steps: S1: Construct the transfer matrix between the excitation input location and the response measurement point. ; S2: Measuring the vibration response signal of the monitoring mechanical structure caused by dynamic excitation. ; S3: Based on Bayesian model and structural vibration response signals Establish a MAP estimator: ； Where x is the dynamic excitation; Let x be the posterior probability of the vibration response signal y given the vibration response signal y. Let be the likelihood function. For prior probability, The maximum posterior probability of the input dynamic stimulus x; In the Bayesian model, the noise vector is a Gaussian white noise vector, i.e. The dynamic excitation follows a generalized Gaussian distribution, that is... Furthermore, all elements of the Gaussian white noise vector follow an independent and identically distributed pattern. ； in, And it is a constant; This represents the Gamma distribution, which is the prior distribution introduced for the parameters of the Bayesian model; p These are the morphological parameters of the generalized Gaussian distribution; Will , Substituting this into the MAP estimator, the MAP estimator can be further expressed as: ； In the formula, For regularization parameters, ; S4: Use the near-point gradient descent algorithm to solve the MAP estimator and obtain the dynamic excitation time history of the mechanical structure.

2. The mechanical structure excitation identification method based on MAP estimator according to claim 1, characterized in that: Step S1 includes: S1.1: Obtain the frequency response function between the excitation input position and the response measurement position. ; S1.2: The frequency response function obtained in step 1.1 Performing an inverse fast Fourier transform yields the corresponding unit impulse response function. And in conjunction with the sampling interval, the unit impulse response function The transfer matrix is ​​obtained after discretization. H .

3. The mechanical structure excitation identification method based on MAP estimator according to claim 1, characterized in that: In step S2, a vibration sensor is used to measure the vibration response signal of the monitoring mechanical structure caused by dynamic excitation. y .

4. The mechanical structure excitation identification method based on MAP estimator according to claim 1, characterized in that: Solving the MAP estimator using the near-point gradient descent algorithm includes: S4.1: Initialization: Set the initial values ​​for the excitation iterations. intermediate variables Regularization parameters Generalized Gaussian distribution topographic parameters ; The hyperparameter pairs can be selected by generalized cross-validation; iteration step size ,in Indicates taking the matrix Maximum eigenvalue; Iteration termination threshold Number of iterations objective function ; S4.2: Update intermediate variables : ； S4.3: Update the current iteration value : ； in, The nearest neighbor operator is specifically expressed as: ； in The solution to the first-order optimal condition is expressed by the following equation: ； In the formula For the first k The iteration vector The i One element; b For variables.

5. A mechanical structure excitation identification method based on a MAP estimator according to claim 4, characterized in that: The solution is obtained through an iterative algorithm, and the specific steps are as follows: enter: ; initialization: when hour, ; if ,otherwise, The solution generated by the fixed-point iteration: ； when hour, The solution generated by the Newton-Raphson iteration: Output ; S4.4: Determine the iteration termination condition according to the following formula: ； in, J This represents the total number of iterations. T p This refers to the threshold during the threshold iteration process. y i To measure the response vector The i One element; Represents the infinite norm; j This represents the current iteration number; If the current solution If the above iteration termination condition is met, then the iteration calculation terminates and the current solution is output. The result is used as the incentive recognition result; otherwise, return to step S4.2 to continue the iterative calculation until the above iteration termination condition is met, or the maximum number of iterations is specified in advance. .

Images