Impact response spectrum time domain acceleration impact signal generation method based on wavelet superposition

Abstract

The invention discloses an impact response spectrum time domain acceleration impact signal generation method based on wavelet superposition. The method comprises the following steps: calculating a wavelet group signal; calculating the fitness value of each particle individual; calculating a global optimal position; optimizing wavelet group signal parameters; and calculating a time domain acceleration impact signal. According to the method, in the process of generating the time domain acceleration impact signal through the impact response spectrum, the half-sine parameter does not need to be rounded in the optimization process, and precision loss caused by rounded is avoided. The method does not depend on actual experimental data, and can generate the time-domain acceleration impact signal from the impact response spectrum under the condition of no experimental data support. According to the method, the constraint of displacement and speed is introduced in the process of generating the time domain acceleration impact signal by the impact response spectrum, so that the generated acceleration impact signal not only can accurately fit the target impact response spectrum, but also can effectively control the maximum displacement and the maximum speed of the impact signal, and the reliability and the safety of an impact vibration experiment are ensured.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of impact vibration test, and particularly relates to a shock response spectrum time-domain acceleration shock signal generation method based on wavelet superposition. BACKGROUND

[0002] The shock response spectrum reproduction test is an important means for checking the original performance capability of a test piece under an impact vibration environment, and is widely applied in the fields of aviation and aerospace. The shock response spectrum describes the maximum response of a single degree of freedom system with different natural frequencies under a specific impact excitation, and can comprehensively reflect the influence of an impact signal on the system. However, only a standard shock response spectrum curve is provided in the existing standard documents, and a time-domain acceleration shock signal meeting the spectrum requirement is not given. In the existing method for generating an acceleration shock signal from a shock response spectrum, a time-domain acceleration shock signal is mainly obtained from a large number of impact vibration tests, and is screened, superposed or adjusted in parameters so as to meet the requirement of the target shock response spectrum. When experimental data cannot be obtained in advance, an acceleration shock signal cannot be generated from the shock response spectrum. In the existing impact vibration standard, only the error spectrum is constrained, and no explicit constraint is proposed for the displacement signal and the velocity signal in the impact vibration, which leads to problems such as displacement or velocity overrun in the test process.

[0003] In order to overcome the dependence on the test data, in recent years, a time-domain acceleration waveform generation method based on wavelet superposition is proposed, an acceleration waveform matched with the target shock response spectrum is constructed by optimizing parameters of multiple wavelets, and the efficiency of the acceleration signal generation is significantly improved. In order to ensure that the acceleration of each wavelet is zero at the starting and ending time, the half-sine number of the wavelet is usually set to an odd number. However, in the optimization process of the traditional algorithm, the matching precision between the generated signal and the target spectrum is reduced by means of approximate rounding and the like.

[0004] Chinese Patent CN201710158663.X discloses a shock response spectrum time-domain waveform matching method based on a genetic algorithm, and proposes a method for synthesizing a time-domain acceleration signal meeting the requirement of a target shock response spectrum by optimizing wavelet parameters by means of the genetic algorithm. In the method, the normalized half-sine parameter and the half-sine parameter are processed by rounding, which leads to a decrease in the optimization precision of the half-sine number parameter, and further affects the precision of the synthesized acceleration signal. SUMMARY

[0005] To solve the above problems in the prior art, the present application aims to design a shock response spectrum time-domain acceleration shock signal generation method based on wavelet superposition, which can improve the precision of the synthesized acceleration signal.

[0006] In order to achieve the above object, the technical solution of the present invention is as follows: a method for generating a time-domain acceleration impact signal of an impact response spectrum based on wavelet superposition, comprising the following steps:

[0007] A. Initialization

[0008] A1. Setting the target shock response spectrum

[0009] Define the target shock response spectrum as R t , R t For an n×2 matrix, the expression is:

[0010] R t =[H,L]

[0011] Where H represents the frequency vector of the target impact response spectrum, H=[h1,h2,...,h n ] T is an n×1 column vector with the unit of Hz; L represents the amplitude vector of the target impact response spectrum, L=[l1,l2,...,l n ] T , is an n×1 column vector with unit of m / s 2 ; n is the number of points in the target spectrum, and the superscript T represents the vector transpose. Parameters h1,h2,...,h n ,l1,l2,...,l n Set up on site by engineers.

[0012] A2. Define the parameters of the wavelet superposition optimization algorithm

[0013] A2-1. Set basic parameters

[0014] The basic parameters include the number of particles N0, the maximum number of iterations T max and threshold ε, where N0, T max are all integers greater than 100; ε takes a value in the range of 0 to 1. The basic parameters are all set on-site by engineers.

[0015] A2-2. Setting the global optimal position

[0016] Define the global optimal position matrix GB(p) at the pth iteration. GB(p) is an n×3 matrix. The specific expression is:

[0017]

[0018] Among them, the element gb in the i-th column and j-th row i,j (p) When p=0, the initial value is 0, that is, gb i,j (0) = 0; p is the number of iterations, i = 1, 2, 3, j = 1, 2, ..., n, p = 1, 2, ..., T max .

[0019] Define matrix PB (k) (p), PB (k) (p) is an n x 3 matrix, and its specific expression is:

[0020]

[0021] where the element in the i-th column and the j-th row is The initial value is 0 when p = 0, that is, k is the particle number, k = 1, 2,..., N0.

[0022] When p = 0, define the parameter and Y g The initial values of (p) and (p) are 0 and 1, respectively, that is, Y g (0) = 1.

[0023] A2-3, set the half-sine parameter matrix

[0024] Define the half-sine parameter value matrix of the k-th particle in the p-th iteration as Δ (k) (p), the motion trend matrix of the k-th particle in the p-th iteration is VN (k) (p), the individual memory matrix of the k-th particle in the p-th iteration is PBN (k) (p), the global peak state matrix in the p-th iteration is GBN(p); matrix Δ (k) (p), VN (k) (p), PBN (k) (p), GBN(p) are all n x L M matrices, and the elements of matrices Δ (k) (p), PBN (k) (p) and GBN(p) are 0 or 1, and the parameter L M is calculated as follows:

[0025]

[0026] where the parameters L min , L max are positive odd numbers, and need to satisfy L max > L min , L min ≥ 3, which are set by engineers on site.

[0027] The expression of matrix Δ (k) (p) is:

[0028]

[0029] where the element in the o-th column and the j-th row is The initial value is 0 when p = 0, that is,

[0030] matrix VN (k) The expression of (p) is:

[0031]

[0032] wherein the element in the oth column and jth row is The initial value when p=0 is a random number in the interval [VM min ,VM max ]; the parameter VM min is in the range of -6 to -4, and the parameter VM max is in the range of 4 to 6, which are set by engineers on site.

[0033] matrix PBN (k) The expression of (p) is:

[0034]

[0035] wherein the element in the oth column and jth row is The initial value when p=0 is 0, that is

[0036] The expression of matrix GBN(p) is:

[0037]

[0038] wherein the element in the oth column and jth row is gbn o,j The initial value when p=0 is 0, that is gbn o,j (0)=0.

[0039] A2-4, setting particle position matrix

[0040] The kth particle position matrix in the pth iteration is defined as Z (k) (p), and Z (k) (p) is an n×3 matrix, and the specific expression is:

[0041]

[0042] wherein the initial value of the element in the 1st column and jth row when p=0 is a random number in the interval [U min ,U max ]; the element in the 2nd column and jth row is a random number in the interval [0,t m ); and the initial value of the element in the 3rd column and jth row is a random positive odd number in the interval [L min ,L max ].min , U max are positive numbers, and satisfy U max > U min , t m is the acceleration signal duration, in seconds, t m is a positive number. Parameters U max , U min , t m are set by engineers on site.

[0043] A2-5, set particle motion step length matrix

[0044] Define the kth particle motion step length matrix V (k) (p) in the pth iteration (k) (p) is an n x 3 matrix, and the expression is:

[0045]

[0046] where, when p = 0, the first column, jth row element initial value is a random number in the interval [VA min , VA max ]; the second column, jth row element initial value is a random number in the interval [Vt min , Vt max ]; the third column element initial value is 0; Parameters VA min , Vt min are negative numbers, VA max , Vt max are positive numbers, and are set by engineers on site.

[0047] B, calculate wavelet group signal

[0048] In the pth iteration, calculate the acceleration signal amplitude from the kth particle position matrix Z (k) (p), and the calculation formula is:

[0049]

[0050] where, M (k) (p) is the amplitude sequence of the kth particle in the pth iteration, in m / s 2 , which is an nt x 1 column vector; h j is the jth element of the target spectral frequency vector H set in step A-1; k q ∈ {1, 2,..., nt}, and the calculation formula of parameter nt is:

[0051]

[0052] Ts is the sampling time in seconds, which is set by the engineer on site, and <·> represents the upward rounding operation.

[0053] C, calculate the fitness value of each particle individual C1, calculate the spectrum error

[0054] is the acceleration signal amplitude sequence of the kth particle at the pth iteration (k) (p), according to the impact response spectrum analysis method in the national standard GB / T 29716.4-2018 "Mechanical vibration and shock signal processing Part 4: Impact response spectrum analysis", calculate the amplitude of the jth frequency point of the target spectrum in m / s 2 , calculate the error spectrum θ according to the following formula (k) (p) in m / s 2 :

[0055]

[0056] wherein, l j is the jth element of the target spectrum amplitude vector L set in step A-1; |·| represents taking the absolute value.

[0057] C2, calculate the maximum speed value

[0058] is the acceleration signal amplitude sequence of the kth particle at the pth iteration (k) (p) is calculated by integration to obtain the corresponding velocity amplitude sequence U (k) (p), which is an nt×1 column vector. Take the element with the largest absolute value in U (k) (p) as the maximum speed value of the current iteration in m / s.

[0059] C3, calculate the maximum displacement value

[0060] is the velocity sequence of the kth particle at the pth iteration (k) (p) is calculated by integration to obtain the corresponding displacement amplitude sequence X (k) (p), which is an nt×1 column vector, take the element with the largest absolute value in X (k) (p) as the maximum displacement value of the current iteration in meters.

[0061] C4, calculate the fitness

[0062] is the fitness value of the kth particle at the pth iteration is calculated by the following formula:

[0063]

[0064] wherein, parameter θmax , U max , X max are positive numbers, parameter θ max is in m / s 2 , U max is in m / s, X max is in meters, parameters α, β, γ are also positive numbers and satisfy α + β + γ = 1, parameter θ max , U max , X max , α, β, γ are set by engineers on site.

[0065] D, calculate the global optimal position

[0066] The parameters of the kth particle in the pth iteration The calculation formula is:

[0067]

[0068] The matrix PB (k) (p) The calculation formula is:

[0069]

[0070] Parameter Y g (p) The calculation formula is:

[0071]

[0072] The calculation formula of the global optimal position GB(p) in the pth iteration is:

[0073]

[0074] E, judge whether the global optimal value Y g (p) is less than the set threshold ε, that is, Y g (p) < ε; or whether the iteration number p is greater than the maximum iteration number T max , that is, p > T max . If satisfied, go to step G; otherwise, let p = p + 1 and go to step F.

[0075] F, optimize the wavelet group signal parameters F1, calculate the weight

[0076] The calculation formula of the weight w(p) in the pth iteration is:

[0077]

[0078] Wherein, parameters w max , w min take values in the range of 0-1, and w max > w min ; parameter wmax , w min Set by engineers on site.

[0079] F2, calculate learning factors

[0080] At the pth iteration, the calculation formulas of learning factors η1(p) and η2(p) are respectively:

[0081]

[0082]

[0083] wherein, parameter η max , η min is in the range of 0~3, and η max > η min , parameter η max , η min is set by engineers on site.

[0084] F3, calculate the motion step matrix and position matrix of particles

[0085] At the pth iteration, the calculation formula of the motion step matrix of the kth particle is:

[0086] V (k) (p) = w(p)V (k) (p-1) + q1η1(p)(PB (k) (p-1) - Z (k) (p-1))

[0087] + q2η2(p)(GB(p-1) - Z (k) (p-1))

[0088] wherein, parameter q1, q2 is in the range of 0~1, and is set by engineers on site.

[0089] Assign all elements of the 3rd column of matrix V (k) (p) to 0, and calculate the position matrix of the kth particle at the pth iteration according to the following formula:

[0090] Z (k) (p) = Z (k) (p-1) + V (k) (p)

[0091] F4, calculate the half-sine parameter matrix

[0092] The calculation formulas of matrix Δ (k) (p), PBN (k) (p) and GBN(p) are respectively:

[0093]

[0094] matrix VN (k) (p) = w(p) VN

[0095] VN (k) (p) = w(p) VN (k) (p-1) + q1η1(p)(PBN (k) (p-1) - Δ (k) (p-1)

[0096] + q2η2(p)(GBN(p-1) - Δ (k) (p-1)

[0097] matrix VN (k) (p) is calculated from matrix VN (k) (p), and the calculation formula is:

[0098]

[0099] F5, recalculate the particle position matrix

[0100] At the pth iteration, the kth particle position matrix Z (k) (p) is the element in the 3rd column and jth row The calculation formula is:

[0101]

[0102] A new particle position matrix Z (k) (p) is obtained, and the process goes to step B.

[0103] G, calculate the time-domain acceleration shock signal

[0104] The time-domain acceleration shock signal M f is calculated from matrix GBN(p), and the calculation formula is:

[0105]

[0106] Compared with the prior art, the present application has the following beneficial effects:

[0107] 1. In the process of generating a time-domain acceleration shock signal from a shock response spectrum, the present application does not need to round the half-sine parameters in the optimization process, avoiding the precision loss caused by rounding.

[0108] 2. In the process of generating a time-domain acceleration shock signal from a shock response spectrum, the present application does not depend on actual experimental data, and can generate a time-domain acceleration shock signal from a shock response spectrum without experimental data support.

[0109] 3. The application introduces displacement and velocity constraints in the process of generating time-domain acceleration shock signals from shock response spectrum, so that the generated acceleration shock signals can not only accurately fit the target shock response spectrum, but also effectively control the maximum displacement and maximum velocity of the shock signal, ensuring the reliability and safety of the shock vibration experiment. BRIEF DESCRIPTION OF DRAWINGS

[0110] Figure 1 is a flowchart of the application. DETAILED DESCRIPTION

[0111] The application will be further described below in conjunction with the drawings. According to the flowchart shown, the embodiments of the application are as follows: Figure 1

[0112] Take the target shock response spectrum of 10Hz-100Hz as an example, where 10Hz-20Hz is the rising spectrum, and the amplitude increases from 8m / s 2 to 18m / s 2 ; 20Hz-100Hz is a flat spectrum, and the amplitude remains 18m / s 2 unchanged. When the time-domain acceleration shock signal is generated by using the traditional algorithm, the error spectrum can be controlled within 6.5m / s 2 , the maximum velocity of the shock signal reaches 3.71×10 -1 m / s, and the maximum displacement reaches 4.61×10 -3 m. When the time-domain acceleration shock signal is generated by using the method of the application, the error spectrum can be controlled within 4.49m / s 2 , and the maximum velocity of the shock signal can be controlled within 3.38×10 -1 m / s, and the maximum displacement can be controlled within 2.75×10 -3 m. By using the method of the application, the maximum displacement and maximum velocity of the shock signal are effectively controlled, and the spectrum error is significantly reduced.

[0113] The application is not limited to the embodiments, and any equivalent concept or change within the technical scope disclosed by the application is included in the protection scope of the application.​

Claims

1. A method for generating time-domain acceleration shock signals from shock response spectrum based on wavelet superposition, characterized by: The following steps are involved: A. Initialization A1. Setting the target shock response spectrum Define the target shock response spectrum as R t , R t For an n×2 matrix, the expression is: R t =[H,L] Where H represents the frequency vector of the target impact response spectrum, H=[h1,h2,...,h n ] T is an n×1 column vector with the unit of Hz; L represents the amplitude vector of the target impact response spectrum, L=[l1,l2,...,l n ] T , is an n×1 column vector with unit of m / s 2 ; n is the number of points in the target spectrum, the superscript T represents the vector transpose; the parameters h1,h2,...,h n ,l1,l2,...,l n Set up on site by engineers; A2. Define the parameters of the wavelet superposition optimization algorithm A2-1. Set basic parameters The basic parameters include the number of particles N0, the maximum number of iterations T max and threshold ε, where N0, T max are all integers greater than 100; ε is in the range of 0 to 1. The basic parameters are all set on-site by engineers; A2-2. Setting the global optimal position Define the global optimal position matrix GB(p) at the pth iteration. GB(p) is an n×3 matrix. The specific expression is: Among them, the element gb in the i-th column and j-th row i,j (p) When p=0, the initial value is 0, that is, gb i,j (0) = 0; p is the number of iterations, i = 1, 2, 3, j = 1, 2, ..., n, p = 1, 2, ..., T max ; Define the matrix PB (k) (p), PB (k) (p) is an n×3 matrix, and the specific expression is: Among them, the element in the i-th column and j-th row When p=0, the initial value is 0, that is, k is the number of particles, k = 1, 2, ..., N0; When p = 0, define the parameter and Y g The initial values ​​of (p) are 0 and 1, that is, A2-3. Setting the half-sine parameter matrix Define the half-sine parameter value matrix of the kth particle in the pth iteration as Δ (k) (p), the k-th particle motion trend matrix of the p-th iteration is VN (k) (p), the individual memory matrix of the kth particle in the pth iteration is PBN (k) (p), the p-th iteration global kurtosis matrix GBN(p); the matrix Δ (k) (p), VN (k) (p), PBN (k) (p) and GBN(p) are both n×L M The matrix, matrix Δ (k) (p), PBN (k) (p) and GBN(p) elements are 0 or 1, and the parameter L M The calculation formula is: Among them, the parameter L min 、L max Are all positive odd numbers and must satisfy L max >L min , L min ≥3, set on-site by engineers; Matrix Δ (k) The expression of (p) is: Among them, the element in the oth column and jth row When p=0, the initial value is 0, that is, Matrix VN (k) The expression of (p) is: Among them, the element in the oth column and jth row Initial value when p = 0 For the interval [VM min ,VM max ] Random number within; parameter VM min The value range is -6 to -4, parameter VM max The value range is 4 to 6 and is set on site by engineers; Matrix PBN (k) The expression of (p) is: Among them, the element in the oth column and jth row When p=0, the initial value is 0, that is, The expression of matrix GBN(p) is: Among them, the element gbn in the oth column and jth row o,j (p) When p=0, the initial value is 0, that is, gbn o,j (0)=0; A2-4. Setting the particle position matrix Define the position matrix of the kth particle in the pth iteration as Z (k) (p), Z (k) (p) is an n×3 matrix, and the specific expression is: Among them, when p = 0, the initial value of the element in the first column and the jth row is For the interval [U min ,U max ]; the element in the 2nd column and the jth row For the interval [0,t m ) in the random number; the initial value of the element in the 3rd column and the jth row For the interval [L min ,L max ] is a random positive odd number; parameter U min 、U max are all positive numbers and satisfy U max >U min , t m is the duration of the acceleration signal in seconds, t m Is a positive number; parameter U max 、U min , t m Set up on site by engineers; A2-5. Setting the particle motion step matrix Define the k-th particle motion step matrix V in the p-th iteration (k) (p), V (k) (p) is an n×3 matrix, and the expression is: Among them, when p = 0, the initial value of the element in the first column and the jth row is For the interval [VA min ,VA max ]; the initial value of the element in the 2nd column and the jth row For the interval [Vt min ,Vt max ]; initial value of the third column element is 0; parameter VA min 、Vt min is a negative number, VA max 、Vt max It is a positive number and is set on site by the engineer; B. Calculate the wavelet group signal At the pth iteration, the kth particle position matrix Z (k) (p) Calculate the acceleration signal amplitude using the following formula: Among them, M (k) (p) is the amplitude sequence of the acceleration signal of the kth particle in the pth iteration, in m / s 2 , is a column vector of nt×1; h j Set the jth element of the target spectrum frequency vector H for step A-1; k q ∈{1,2,...,nt}, the calculation formula of the parameter nt is: T s The sampling time is in seconds and is set on site by the engineer. <·> indicates rounding up. C. Calculate the fitness value of each particle C1. Calculate spectral error The kth particle acceleration signal amplitude sequence M at the pth iteration (k) (p) Calculate the amplitude of the jth frequency point of the target spectrum according to the shock response spectrum analysis method in the national standard GB / T 29716.4-2018 "Mechanical vibration and shock signal processing Part 4: Shock response spectrum analysis" The unit is m / s 2 , calculate the error spectrum θ according to the following formula (k) (p), unit is m / s 2 : Among them, l j Set the j-th element of the target spectrum amplitude vector L for step A-1; |·| represents taking the absolute value; C2. Calculate the maximum speed value The kth particle acceleration signal amplitude sequence M at the pth iteration (k) (p) The corresponding velocity amplitude sequence U is obtained by integral calculation (k) (p), is a column vector of nt×1; take U (k) The element with the largest absolute value in (p) is used as the maximum speed value of the current iteration The unit is m / s; C3. Calculate the maximum displacement value The velocity sequence U of the kth particle at the pth iteration is (k) (p) The corresponding displacement amplitude sequence X is obtained by integral calculation (k) (p), is a column vector of nt×1, and takes X (k) The element with the largest absolute value in (p) is used as the maximum displacement value of the current iteration The unit is meter; C4. Calculate fitness The fitness value of the kth particle at the pth iteration is Calculated by the following formula: Among them, the parameter θ max 、U max 、X max are all positive numbers, parameter θ max The unit is m / s 2 , U max The unit is m / s, X max The unit is meter. The parameters α, β, and γ are also positive and satisfy α+β+γ=1. The parameter θ max 、U max 、X max , α, β, and γ are set on site by engineers; D. Calculate the global optimal position Parameters of the kth particle at the pth iteration The calculation formula is: Matrix PB (k) (p) is calculated as follows: Parameter Y g (p) is calculated as follows: The calculation formula for the p-th iteration global optimal position GB(p) is: E. Determine the global optimal value Y g (p) is less than the set threshold ε, that is, Y g (p)<ε; or whether the number of iterations p is greater than the maximum number of iterations T max , that is, p>T max ; If satisfied, go to step G; otherwise set p = p + 1 and go to step F; F. Optimize wavelet group signal parameters F1. Calculate weights At the pth iteration, the calculation formula of weight w(p) is: Among them, the parameter w max 、w min The value ranges from 0 to 1, and w max >w min ; parameter w max 、w min Set up on site by engineers; F2. Calculate the learning factor At the pth iteration, the calculation formulas of the learning factors η1(p) and η2(p) are: Among them, the parameter η max ,η min The value ranges from 0 to 3, and η max >η min , parameter η max ,η min Set up on site by engineers; F3. Calculate the particle's motion step matrix and position matrix At the pth iteration, the calculation formula of the motion step matrix of the kth particle is: V (k) (p)=w(p)V (k) (p-1)+q1η1(p)(PB (k) (p-1)-Z (k) (p-1))+q2η2(p)(GB(p-1)-Z (k) (p-1)) Among them, the parameters q1 and q2 have values ​​in the range of 0 to 1 and are set on site by engineers; The matrix V (k) (p) All elements in the third column are assigned 0, and the position matrix of the kth particle at the pth iteration is calculated according to the following formula: Z (k) (p)=Z (k) (p-1)+V (k) (p) F4. Calculate half-sine parameter matrix Matrix Δ (k) (p), PBN (k) The calculation formulas for (p) and GBN(p) are: Matrix VN (k) (p) is calculated as follows: VN (k) (p)=w(p)VN (k) (p-1)+q1η1(p)(PBN (k) (p-1)-Δ (k) (p-1))+q2η2(p)(GBN(p-1)-Δ (k) (p-1)) By matrix VN (k) (p) Calculate the matrix Δ (k) (p), calculated as: F5, recalculate particle position matrix At the pth iteration, the kth particle position matrix Z (k) (p) The element in the 3rd column and the jth row The calculation formula is: Get the new particle position matrix Z (k) (p), go to step B; G. Calculate the time domain acceleration impact signal The time domain acceleration impact signal M is calculated by the matrix GB(p) f , the calculation formula is:

Images