CN116186495A - Structural parameter solving method based on complete modal decomposition and random forest response surface fitting - Google Patents

Structural parameter solving method based on complete modal decomposition and random forest response surface fitting Download PDF

Info

Publication number
CN116186495A
CN116186495A CN202211090323.5A CN202211090323A CN116186495A CN 116186495 A CN116186495 A CN 116186495A CN 202211090323 A CN202211090323 A CN 202211090323A CN 116186495 A CN116186495 A CN 116186495A
Authority
CN
China
Prior art keywords
frequency
signal
modal
entropy
follows
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211090323.5A
Other languages
Chinese (zh)
Inventor
马磊明
姜斌
肖玲斐
于乾坤
郭勤涛
陆宁云
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202211090323.5A priority Critical patent/CN116186495A/en
Publication of CN116186495A publication Critical patent/CN116186495A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/27Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00Machine learning
    • G06N20/10Machine learning using kernel methods, e.g. support vector machines [SVM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/004Artificial life, i.e. computing arrangements simulating life
    • G06N3/006Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/06Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Complex Calculations (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a structural parameter solving method based on complete modal decomposition and random forest response surface fitting; firstly, obtaining the structure self-vibration frequency and the main frequency number estimated value; secondly, obtaining the minimum frequency difference, designing a frequency overlapping criterion, combining the VMD to determine the number K of decomposition components, and respectively carrying out noise reduction treatment on K+1 components including the obtained K IMF components and residual items; then, the signals after noise reduction and reconstruction are subjected to complete division by utilizing a variable speed GWO algorithm and combining with an FIR band-pass filter, an objective function is set by comprehensively considering a correlation coefficient, a mutual approximation entropy and an envelope spectrum entropy, weights of different factors are determined, optimal division frequencies are obtained, and a damping ratio is solved based on each-order single-mode damping signals; and finally, fitting a response surface by using a random forest and solving the elastic modulus and the poisson ratio of the structure by combining the modal frequencies. The method can identify the mode structure parameters with high precision in the environments with strong noise and dense modes, and improves the accuracy, robustness, reliability and practicability of the algorithm.

Description

Structural parameter solving method based on complete modal decomposition and random forest response surface fitting
Technical Field
The invention relates to the field of mechanical structure parameter identification, in particular to a structural parameter solving method based on complete modal decomposition and random forest response surface fitting.
Background
Modal parameter identification is an important means for researching the inherent dynamic characteristics of a system, and modal parameters mainly comprise modal frequency, damping ratio, elastic modulus, poisson ratio and the like. Common modal parameter identification methods include time domain methods (Ibrahim time domain (ITD), feature system implementation (ERA) and random subspace (SSI)), frequency domain methods (peak picking method and frequency domain decomposition method), time domain methods (empirical mode decomposition (EMD), singular Spectrum Analysis (SSA), analytical Mode Decomposition (AMD), variational Mode Decomposition (VMD) and SVD), and the like. The method has obvious effect on processing stable signal modal parameter identification, but has poor effect on processing strong noise and dense modal signals. How to overcome the defects of the existing algorithm and improve the identification precision of the modal structure parameters is the key point of future research.
The complete modal decomposition by combining the time-frequency domain method with the optimization algorithm can obviously improve the modal frequency and damping ratio recognition accuracy in consideration of strong noise and dense modal signals; parameters such as the elastic modulus of the structure, the Poisson's ratio and the like can be accurately obtained by adopting a random forest response surface fitting algorithm.
Disclosure of Invention
In view of the above, the invention aims to provide a structural parameter solving method based on complete modal decomposition and random forest response surface fitting, which can identify modal structural parameters with high precision in environments with strong noise and dense modalities, and improves accuracy, robustness, reliability and practicability of an algorithm.
In order to achieve the above purpose, the present invention adopts the following technical scheme:
a structural parameter solving method based on complete modal decomposition and random forest response surface fitting comprises the following steps:
s1, acquiring a structure main frequency number estimated value and a self-vibration frequency estimated value, and estimating a frequency division range and a minimum frequency difference of adjacent frequencies according to a spectrogram;
s2, constructing a frequency overlapping criterion and a threshold value thereof based on the minimum frequency difference of the adjacent frequencies obtained in the step S1, and determining a VMD decomposition optimal component value K and a main frequency number according to the frequency overlapping criterion;
s3, decomposing the original modal signal based on the optimal component value K obtained in the step S2, respectively carrying out noise reduction treatment on K+1 signal components including K IMF components and residual items obtained by VMD decomposition, and carrying out reconstruction on the K+1 signal components after noise reduction to obtain a noise reduction reconstruction signal;
S4, the noise reduction reconstruction signal obtained in the step S3 is completely divided by combining a variable speed GWO algorithm with an FIR band-pass filter, the division value of the band-pass filter is continuously adjusted by utilizing a variable speed GWO algorithm, the optimal division frequency is obtained, the correlation coefficient, the mutual approximation entropy and the envelope spectrum entropy are comprehensively considered by an optimization objective function, and the weights of different factors are determined by utilizing an entropy weight method, so that the optimal division frequency is obtained;
s5, acquiring a single-mode complete attenuation signal based on the optimal segmentation frequency obtained in the step S4, and solving a mode frequency and a damping ratio;
and S6, acquiring the elastic modulus and the Poisson' S ratio of the structure by using a random forest fit response surface based on the modal frequency obtained in the step S5.
Preferably, the step S1 includes:
s101, acquiring a structure main frequency number estimated value by combining an autoregressive power spectrum AR model with an automatic peak value picking algorithm
Figure BDA0003836752540000021
And the natural frequency estimate +.>
Figure BDA0003836752540000022
Autoregressive power spectrum P of signal y (n) AR (e ) Calculated from the following formula:
Figure BDA0003836752540000023
wherein p is R The order of the AR model; n=1, 2, …, N representing the nth data point; n is the signal data length;
Figure BDA0003836752540000024
is the variance of y (n); />
Figure BDA0003836752540000025
Corresponding model parameters; r is (r) R (k R ) Is an autocorrelation function;
step S102, estimating a frequency division range according to the spectrogram
Figure BDA0003836752540000026
And minimum frequency difference between adjacent frequencies
Figure BDA0003836752540000027
The following respectively satisfy:
Figure BDA0003836752540000028
Figure BDA0003836752540000029
in the method, in the process of the invention,
Figure BDA00038367525400000210
and->
Figure BDA00038367525400000211
Respectively estimating values of frequencies of the ith order and the jth order; min {.cndot }, represents taking the minimum value.
Preferably, the step S2 includes:
step S201, the VMD decomposition process is as follows:
in the case of constraint that the sum of the eigenmode functions is equal to the original input signal y (t), k eigenmode functions u are sought k (t) minimizing the sum of the estimated bandwidths of each eigenmode function; decomposing for the original input signal:
Figure BDA00038367525400000212
wherein y (t) is a multicomponent signal to be decomposed; u (u) k (t) is a single component signal obtained after VMD decomposition; a is that k (t) is u k (t) an instantaneous amplitude;
Figure BDA00038367525400000213
is u k The instantaneous phase of (t); the constraint variant model expression for VMD decomposition is as follows:
Figure BDA0003836752540000031
in the formula, { u k }={u 1 ,...,u K Sum { omega } k }={ω 1 ,...,ω K Each of which is a sequence of all modes and its center frequency; delta (t) is the dirac distribution;
Figure BDA0003836752540000032
is the derivative with respect to time; * Is convolution operation; />
Figure BDA0003836752540000033
L as vector 2 A norm;
Figure BDA0003836752540000034
is the instantaneous frequency; VMD decomposes the last term into residual terms, satisfying:
Figure BDA0003836752540000035
wherein ru (t) is a residual term; k is VMD decomposition modal number;
step S202, based on the adjacent frequency minimum frequency difference obtained in step S1
Figure BDA0003836752540000036
The frequency overlap criterion and its threshold are constructed as follows:
Figure BDA0003836752540000037
Wherein omega is l Representing the center frequency of the first component obtained by VMD decomposition; r represents the center frequency overlap ratio; setting a threshold value of a frequency overlapping criterion, if the change of r is larger than the threshold value, judging that modal splitting occurs, otherwise, judging that VMD splitting meets the requirement, and selecting the maximum value in the number K of all components meeting the threshold value as the number of optimal VMD splitting components;
step S203, determining, based on step S202, that the main frequency is:
n * =K+1 (8)。
preferably, the step S3 includes:
step S301, decomposing the original modal signal based on the number K of the optimal VMD decomposition components obtained in the step S2, and respectively carrying out noise reduction processing on K+1 components including K IMF components and residual items obtained by VMD decomposition by utilizing SVD;
step S302, selecting the SVD noise reduction effective rank as the order2 times of the number of main frequencies in the source vibration signal, and reconstructing K+1 signal components after noise reduction to obtain y f (t)。
Preferably, the step S4 includes:
s401, performing complete division on the noise reduction reconstruction signal obtained in the step S3 by combining a variable speed GWO algorithm with an FIR band-pass filter; structural main frequency n obtained based on step S2 * And corresponding natural frequency estimation value
Figure BDA0003836752540000038
Establishing a speed change GWO optimization parameter as +.>
Figure BDA0003836752540000039
Wherein D is i =[d i1 ,d i2 ]For the left and right division frequencies of the i-th order mode, i=1, 2, …, n * -1;d i2 -d i1 Dividing the bandwidth for the boundary frequency of the ith order mode; the left and right cut-off frequencies of the FIR band-pass filter determined based on the variable speed GWO algorithm satisfy:
Figure BDA0003836752540000041
based on (9) according to the previous n respectively * -the order 1 left and right cut-off frequencies are filtered with FIR band-pass filters and corresponding optimal left and right cut-off frequencies are determined according to an optimized objective function of the variable speed GWO algorithm;
before n * -1 st order cut-off frequency and corresponding single mode time domain signal after determination * The order Shan Motai time domain signal is calculated as follows:
Figure BDA0003836752540000042
step S402, optimizing an objective function by utilizing a variable speed GWO algorithm, comprehensively considering a correlation coefficient, a mutual approximation entropy and an envelope spectrum entropy, determining weights of different factors by utilizing an entropy weight method, and obtaining an optimal segmentation frequency, wherein the correlation calculation process is as follows:
the velocity and location component update formulas are as follows:
Figure BDA0003836752540000043
in the method, in the process of the invention,
Figure BDA0003836752540000044
is the moving speed of the gray wolf population; />
Figure BDA0003836752540000045
Is an inertial factor; c VSGWO1 ,c VSGWO2 And c VSGWO3 Learning factors, respectively; r is (r) VSGWO1 =random(0,1),r VSGWO2 =random (0, 1) and r VSGWO3 =random (0, 1) is a random number, respectively; />
Figure BDA0003836752540000046
Is the current position of the gray wolf; t is t GWO Is the number of iterations; />
Figure BDA0003836752540000047
And->
Figure BDA0003836752540000048
Intermediate auxiliary variables respectively;
optimizing objective function setting and comprehensively considering correlation coefficient, mutual approximation entropy and envelope spectrum entropy, determining different factor weights by utilizing entropy weight method, obtaining optimal segmentation frequency, and reconstructing signal y for noise reduction f The ith single frequency signal y of (t) i (t) calculating the entropy of its envelope spectrum and its signal y respectively ri (t)=y f (t)-y i Correlation coefficient and mutual approximation entropy of (t), i=1, 2, …, n *
Preferably, the reconstructed signal y is reconstructed for noise reduction f The ith single frequency signal y of (t) i (t) calculating the entropy of its envelope spectrum and its signal y respectively ri (t)=y f (t)-y i The correlation calculation process of the correlation coefficient and the mutual approximation entropy of (t) is as follows:
(1) The correlation coefficient calculation process is as follows:
Figure BDA0003836752540000049
where Cov (·) represents the covariance of the variables; var [. Cndot. ] is the variance of the signal; re (·) is the correlation coefficient of the signal; the larger the correlation coefficient is, the larger the correlation between the two signals is;
(2) The mutual approximation entropy calculation process is as follows:
sequence y i (t) and y j (t)(i=1,2,...,n * ;j=1,2,...,n * The method comprises the steps of carrying out a first treatment on the surface of the i+.j) are normalized separately, i.e.:
y i * (i n )=[y i (i n )-mean(y i )]/SD(y i ) (13)
y j * (i m )=[y j (i m )-mean(y j )]/SD(y j ) In the formula (14), SD (. Cndot.) represents the standard deviation of the variables; i.e n And i m Respectively represent the ith n And i m Data points, and i n =1,2,...,N;i m =1,2,...,N;i n ≠i m The method comprises the steps of carrying out a first treatment on the surface of the mean (·) represents averaging the variables; then a similar tolerance of ry=0.2 COV (y i ,y j ) Calculate a new sequence { y } i * (i n ),y j * (i m ) Mutual approximation entropy of { as original sequence { y } i (i n ),y j (i m ) Mutual approximation entropy;
for a time series containing N data points, an N binary distance matrix D is calculated t ,D t Is the ith of (2) n Line i m Column elements are written as
Figure BDA0003836752540000051
The expression is as follows:
Figure BDA0003836752540000052
using matrix D t The element in (2) is calculated to obtain the component association coefficient
Figure BDA0003836752540000053
And->
Figure BDA0003836752540000054
Figure BDA0003836752540000055
Figure BDA0003836752540000056
Adding the results of the diagonal sums of each row based on the results of the formulas (16) and (17), and averaging to obtain
Figure BDA0003836752540000057
And->
Figure BDA0003836752540000058
The further resulting mutual approximation entropy formula is as follows:
Figure BDA0003836752540000059
wherein CrossAPen (m, ry) is the mutual approximation entropy;
(3) Calculating the component y obtained by signal decomposition i (t) envelope spectrum entropy;
(4) And determining weights of different factors by using an entropy weight method.
Preferably, the specific step of determining weights of different factors by using the entropy weight method comprises the following steps:
firstly, the evaluation index is subjected to homodromous, standardized and a standardized matrix χ= { X is obtained 1 ,X 2 ,...,X α Sum of weight vector omega τ ={W 1 ,W 2 ,...,W β }, whereinAlpha and beta are the dimension of the data to be evaluated and the number of evaluation indexes respectively;
secondly, obtaining a matrix by carrying out homodromous on index attributes in the original data set to obtain E, constructing a weighted normalized matrix Z, wherein Z is formed by E and omega τ After multiplication, the following products are obtained:
Figure BDA00038367525400000510
wherein:
Figure BDA0003836752540000061
is j th ω The weights of the attributes are determined by an entropy weight method, and the specific process of the entropy weight method is as follows:
calculating each element in the probability matrix P
Figure BDA0003836752540000062
The following are provided:
Figure BDA0003836752540000063
the entropy of each index is calculated separately as follows:
Figure BDA0003836752540000064
the weight coefficients are calculated as follows:
Figure BDA0003836752540000065
then, the distances between the indexes of each sample and the positive and negative ideal solutions are calculated, namely the proximity degree between each evaluation object and the optimal scheme and the worst scheme is calculated:
The positive ideal solution calculation formula is as follows:
Figure BDA0003836752540000066
the negative ideal solution calculation formula is as follows:
Figure BDA0003836752540000067
the distance calculation formula of each sample index from the forward solution is as follows:
Figure BDA0003836752540000068
the distance calculation formula of each sample index distance negative solution is as follows:
Figure BDA0003836752540000069
finally, calculating the closeness degree of each evaluation object and the optimal scheme
Figure BDA00038367525400000610
And according to->
Figure BDA00038367525400000611
Sequencing the sizes to give an evaluation result;
Figure BDA00038367525400000612
wherein:
Figure BDA0003836752540000071
Figure BDA0003836752540000072
a closer to 1 indicates a better sample score.
Preferably, the step S5 includes:
step S501, acquiring a complete modal damping signal y based on the optimal segmentation frequency obtained in the step S4 i (t),i=1,2,…,n * The envelope spectrum is obtained by utilizing Hilbert transformation, and meanwhile, the cross-correlation function calculation process based on the NExT method is as follows:
when externally excited
Figure BDA0003836752540000073
Scope system k d At the point of point, ith p Point and j p The cross-correlation function for the point correspondence is calculated as follows:
Figure BDA0003836752540000074
in the method, in the process of the invention,
Figure BDA0003836752540000075
is of structure i p The r of the measuring point d A mode shape of the order mode; />
Figure BDA0003836752540000076
Representing the excitation point k d And modal order r d A related constant term; />
Figure BDA0003836752540000077
Is of structure j p The s < th > of the measuring point d A mode shape of the order mode; />
Figure BDA0003836752540000078
Representing the excitation point k d And modal order s d A related constant term; />
Figure BDA0003836752540000079
Is the r d Characteristic values of the order vibration mode; />
Figure BDA00038367525400000710
Is the s th d Characteristic values of the order vibration mode; n (N) d Is the degree of freedom of the system; p is p d 、q d And τ both represent system delay; e [. Cndot.]A cross-correlation function representing the variables;
Step S502, the modal frequency and the damping ratio are obtained based on the envelope of the least square fitting cross-correlation function.
Preferably, the step S6 includes:
the method comprises the steps that a data set theta is formed by mu groups of structure simulation results under different elastic moduli and Poisson ratios, wherein the structure simulation results are obtained based on finite element software Patran, the data set theta is fitted by using a random forest algorithm with the elastic moduli and Poisson ratios as abscissa and the modal frequencies as ordinate, and the structure elastic moduli and Poisson ratios are obtained by combining experimental modal results;
step S601, fitting a response surface by using a random forest, wherein the specific process is as follows:
(1) Coding mu group of model simulation result data Θ calculated based on finite element software Patran, and taking the data Θ as input of a random forest algorithm;
(2) Randomly replacing the extracted training set and the test set from the data set theta by using a Bootstrap method; the training set is taken as a root node to split and train a plurality of regression trees;
(3) Splitting nodes in the regression tree according to the minimum mean square error principle until the nodes cannot be split further, repeating the operation to generate a plurality of regression trees, and forming the trees into a forest
Figure BDA00038367525400000711
Wherein (1)>
Figure BDA00038367525400000712
Is a random quantity; ρ is the sample; t is the number of regression trees;
(4) A series of independent decision trees with the same distribution are constructed, voting is carried out, the final category of the sample is determined, and the obtained prediction result is as follows:
Figure BDA00038367525400000713
wherein: v (ρ) is a random forest combination classification prediction result;
Figure BDA00038367525400000714
classifying models for decision trees;Y T Output variables for random forests; i (·) is an indication function; />
Figure BDA00038367525400000715
To maximize the variable Y T Is a value of (2);
(5) Fitting a response surface based on a random forest model by taking the elastic modulus of the structure and the poisson ratio as abscissa and the modal frequency as ordinate;
step S602, based on the modal frequency obtained in step S5, the response surface fitted in step S601 and experimental modal results, the elastic modulus of the structure and the Poisson' S ratio are obtained.
The beneficial effects of the invention are as follows:
1. the invention can avoid frequency spectrum leakage by utilizing the autoregressive power spectrum AR model, can effectively process noise-containing signals or nonstationary signals, and can obtain accurate structure main frequency and self-vibration frequency estimation values by combining an automatic peak value pickup algorithm;
2. constructing a frequency overlapping criterion and a threshold value thereof by utilizing the minimum frequency difference of adjacent frequencies, and simultaneously determining an optimal component value by combining the VMD and the frequency overlapping criterion, thereby effectively solving the problem of mode order determination and having high recognition precision for dense modes in a strong interference environment;
3. According to the invention, the original modal signal is decomposed based on the number of the extremely optimal components, all IMF components and residual items obtained by VMD decomposition are subjected to noise reduction processing by utilizing SVD respectively, so that the phenomenon of modal loss caused by the existence of VMD decomposition residual items can be effectively avoided, and meanwhile, the noise reduction of each component can be realized on the premise of ensuring the original characteristics of the modal signal, so that the energy loss caused by the noise reduction processing is reduced to the greatest extent, and the noise reduction effect of an algorithm is improved;
4. the speed component of the traditional particle swarm optimization algorithm is introduced into the gray wolf optimization algorithm to form the variable speed GWO optimization algorithm, the advantages of strong local searching capability of the gray wolf optimization algorithm and the advantages of high convergence speed and strong global searching capability of the particle swarm optimization algorithm are combined, meanwhile, the comprehensive consideration of the correlation coefficient, the mutual approximation entropy and the envelope spectrum entropy is set for the optimization objective function, different factor weights are determined by utilizing the entropy weight method, and the performance and the optimization result of the optimization algorithm can be remarkably improved.
5. The invention utilizes the speed change GWO algorithm to combine with the FIR band-pass filter to complete division of the single-mode cut-off frequency, so that an optimized modal frequency division result can be obtained, the dense modal identification precision is effectively improved, and an accurate single-mode attenuation vibration signal is obtained;
6. According to the invention, the random forest fitting response surface is utilized to obtain the structural elastic modulus and poisson ratio, the random forest is not easy to be subjected to the fitting phenomenon, the capability of processing high-dimensional data and unbalanced data is strong, and the higher response surface fitting precision can be obtained.
Drawings
Fig. 1 is a flow chart of a structural parameter solving method based on complete modal decomposition and random forest response surface fitting provided in embodiment 1;
FIG. 2 is a bearing block used in example 1, wherein a is a bearing block physical diagram, and b is a bearing block modal test geometric modeling diagram;
FIG. 3 is a plot of the original time domain vibration signal and power of the modal test provided in example 1, wherein a is the original time domain vibration signal and b is the AR autoregressive power plot;
fig. 4 is a VMD component number influence diagram provided in embodiment 1, wherein a represents a VMD component number versus center frequency repetition rate maximum influence diagram, and b represents a VMD component number versus energy loss coefficient influence diagram;
FIG. 5 is an AR autoregressive power spectrum obtained based on the VMD and in combination with the optimal component numbers provided in example 1;
FIG. 6 is a graph of the SVD-based noise reduction reconstructed time domain vibration signal and power provided in example 1, wherein a is the SVD noise reduction reconstructed time domain signal and b is the AR autoregressive power spectrum;
FIG. 7 is the optimized objective function value based on the shift GWO provided in example 1;
fig. 8 is a first-order modal damping signal and a power spectrum obtained based on the optimal split frequency provided in embodiment 1, wherein a is a first-order modal time-domain damping signal, and b is a first-order modal AR autoregressive power spectrum;
fig. 9 is a second-order modal attenuation signal and a power spectrum obtained based on the optimal division frequency provided in embodiment 1, wherein a is a second-order modal time-domain attenuation signal, and b is a second-order modal AR autoregressive power spectrum;
fig. 10 is a third-order modal damping signal and a power spectrum obtained based on the optimal split frequency provided in embodiment 1, wherein a is a third-order modal time-domain damping signal, and b is a third-order modal AR autoregressive power spectrum;
fig. 11 is a fourth-order mode attenuation signal and a power spectrum obtained based on the optimal division frequency provided in embodiment 1, wherein a is a fourth-order mode time domain attenuation signal, and b is a fourth-order mode AR autoregressive power spectrum;
fig. 12 is a fifth-order mode attenuation signal and a power spectrum obtained based on the optimal split frequency provided in embodiment 1, wherein a is a fifth-order mode time domain attenuation signal, and b is a fifth-order mode AR autoregressive power spectrum;
Fig. 13 is a sixth-order modal damping signal and power spectrum obtained based on the optimal split frequency provided in embodiment 1, wherein a is a sixth-order modal time-domain damping signal, and b is a sixth-order modal AR autoregressive power spectrum;
fig. 14 is a plot of the random forest response surface fitting results provided in example 1.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are some embodiments of the present invention, but not all embodiments of the present invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
Example 1
Referring to fig. 1-14, the present embodiment provides a structural parameter solving method based on complete modal decomposition and random forest response surface fitting, and the flow of the method is shown in fig. 1, in the present embodiment, the free modal test of the bearing seat is taken as an example to verify the algorithm performance, fig. 2 a shows a bearing seat physical diagram, and fig. 2 b shows a bearing seat modal test geometric modeling diagram. The strategy specifically comprises the following steps:
S1, acquiring a structure main frequency number estimated value by combining an autoregressive power spectrum AR model with an automatic peak value picking algorithm
Figure BDA0003836752540000091
And the natural frequency estimate +.>
Figure BDA0003836752540000092
Estimating the frequency division range from the spectrogram>
Figure BDA0003836752540000093
And minimum frequency difference between adjacent frequencies->
Figure BDA0003836752540000101
Specifically, in this embodiment, the step S1 includes:
s101, acquiring a structure main frequency number estimated value by combining an autoregressive power spectrum AR model with an automatic peak value picking algorithm
Figure BDA0003836752540000102
And the natural frequency estimate +.>
Figure BDA0003836752540000103
Autoregressive power spectrum P of signal y (n) AR (e ) Can be calculated by the following formula:
Figure BDA0003836752540000104
wherein p is R The order of the AR model; n (n=1, 2,., N) is the nth data point; n is the signal data length;
Figure BDA0003836752540000105
representing the variance of the white noise sequence; />
Figure BDA0003836752540000106
Corresponding model parameters; r is (r) R (k R ) Is an autocorrelation function; the original time domain signal of the modal test and the AR autoregressive power spectrum thereof are shown in FIG. 3; estimated dominant frequency +.>
Figure BDA0003836752540000107
The estimated values of the self-vibration frequencies are respectively
Figure BDA0003836752540000108
Figure BDA0003836752540000109
And
Figure BDA00038367525400001010
step S102, estimating a frequency division range according to the spectrogram
Figure BDA00038367525400001011
And minimum frequency difference between adjacent frequencies
Figure BDA00038367525400001012
The following respectively satisfy:
Figure BDA00038367525400001013
Figure BDA00038367525400001014
in the method, in the process of the invention,
Figure BDA00038367525400001015
and->
Figure BDA00038367525400001016
Respectively estimating values of frequencies of the ith order and the jth order; min { · } represents taking the minimum value; the present embodiment can obtain the minimum frequency difference of adjacent frequencies +.>
Figure BDA00038367525400001017
/>
Step S2, obtaining the minimum frequency difference of adjacent frequencies based on the step S1
Figure BDA00038367525400001018
Constructing a frequency overlap criterion and a threshold value thereof, determining a VMD decomposition optimal component value K by utilizing the VMD and the frequency overlap criterion, and determining a main frequency value n * =K+1;
Specifically, in this embodiment, the step S2 includes:
step S201, the VMD decomposition process is as follows:
in the case of constraint that the sum of the eigenmode functions is equal to the original input signal y (t), k eigenmode functions u are sought k (t) minimizing the sum of the estimated bandwidths of each eigenmode function; decomposing for the original input signal:
Figure BDA00038367525400001019
wherein y (t) is a multicomponent signal to be decomposed; u (u) k (t) is a single component signal obtained after VMD decomposition; a is that k (t) is u k (t) an instantaneous amplitude;
Figure BDA00038367525400001020
is u k The instantaneous phase of (t); the constraint variant model expression for VMD decomposition is as follows:
Figure BDA0003836752540000111
in the formula, { u k }={u 1 ,...,u K Sum { omega } k }={ω 1 ,...,ω K Respectively is } isA modal sequence and a center frequency thereof; delta (t) is the dirac distribution;
Figure BDA0003836752540000112
is the derivative with respect to time; * Is convolution operation; I.I 2 L as vector 2 A norm;
Figure BDA0003836752540000113
is the instantaneous frequency; VMD decomposes the last term into residual terms, satisfying:
Figure BDA0003836752540000114
wherein ru (t) is a residual term; k is VMD decomposition modal number;
step S202, obtaining the minimum frequency difference of adjacent frequencies based on step S1
Figure BDA0003836752540000115
The frequency overlap criterion and its threshold are constructed as follows:
Figure BDA0003836752540000116
Wherein omega is l Representing the center frequency of the VMD decomposition first component; r represents the center frequency overlap ratio; setting a threshold value of a frequency overlapping criterion as 1, if the change of r is greater than 1, judging that modal splitting occurs, otherwise, judging that VMD splitting meets the requirement, and selecting the maximum value in the number K of all components meeting the threshold value as the optimal VMD splitting component number; as can be seen from fig. 4 a, in this embodiment, the optimal decomposition result of the VMD can be obtained by taking k=5, and the relevant values are shown in table 1;
table 1 center frequency repetition rate maximum
Component number 2 3 4 5 6 7 8
Repetition rate 0 1 1 1 136 88 884
Component number 9 10 11 12 13 14 15
Repetition rate 401 7 99 1014 10038 944 30
Component number 16 17 18 19 20
Repetition rate 935 3858 8667 961 78
Fig. 4 b shows a graph of energy loss coefficient variation with the number of components; the corresponding results are shown in table 2:
TABLE 2 energy loss coefficient
Component number 2 3 4 5 6 7 8
Loss coefficient 0.0745 0.0502 0.0307 0.0124 0.0121 0.0076 0.0075
Component number 9 10 11 12 13 14 15
Loss coefficient 0.0064 0.0060 0.0045 0.0041 0.0040 0.0031 0.0029
Component number 16 17 18 19 20
Loss coefficient 0.0024 0.0022 0.0023 0.0018 0.0017
Step S203, determining, based on step S202, that the main frequency is:
n * =K+1 (8)
s3, decomposing the original modal signal based on the number K of the optimal components obtained in the step S2, respectively carrying out noise reduction treatment on K+1 components including K IMF components and residual terms obtained by VMD decomposition by SVD, selecting the order of SVD noise reduction effective rank to be 2 times of the number of main frequencies in the source vibration signal, and reconstructing the K+1 signal components after noise reduction to obtain y f (t);
Specifically, in this embodiment, the step S3 includes:
in step S301, the original modal signal is decomposed based on the number K of the optimal components obtained in step S2, fig. 5 shows an AR autoregressive power spectrum obtained based on the VMD and combined with the number K of the optimal components, as can be seen from comparing fig. 5 with fig. 3 b, the VMD causes modal loss due to the existence of the residual term, which brings a larger error to the recognition result of the modal parameter, and the modal loss result is shown in table 3:
table 3 comparison of Modal loss results
Figure BDA0003836752540000121
In order to avoid the phenomenon of modal loss, K+1 components including K IMF components and residual terms obtained by VMD decomposition are subjected to noise reduction processing by SVD in the subsequent processing, and the basic principle of SVD noise reduction is as follows:
for the kth noise-containing component u of the K+1 components k (t) performing phase space reconstruction, constructing the phase space reconstruction as a p×q order Hankel matrix:
Figure BDA0003836752540000122
wherein n=p+q-1, and p is not less than q. For H m Singular value decomposition is carried out to obtain:
H m =U∑V T (10)
wherein U and V are unitary matrices of p×p and q×q, respectively; sigma is p×q diagonal matrix, satisfying the following equation:
∑=diag(λ 12 ,...,λ s ) (11)
wherein lambda is 12 ,...,λ s For matrix H m And lambda is a singular value of 1 ≥λ 2 ≥...≥λ s ≥0;
Before using r s The larger order singular values represent the true signal, followed by s-r s The smaller order singular values represent the noise signal; performing inverse operation on matrix components corresponding to the singular values of the previous r-order to obtain H m Is the best approximation matrix of (a)
Figure BDA0003836752540000123
Matrix->
Figure BDA0003836752540000124
Adding and averaging the anti-diagonal elements in the two to obtain a noise-reduced signal;
step S302, selecting the order of SVD noise reduction effective rank to be 2 times of the number of main frequencies in the source vibration signal, and reconstructing the K+1 noise reduced signal components to obtain y f (t); fig. 6 a shows a reconstructed time domain signal based on SVD noise reduction, fig. 6 b shows an AR autoregressive power spectrum corresponding to a in fig. 6, and comparison with b in fig. 5 and 3 shows that, by comprehensively considering residual signals for VMDs based on optimal component numbers, a mode loss phenomenon can be effectively avoided and a better noise reduction effect can be obtained, and comparison results are shown in table 4:
table 4 comparison of modal loss results
Figure BDA0003836752540000131
S4, the noise reduction reconstruction signal obtained in the step S3 is completely divided by combining a variable speed GWO algorithm with an FIR band-pass filter, the division value of the band-pass filter is continuously adjusted by utilizing a variable speed GWO algorithm, the optimal division frequency is obtained, the correlation coefficient, the mutual approximation entropy and the envelope spectrum entropy are comprehensively considered by an optimization objective function, and the weights of different factors are determined by utilizing an entropy weight method, so that the optimal division frequency is obtained;
specifically, in this embodiment, the step S4 includes:
step S402, the noise reduction reconstruction signal obtained in the step S3 is completely divided by combining a variable speed GWO algorithm with an FIR band-pass filter; based on step S2, a structure main frequency of n can be obtained * And a natural frequency estimate
Figure BDA0003836752540000132
Establishing a speed change GWO optimization parameter as +.>
Figure BDA0003836752540000133
Wherein D is i =[d i1 ,d i2 ](i=1,2,....,n * -1) left and right segmentation frequencies for the i-th order mode; d, d i2 -d i1 Dividing the bandwidth for the boundary frequency of the ith order mode; the FIR filter left and right cut-off frequencies determined based on the shift GWO algorithm satisfy:
Figure BDA0003836752540000134
based on (12) according to the previous n respectively * The-1 order left and right cut-off frequencies are filtered by using an FIR band-pass filter and the corresponding optimal left and right cut-off frequencies are determined according to the optimized objective function of the speed change GWO algorithm, and the FIR band-pass filter expression is as follows:
Figure BDA0003836752540000135
where y (n) is the filter input; y is r (n) is the filter output; h (k) is a filter coefficient; n (N) J Is the filter order;
before n * -1 st order cut-off frequency and corresponding single mode time domain signal, n < th > after determination * The order Shan Motai time domain signal is calculated as follows:
Figure BDA0003836752540000136
step S402, optimizing an objective function by utilizing a variable speed GWO algorithm, comprehensively considering a correlation coefficient, a mutual approximation entropy and an envelope spectrum entropy, determining weights of different factors by utilizing an entropy weight method, and obtaining an optimal segmentation frequency, wherein the correlation calculation process is as follows:
the speed change GWO algorithm combines the advantages of strong local searching capability of the wolf optimization algorithm and the advantages of high convergence speed and strong global searching capability of the particle swarm optimization algorithm, and introduces the speed component of the traditional particle swarm optimization algorithm into the wolf optimization algorithm to form the speed change wolf optimization algorithm, and the fused speed and position component update formula is as follows:
Figure BDA0003836752540000141
In the method, in the process of the invention,
Figure BDA0003836752540000142
is the moving speed of the gray wolf population; />
Figure BDA0003836752540000143
Is an inertial factor; c VSGWO1 ,c VSGWO2 And c VSGWO3 Learning factors, respectively; r is (r) VSGWO1 =random(0,1),r VSGWO2 =random (0, 1) and r VSGWO3 =random (0, 1) is a random number, respectively;
Figure BDA0003836752540000144
is the current position of the gray wolf; t is t GWO Is the number of iterations; />
Figure BDA0003836752540000145
And->
Figure BDA0003836752540000146
Intermediate auxiliary variables respectively;
optimization objective function optimization process based on speed change GWO is shown in fig. 7, and optimization values of different modes are shown in table 5:
table 5 objective function optimization values for different modalities
Mode order 1 2 3 4 5 6
Objective function value 1.8526e-09 1.0567e-06 1.3494e-09 3.4882-09 3.2264-08 3.7036e-09
Optimizing objective function setting and comprehensively considering correlation coefficient, mutual approximation entropy and envelope spectrum entropy, determining different factor weights by utilizing entropy weight method, obtaining optimal segmentation frequency, and reconstructing signal y for noise reduction f The ith single frequency signal y of (t) i (t)(i=1,2,....,n * ) Respectively calculating the envelope spectrum entropy and the signal y thereof ri (t)=y f (t)-y i The correlation coefficient and the mutual approximation entropy of (t) are calculated as follows: the correlation coefficient calculation process of (1) is as follows:
Figure BDA0003836752540000147
wherein Cov (·) is the signal covariance; var [. Cndot. ] is the variance of the signal; re (·) is the correlation coefficient of the signal; the larger the correlation coefficient is, the larger the correlation between the two signals is;
(2) The mutual approximation entropy calculation process is as follows:
if sequence y i (t) and y j (t)(i=1,2,...,n * ;j=1,2,…,n * The method comprises the steps of carrying out a first treatment on the surface of the i.noteq.j), the sequence y is preferably chosen to be larger in magnitude i (t) and y j (t 0 is normalized separately, namely:
y i * (i n )=[y i (i n )-mean(y i )]/SD(y i ) (17)
y j * (i m )=[y j (i m )-mean(y j )]/SD(y j ) In the formula (18), i n And i m Respectively represent the ith n And i m Data points, and i n =1,2,...,N;i m =1,2,...,N;i n ≠i m The method comprises the steps of carrying out a first treatment on the surface of the mean (·) represents averaging the variables; then the similar tolerance ry=0.2 COV (y i ,y j ) Calculate a new sequence { y } i * (i n ),y j * (i m ) Mutual approximation entropy of { as original sequence { y } i (i n ),y j (i m ) Mutual approximation entropy;
for a time series containing N data points, an N binary distance matrix D is calculated t ,D t Is the ith of (2) n Line i m Column elements are written as
Figure BDA0003836752540000151
The expression is as follows:
Figure BDA0003836752540000152
using matrix D t The element in the model can be conveniently calculated to obtain the component association coefficient
Figure BDA0003836752540000153
And->
Figure BDA0003836752540000154
Figure BDA0003836752540000155
/>
Figure BDA0003836752540000156
The result of each row of oblique line summation is added up based on the results of the formulas (20) and (21) to obtain an average value
Figure BDA0003836752540000157
And->
Figure BDA0003836752540000158
Further available mutual approximation entropy formulas are as follows:
Figure BDA0003836752540000159
wherein CrossAPen (m, ry) is the mutual approximation entropy;
(3) Calculating the component y obtained by signal decomposition i (t)(i=1,2,...,n * ) Envelope spectrum entropy;
(4) Determining weights of different factors by using entropy weight method
Firstly, the evaluation index is subjected to homodromous, standardized and a standardized matrix χ= { X is obtained 1 ,X 2 ,...,X α Sum of weight vector omega τ ={W 1 ,W 2 ,...,W β -wherein α and β are the dimension of the data to be evaluated and the number of evaluation indices, respectively;
and secondly, carrying out homodromous on index attributes in the original data set to obtain a matrix to obtain E, and constructing a weighted normalized matrix Z. Z is composed of E and omega τ After multiplication, the following products are obtained:
Figure BDA00038367525400001510
wherein:
Figure BDA00038367525400001511
is j th ω The weight (importance degree) of each attribute is determined by an entropy weight method, and the specific process of the entropy weight method is as follows:
calculating each element in the probability matrix P
Figure BDA00038367525400001512
The following are provided:
Figure BDA00038367525400001513
the entropy of each index is calculated separately as follows:
Figure BDA0003836752540000161
the weight coefficients are calculated as follows:
Figure BDA0003836752540000162
then, the distances between the indexes of each sample and the positive and negative ideal solutions are calculated, namely the proximity degree between each evaluation object and the optimal scheme and the worst scheme is calculated:
the positive ideal solution calculation formula is as follows:
Figure BDA0003836752540000163
the negative ideal solution calculation formula is as follows:
Figure BDA0003836752540000164
the distance calculation formula of each sample index from the forward solution is as follows:
Figure BDA0003836752540000165
the distance calculation formula of each sample index distance negative solution is as follows:
Figure BDA0003836752540000166
finally, calculating the closeness degree of each evaluation object and the optimal scheme
Figure BDA0003836752540000167
And according to->
Figure BDA0003836752540000168
Sequencing the sizes to give an evaluation result;
Figure BDA0003836752540000169
wherein:
Figure BDA00038367525400001610
Figure BDA00038367525400001611
a closer to 1 indicates a better sample score; the sample index weights are shown in table 6:
table 6 index weights of samples
Sample index Correlation coefficient Mutual approximation entropy Envelope spectrum entropy
Weighting of 0.6732 0.1245 0.2023
Step S5, acquiring a complete modal damping signal y based on the optimal segmentation frequency obtained in the step S4 i (t)(i=1,2,....,n * ) Obtaining an envelope spectrum by utilizing Hilbert transformation, and simultaneously obtaining a modal frequency and a damping ratio based on a natural excitation technique (NExT) method and a least square method;
Specifically, in this embodiment, the step S5 includes:
step S501, acquiring a complete modal damping signal y based on the optimal segmentation frequency obtained in the step S4 i (t)(i=1,2,....,n * ) The envelope spectrum is obtained by utilizing Hilbert transformation, and meanwhile, the cross-correlation function calculation process based on the NExT method is as follows:
when externally excited
Figure BDA0003836752540000171
Scope system k d When the point is reached, the corresponding cross-correlation function of the ith point and the jth point is calculated as follows:
Figure BDA0003836752540000172
in the method, in the process of the invention,
Figure BDA0003836752540000173
is of structure i p The r of the measuring point d A mode shape of the order mode; />
Figure BDA0003836752540000174
Representing the excitation point k d And modal order r d A related constant term; />
Figure BDA0003836752540000175
Is of structure j p The s < th > of the measuring point d A mode shape of the order mode; />
Figure BDA0003836752540000176
Representing the excitation point k d And modal order s d A related constant term; />
Figure BDA0003836752540000177
Is the r d Characteristic values of the order vibration mode; />
Figure BDA0003836752540000178
Is the s th d Characteristic values of the order vibration mode; n (N) d Is the degree of freedom of the system; p is p d 、q d And τ both represent system delay; e [. Cndot.]A cross-correlation function representing the variables; the acquired complete modal attenuation signals of each order and the corresponding AR autoregressive power spectrum thereof are shown in figures 8-13; the optimal division frequencies of each order are shown in table 7:
TABLE 7 optimal fractional frequency Bandwidth for each order
Mode order 1 2 3 4 5
Frequency bandwidth 10.4323 12.0741 14.6987 36.4313 26.1941
Left cut-off frequency 290.3119 511.1822 539.9800 1336.0760 1744.3150
Right cut-off frequency 300.7442 523.2563 554.6787 1372.5073 1770.5091
Step S502, obtaining modal frequency and damping ratio based on the envelope of the least square fitting cross-correlation function;
The identification results of the optimal modal frequencies and damping ratios of each order are shown in table 8:
table 8 results of modal frequency identification for each order
Mode order 1 2 3 4 5 6
Modal frequencies 295.8021 517.2863 548.3121 1354.8045 1757.4623 1891.2932
The damping ratio recognition results and errors of each stage are shown in table 9:
table 9 damping ratio identification results for each stage
Figure BDA0003836752540000179
S6, acquiring the elastic modulus and the Poisson' S ratio of the structure by using a random forest fit response surface based on the modal frequency obtained in the step S5;
specifically, in the present embodiment, this step S6:
fitting a mu group of mode simulation result data set theta obtained based on Patran calculation, and matching by combining an experimental mode result to obtain a structural elastic modulus and a Poisson ratio;
step S601, fitting a response surface by using a random forest, wherein the specific process is as follows:
(1) Coding mu group of model simulation result data Θ obtained based on Patran calculation, and taking the data Θ as input of a random forest algorithm;
(2) Randomly replacing and extracting 70% of the data set theta to be a training set and 30% to be a test set by using a Bootstrap method; the training set is taken as a root node to split and train a plurality of regression trees;
(3) Splitting nodes in the regression tree according to the minimum mean square error principle until the nodes cannot be split further, and repeating the operationGenerating multiple regression trees, forming forest from the trees
Figure BDA0003836752540000181
Wherein (1)>
Figure BDA0003836752540000182
Is a random quantity; ρ is the sample; t is the number of regression trees;
(4) A series of independent decision trees with the same distribution are constructed, voting is carried out, the final category of the sample is determined, and the obtained prediction result is as follows:
Figure BDA0003836752540000183
wherein: v (ρ) is a random forest combination classification prediction result;
Figure BDA0003836752540000184
classifying the model for a decision tree; y is Y T Output variables for random forests; i (·) is an indication function; />
Figure BDA0003836752540000185
To maximize the variable Y T Is a value of (2);
(5) Fitting a response surface based on a random forest model by taking the elastic modulus of the structure and the poisson ratio as abscissa and the modal frequency as ordinate; the response surface fitting results are shown in fig. 14, and the response surface fitting accuracy of each order is shown in table 10:
table 10 accuracy of the fitting of the response surfaces of each order
Mode order 1 2 3 4 5 6
Root mean square error RMSE 0.0015 0.0008 0.0013 0.0002 0.0001 0.0002
Fitting accuracy R2 0.9847 0.9864 0.9760 0.9119 0.9951 0.9924
Step S602, based on the modal frequency obtained in step S5, the response surface fitted in step S601 and experimental modal results, the elastic modulus of the structure and the Poisson' S ratio are obtained.
Supplementary data
The present invention is not described in detail in the present application, and is well known to those skilled in the art.
The foregoing describes in detail preferred embodiments of the present invention. It should be understood that numerous modifications and variations can be made in accordance with the concepts of the invention by one of ordinary skill in the art without undue burden. Therefore, all technical solutions which can be obtained by logic analysis, reasoning or limited experiments based on the prior art by the person skilled in the art according to the inventive concept shall be within the scope of protection defined by the claims.

Claims (9)

1. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting is characterized by comprising the following steps of:
s1, acquiring a structure main frequency number estimated value and a self-vibration frequency estimated value, and estimating a frequency division range and a minimum frequency difference of adjacent frequencies according to a spectrogram;
s2, constructing a frequency overlapping criterion and a threshold value thereof based on the minimum frequency difference of the adjacent frequencies obtained in the step S1, and determining a VMD decomposition optimal component value K and a main frequency number according to the frequency overlapping criterion;
s3, decomposing the original modal signal based on the optimal component value K obtained in the step S2, respectively carrying out noise reduction treatment on K+1 signal components including K IMF components and residual items obtained by VMD decomposition, and carrying out reconstruction on the K+1 signal components after noise reduction to obtain a noise reduction reconstruction signal;
s4, the noise reduction reconstruction signal obtained in the step S3 is completely divided by combining a variable speed GWO algorithm with an FIR band-pass filter, the division value of the band-pass filter is continuously adjusted by utilizing a variable speed GWO algorithm, the optimal division frequency is obtained, the correlation coefficient, the mutual approximation entropy and the envelope spectrum entropy are comprehensively considered by an optimization objective function, and the weights of different factors are determined by utilizing an entropy weight method, so that the optimal division frequency is obtained;
S5, acquiring a single-mode complete attenuation signal based on the optimal segmentation frequency obtained in the step S4, and solving a mode frequency and a damping ratio;
and S6, acquiring the elastic modulus and the Poisson' S ratio of the structure by using a random forest fit response surface based on the modal frequency obtained in the step S5.
2. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S1 comprises:
s101, acquiring a structure main frequency number estimated value by combining an autoregressive power spectrum AR model with an automatic peak value picking algorithm
Figure FDA0003836752530000011
And the natural frequency estimate +.>
Figure FDA0003836752530000012
Autoregressive power spectrum P of signal y (n) AR (e ) Calculated from the following formula:
Figure FDA0003836752530000013
wherein p is R The order of the AR model; n=1, 2, …, N representing the nth data point; n is the signal data length;
Figure FDA0003836752530000014
is the variance of y (n); />
Figure FDA0003836752530000015
Corresponding model parameters; r is (r) R (k R ) Is an autocorrelation function;
step S102, estimating a frequency division range according to the spectrogram
Figure FDA0003836752530000016
And minimum frequency difference between adjacent frequencies->
Figure FDA0003836752530000017
The following respectively satisfy:
Figure FDA0003836752530000021
Figure FDA0003836752530000022
in the method, in the process of the invention,
Figure FDA0003836752530000023
and->
Figure FDA0003836752530000024
Respectively estimating values of frequencies of the ith order and the jth order; min {.cndot }, represents taking the minimum value.
3. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S2 comprises:
Step S201, the VMD decomposition process is as follows:
in the case of constraint that the sum of the eigenmode functions is equal to the original input signal y (t), k eigenmode functions u are sought k (t) minimizing the sum of the estimated bandwidths of each eigenmode function; decomposing for the original input signal:
Figure FDA0003836752530000025
wherein y (t) is a multicomponent signal to be decomposed; u (u) k (t) is a single component signal obtained after VMD decomposition; a is that k (t) is u k (t) an instantaneous amplitude;
Figure FDA0003836752530000026
is u k The instantaneous phase of (t); the constraint variant model expression for VMD decomposition is as follows:
Figure FDA0003836752530000027
in the formula, { u k }={u 1 ,...,u K Sum { omega } k }={ω 1 ,...,ω K Respectively all mode sequences and the mode sequencesA heart frequency; delta (t) is the dirac distribution;
Figure FDA0003836752530000028
is the derivative with respect to time; * Is convolution operation; />
Figure FDA0003836752530000029
L as vector 2 A norm; />
Figure FDA00038367525300000210
Is the instantaneous frequency; VMD decomposes the last term into residual terms, satisfying:
Figure FDA00038367525300000211
wherein ru (t) is a residual term; k is VMD decomposition modal number;
step S202, based on the adjacent frequency minimum frequency difference obtained in step S1
Figure FDA00038367525300000212
The frequency overlap criterion and its threshold are constructed as follows:
Figure FDA00038367525300000213
wherein omega is l Representing the center frequency of the first component obtained by VMD decomposition; r represents the center frequency overlap ratio; setting a threshold value of a frequency overlapping criterion, if the change of r is larger than the threshold value, judging that modal splitting occurs, otherwise, judging that VMD splitting meets the requirement, and selecting the maximum value in the number K of all components meeting the threshold value as the number of optimal VMD splitting components;
Step S203, determining, based on step S202, that the main frequency is:
n * =K+1 (8)。
4. the structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S3 includes:
step S301, decomposing the original modal signal based on the number K of the optimal VMD decomposition components obtained in the step S2, and respectively carrying out noise reduction processing on K+1 components including K IMF components and residual items obtained by VMD decomposition by utilizing SVD;
step S302, selecting the order of SVD noise reduction effective rank to be 2 times of the number of main frequencies in the source vibration signal, and reconstructing the K+1 noise reduced signal components to obtain y f (t)。
5. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S4 comprises:
s401, performing complete division on the noise reduction reconstruction signal obtained in the step S3 by combining a variable speed GWO algorithm with an FIR band-pass filter; structural main frequency n obtained based on step S2 * And corresponding natural frequency estimation value
Figure FDA0003836752530000031
Establishing a speed change GWO optimization parameter as +.>
Figure FDA0003836752530000037
Wherein D is i =[d i1 ,d i2 ]For the left and right division frequencies of the i-th order mode, i=1, 2, …, n * -1;d i2 -d i1 Dividing the bandwidth for the boundary frequency of the ith order mode; the left and right cut-off frequencies of the FIR band-pass filter determined based on the variable speed GWO algorithm satisfy: / >
Figure FDA0003836752530000032
Based on (9) according to the previous n respectively * The cut-off frequency of the-1 order is filtered by an FIR band-pass filter and is changed according to the changeThe optimization objective function of the speed GWO algorithm determines the corresponding optimal left and right cut-off frequencies;
before n * -1 st order cut-off frequency and corresponding single mode time domain signal after determination * The order Shan Motai time domain signal is calculated as follows:
Figure FDA0003836752530000033
step S402, optimizing an objective function by utilizing a variable speed GWO algorithm, comprehensively considering a correlation coefficient, a mutual approximation entropy and an envelope spectrum entropy, determining weights of different factors by utilizing an entropy weight method, and obtaining an optimal segmentation frequency, wherein the correlation calculation process is as follows:
the velocity and location component update formulas are as follows:
Figure FDA0003836752530000034
in the method, in the process of the invention,
Figure FDA0003836752530000035
is the moving speed of the gray wolf population; />
Figure FDA0003836752530000036
Is an inertial factor; c VSGWO1 ,c VSGWO2 And c VSGWO3 Learning factors, respectively; r is (r) VSGWO1 =random(0,1),r VSGWO2 =random (0, 1) and r VSGWO3 =random (0, 1) is a random number, respectively;
Figure FDA0003836752530000041
is the current position of the gray wolf; t is t GWO Is the number of iterations; />
Figure FDA0003836752530000042
And->
Figure FDA0003836752530000043
Intermediate auxiliary variables respectively;
optimizing objective function setting and comprehensively considering correlation coefficient, mutual approximation entropy and envelope spectrum entropy, determining different factor weights by utilizing entropy weight method, obtaining optimal segmentation frequency, and reconstructing signal y for noise reduction f The ith single frequency signal y of (t) i (t) calculating the entropy of its envelope spectrum and its signal y respectively ri (t)=y f (t)-y i Correlation coefficient and mutual approximation entropy of (t), i=1, 2, …, n *
6. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 5, wherein the structural parameter solving method is used for denoising reconstruction signal y f The ith single frequency signal y of (t) i (t) calculating the entropy of its envelope spectrum and its signal y respectively ri (t)=y f (t)-y i The correlation calculation process of the correlation coefficient and the mutual approximation entropy of (t) is as follows:
(1) The correlation coefficient calculation process is as follows:
Figure FDA0003836752530000044
where Cov (·) represents the covariance of the variables; var [. Cndot. ] is the variance of the signal; re (·) is the correlation coefficient of the signal; the larger the correlation coefficient is, the larger the correlation between the two signals is;
(2) The mutual approximation entropy calculation process is as follows:
sequence y i (t) and y j (t)(i=1,2,...,n * ;j=1,2,...,n * The method comprises the steps of carrying out a first treatment on the surface of the i+.j) are normalized separately, i.e.:
y i * (i n )=[y i (i n )-mean(y i )]/SD(y i ) (13)
y j * (i m )=[y j (i m )-mean(y j )]/SD(y j ) (14)
wherein S isD (·) represents the standard deviation of the variables; i.e n And i m Respectively represent the ith n And i m Data points, and i n =1,2,...,N;i m =1,2,...,N;i n ≠i m The method comprises the steps of carrying out a first treatment on the surface of the mean (·) represents averaging the variables; then a similar tolerance of ry=0.2 COV (y i ,y j ) Calculate a new sequence { y } i * (i n ),y j * (i m ) Mutual approximation entropy of { as original sequence { y } i (i n ),y j (i m ) Mutual approximation entropy;
for a time series containing N data points, an N binary distance matrix D is calculated t ,D t Is the ith of (2) n Line i m Column elements are written as
Figure FDA0003836752530000049
The expression is as follows:
Figure FDA0003836752530000045
using matrix D t The element in (2) is calculated to obtain the component association coefficient
Figure FDA0003836752530000046
And->
Figure FDA0003836752530000047
Figure FDA0003836752530000048
Figure FDA0003836752530000051
Adding the results of the diagonal sums of each row based on the results of the formulas (16) and (17), and averaging to obtain
Figure FDA0003836752530000052
And->
Figure FDA0003836752530000053
The further resulting mutual approximation entropy formula is as follows:
Figure FDA0003836752530000054
wherein CrossAPen (m, ry) is the mutual approximation entropy;
(3) Calculating the component y obtained by signal decomposition i (t) envelope spectrum entropy;
(4) And determining weights of different factors by using an entropy weight method.
7. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 6, wherein the specific step of determining weights of different factors by using an entropy weight method comprises the following steps:
firstly, the evaluation index is subjected to homodromous, standardized and a standardized matrix χ= { X is obtained 1 ,X 2 ,...,X α Sum of weight vector omega τ ={W 1 ,W 2 ,...,W β -wherein α and β are the dimension of the data to be evaluated and the number of evaluation indices, respectively;
secondly, obtaining a matrix by carrying out homodromous on index attributes in the original data set to obtain E, constructing a weighted normalized matrix Z, wherein Z is formed by E and omega τ After multiplication, the following products are obtained:
Figure FDA0003836752530000055
wherein:
Figure FDA0003836752530000056
is j th ω The weights of the attributes are determined by an entropy weight method, and the index weights are determined by the entropy weight methodThe method comprises the following specific processes:
calculating each element in the probability matrix P
Figure FDA0003836752530000057
The following are provided:
Figure FDA0003836752530000058
the entropy of each index is calculated separately as follows:
Figure FDA0003836752530000059
The weight coefficients are calculated as follows:
Figure FDA00038367525300000510
then, the distances between the indexes of each sample and the positive and negative ideal solutions are calculated, namely the proximity degree between each evaluation object and the optimal scheme and the worst scheme is calculated:
the positive ideal solution calculation formula is as follows:
Figure FDA0003836752530000061
the negative ideal solution calculation formula is as follows:
Figure FDA0003836752530000062
the distance calculation formula of each sample index from the forward solution is as follows:
Figure FDA0003836752530000063
the distance calculation formula of each sample index distance negative solution is as follows:
Figure FDA0003836752530000064
finally, calculating the closeness degree of each evaluation object and the optimal scheme
Figure FDA0003836752530000065
And according to->
Figure FDA0003836752530000066
Sequencing the sizes to give an evaluation result;
Figure FDA0003836752530000067
wherein:
Figure FDA0003836752530000068
Figure FDA0003836752530000069
a closer to 1 indicates a better sample score.
8. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S5 comprises:
step S501, acquiring a complete modal damping signal y based on the optimal segmentation frequency obtained in the step S4 i (t),i=1,2,…,n * The envelope spectrum is obtained by utilizing Hilbert transformation, and meanwhile, the cross-correlation function calculation process based on the NExT method is as follows:
when externally excited
Figure FDA00038367525300000617
Scope system k d At the point of point, ith p Point and j p The cross-correlation function for the point correspondence is calculated as follows:
Figure FDA00038367525300000610
in the method, in the process of the invention,
Figure FDA00038367525300000611
is of structure i p The r of the measuring point d A mode shape of the order mode; / >
Figure FDA00038367525300000612
Representing the excitation point k d And modal order r d A related constant term; />
Figure FDA00038367525300000613
Is of structure j p The s < th > of the measuring point d A mode shape of the order mode; />
Figure FDA00038367525300000614
Representing the excitation point k d And modal order s d A related constant term; />
Figure FDA00038367525300000615
Is the r d Characteristic values of the order vibration mode; />
Figure FDA00038367525300000616
Is the s th d Characteristic values of the order vibration mode; n (N) d Is the degree of freedom of the system; p is p d 、q d And τ both represent system delay; e [. Cndot.]A cross-correlation function representing the variables;
step S502, the modal frequency and the damping ratio are obtained based on the envelope of the least square fitting cross-correlation function.
9. The structural parameter solving method based on complete modal decomposition and random forest response surface fitting according to claim 1, wherein the step S6 includes:
the method comprises the steps that a data set theta is formed by mu groups of structure simulation results under different elastic moduli and Poisson ratios, wherein the structure simulation results are obtained based on finite element software Patran, the data set theta is fitted by using a random forest algorithm with the elastic moduli and Poisson ratios as abscissa and the modal frequencies as ordinate, and the structure elastic moduli and Poisson ratios are obtained by combining experimental modal results;
step S601, fitting a response surface by using a random forest, wherein the specific process is as follows:
(1) Coding mu group of model simulation result data Θ calculated based on finite element software Patran, and taking the data Θ as input of a random forest algorithm;
(2) Randomly replacing the extracted training set and the test set from the data set theta by using a Bootstrap method; the training set is taken as a root node to split and train a plurality of regression trees;
(3) Splitting nodes in the regression tree according to the minimum mean square error principle until the nodes cannot be split further, repeating the operation to generate a plurality of regression trees, and forming the trees into a forest
Figure FDA0003836752530000071
Wherein (1)>
Figure FDA0003836752530000072
Is a random quantity; ρ is the sample; t is the number of regression trees;
(4) A series of independent decision trees with the same distribution are constructed, voting is carried out, the final category of the sample is determined, and the obtained prediction result is as follows:
Figure FDA0003836752530000073
wherein: v (ρ) is a random forest combination classification prediction result;
Figure FDA0003836752530000074
classifying the model for a decision tree; y is Y T Output variables for random forests; i (·) is an indication function; />
Figure FDA0003836752530000075
To maximize the variable Y T Is a value of (2);
(5) Fitting a response surface based on a random forest model by taking the elastic modulus of the structure and the poisson ratio as abscissa and the modal frequency as ordinate;
step S602, based on the modal frequency obtained in step S5, the response surface fitted in step S601 and experimental modal results, the elastic modulus of the structure and the Poisson' S ratio are obtained.
CN202211090323.5A 2022-09-07 2022-09-07 Structural parameter solving method based on complete modal decomposition and random forest response surface fitting Pending CN116186495A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211090323.5A CN116186495A (en) 2022-09-07 2022-09-07 Structural parameter solving method based on complete modal decomposition and random forest response surface fitting

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211090323.5A CN116186495A (en) 2022-09-07 2022-09-07 Structural parameter solving method based on complete modal decomposition and random forest response surface fitting

Publications (1)

Publication Number Publication Date
CN116186495A true CN116186495A (en) 2023-05-30

Family

ID=86444837

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211090323.5A Pending CN116186495A (en) 2022-09-07 2022-09-07 Structural parameter solving method based on complete modal decomposition and random forest response surface fitting

Country Status (1)

Country Link
CN (1) CN116186495A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116401983A (en) * 2023-06-07 2023-07-07 湖南泛联新安信息科技有限公司 Automatic mapping method for signals driven by simulated main frequency
CN117708547A (en) * 2024-02-05 2024-03-15 西安热工研究院有限公司 Novel turbine unit vibration signal processing method and system

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116401983A (en) * 2023-06-07 2023-07-07 湖南泛联新安信息科技有限公司 Automatic mapping method for signals driven by simulated main frequency
CN116401983B (en) * 2023-06-07 2023-09-22 湖南泛联新安信息科技有限公司 Automatic mapping method for signals driven by simulated main frequency
CN117708547A (en) * 2024-02-05 2024-03-15 西安热工研究院有限公司 Novel turbine unit vibration signal processing method and system
CN117708547B (en) * 2024-02-05 2024-04-30 西安热工研究院有限公司 Method and system for processing vibration signals of steam turbine unit

Similar Documents

Publication Publication Date Title
CN116186495A (en) Structural parameter solving method based on complete modal decomposition and random forest response surface fitting
JP2020017293A (en) Method for forming dynamic model for machine behavior from detection data
CN104807534B (en) Equipment eigentone self study recognition methods based on on-line vibration data
CN114399032B (en) Method and system for predicting metering error of electric energy meter
CN104737229A (en) Method for transforming input signal
CN112861066B (en) Machine learning and FFT (fast Fourier transform) -based blind source separation information source number parallel estimation method
CN110782041B (en) Structural modal parameter identification method based on machine learning
CN108228978B (en) Xgboost time sequence prediction method combined with complementary set empirical mode decomposition
CN110207689A (en) A kind of pulsar signal denoising and discrimination method based on Wavelet Entropy
CN116629431A (en) Photovoltaic power generation amount prediction method and device based on variation modal decomposition and ensemble learning
Salman et al. Comparative Study of QPSO and other methods in Blind Source Separation
CN112001256B (en) Mixed signal power frequency interference removing method and system
CN117056669A (en) Denoising method and system for vibration signal
Kalantari et al. Time series imputation via l 1 norm-based singular spectrum analysis
CN113177078B (en) Approximate query processing algorithm based on condition generation model
Darojah et al. The training of feedforward neural network using the unscented Kalman filter for voice classification application
CN112399366A (en) Indoor positioning method based on Hankel matrix and WKNN variance extraction
Carvajal et al. An identification method for Errors-in-Variables systems using incomplete data
CN107315713B (en) One-dimensional signal denoising and enhancing method based on non-local similarity
Paulson et al. Integrated distance estimators for linear models applied to some published data sets
Gang et al. Towards automated single channel source separation using neural networks
CN112649869A (en) Reservoir characteristic parameter prediction method and system based on GA-WNN
Weerapong et al. Cluster validity index for big Data based on density discriminant analysis
JP6636969B2 (en) Signal estimation device, method, and program
Chen et al. Identification of both coefficients and orders for ARMAX system

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination