Hemispherical harmonic oscillator damping characteristic online identification method realized by active standing wave precession

Abstract

The invention discloses an online identification method for damping characteristics of a hemispherical harmonic oscillator by utilizing active standing wave precession, and belongs to the technical field of inertia. According to the invention, the problem that the damping characteristic of the hemispherical resonator in the working process of the hemispherical resonator gyroscope cannot be identified online in real time at present is solved. According to the method, the extra driving force applied by the control electrode of the hemispherical resonator gyroscope is used for driving the standing wave of the harmonic oscillator to actively precesde, and the damping characteristic of the harmonic oscillator can be detected in real time in the working process of the hemispherical resonator gyroscope system by resolving the control quantity of each part. The error of the damping characteristic of the hemispherical resonator gyroscope system caused by the change of the external environment and the service life of components is overcome, and the problem that the damping characteristic of the gyroscope is difficult to accurately measure during long-time work is solved. Compared with a traditional off-line measurement and identification method using a rotary table when leaving a factory, the method can stably measure the damping characteristics for a long time. The method can be applied to online identification of the damping characteristics of the hemispherical harmonic oscillator.

Description

Technical Field

[0001] The invention belongs to the field of inertial technology, and in particular relates to an online identification method of a hemispherical resonator damping characteristic realized by using active standing wave precession. Background Art

[0002] The hemispherical resonant gyroscope is a rate integrating gyroscope based on the Coriolis effect. It senses the angular velocity of motion by detecting the precession effect of the hemispherical resonator's vibration standing wave relative to the base. It has the advantages of simple structure, high reliability and wide detection range.

[0003] Since the position of the vibration standing wave in the rate integration mode is not fixed, it has a higher requirement for the circumferential isotropy of the system. However, during the processing and assembly of the hemispherical resonator, it is inevitable that the uneven distribution of damping and other characteristics will occur, affecting the symmetry of the gyro structure and causing measurement errors. Therefore, identifying the damping characteristics of the hemispherical resonator is an important means to further eliminate the damping error and improve the accuracy of the hemispherical resonator gyroscope.

[0004] However, the identification of the damping characteristics of the hemispherical resonant gyroscope is currently still at the stage of turntable experimental measurement, that is, the damping characteristics can only be identified before the hemispherical resonant gyroscope starts working, and it is still impossible to perform real-time online identification of the damping characteristics of the hemispherical resonator during the operation of the hemispherical resonator. Therefore, it is an urgent problem to propose a real-time online identification method for the damping characteristics of the hemispherical resonator. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that it is still impossible to perform real-time online identification of the damping characteristics of a hemispherical resonator during the operation of a hemispherical resonator gyroscope, and to propose an online identification method for the damping characteristics of a hemispherical resonator using active standing wave precession.

[0006] The technical solution adopted by the present invention to solve the above technical problems is: an online identification method of the damping characteristics of a hemispherical resonator using active standing wave precession, the method specifically comprising the following steps:

[0007] Step 1: Establish equations for the first-order derivative of the standing wave energy and the first-order derivative of the orthogonal error of the hemispherical resonant gyroscope, use the rate integral control mode to perform closed-loop control on the hemispherical resonant gyroscope, and obtain the orthogonal control force required by the resonator when the hemispherical resonant gyroscope is in a stable state;

[0008] Step 2: Introduce active driving force f v Drive the standing wave to precess, using the spindle amplitude control quantity d a , orthogonal control quantity d q And the external drive control quantity d v The driving electrode of the hemispherical resonant gyroscope is superimposed to drive the hemispherical resonant gyroscope to vibrate;

[0009] At each data sampling moment during the vibration of the hemispherical resonant gyroscope, the combined signal S, the combined signal R, the main axis amplitude control quantity and the orthogonal control quantity are collected, and the standing wave angle at each sampling moment is calculated based on S and R;

[0010] Step three: Identify the damping distribution characteristics of the hemispherical resonant gyroscope based on the orthogonal control force required by the resonator in step one and the data collected in step two.

[0011] The beneficial effects of the present invention are:

[0012] The method of the present invention utilizes the additional driving force applied by the control electrode of the hemispherical resonant gyroscope to drive the standing wave of the resonator to actively precess, and by solving the control quantity of each part, the damping characteristics of the resonator can be detected in real time during the operation of the hemispherical resonant gyroscope system. It overcomes the errors caused by changes in the external environment and the life of components to the damping characteristics of the hemispherical resonant gyroscope system, and solves the problem that the gyroscope damping characteristics are difficult to accurately measure during long-term operation. Compared with the traditional off-line measurement and identification method using a turntable at the time of leaving the factory, the method of the present invention can stably measure the damping characteristics for a long time. And it can be realized by only requiring a matching control circuit, and has a wider range of applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 It is the equivalent model of hemispherical oscillator mechanics;

[0014] In the figure, represents the intrinsic damping axis x of the hemispherical resonator τ The force on represents the intrinsic damping axis y of the hemispherical resonator τ The force on

[0015] Figure 2 This is the signal flow diagram of the hemispherical resonant gyroscope control loop. DETAILED DESCRIPTION

[0016] Specific implementation method 1: This implementation method describes an online identification method of the damping characteristics of a hemispherical resonator using active standing wave precession, and the method specifically includes the following steps:

[0017] Step 1: Establish equations for the first-order derivative of the standing wave energy and the first-order derivative of the orthogonal error of the hemispherical resonant gyroscope, use the rate integral control mode to perform closed-loop control on the hemispherical resonant gyroscope, and obtain the orthogonal control force required by the resonator when the hemispherical resonant gyroscope is in a stable state;

[0018] Step 2: Introduce active driving force f v Drive the standing wave to precess, using the spindle amplitude control quantity d a , orthogonal control quantity d q And the external drive control quantity dv The driving electrode of the hemispherical resonant gyroscope is superimposed to drive the hemispherical resonant gyroscope to vibrate;

[0019] At each data sampling moment during the vibration of the hemispherical resonant gyro, the combined signal S, the combined signal R, the main axis amplitude control quantity and the orthogonal control quantity are collected, and the standing wave angle at each sampling moment is calculated based on S and R;

[0020] Step three: Identify the damping distribution characteristics of the hemispherical resonant gyroscope based on the orthogonal control force required by the resonator in step one and the data collected in step two.

[0021] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the equations of the first-order derivative of the standing wave energy and the first-order derivative of the orthogonal error of the hemispherical resonant gyroscope are respectively:

[0022]

[0023]

[0024] Where E represents the standing wave energy, Q represents the orthogonality error, represents the damping mismatch, represents the damping parameter, θ represents the standing wave rotation angle, a represents the main standing wave amplitude, q represents the orthogonal wave amplitude, ω represents the resonator vibration frequency, and f represents the standing wave amplitude. qs represents the orthogonal control force, f ac represents the amplitude control force, θ τ represents the resonator damping axis angle, θ ω represents the rigid axis angle of the resonator, and Δω is the frequency decomposition.

[0025] The other steps and parameters are the same as those in the first embodiment.

[0026] In this field, and is defined as follows:

[0027]

[0028] Among them, τ 1 is the main standing wave axis damping time constant, τ 2 is the orthogonal axis damping time constant.

[0029] Specific implementation method three: This implementation method is different from specific implementation method one or two in that the rate integral control mode is used to perform closed-loop control on the hemispherical resonant gyroscope to obtain the orthogonal control force required by the resonator in the stable state of the hemispherical resonant gyroscope; specifically:

[0030] Through the amplitude and orthogonal control loop, the rate integral control mode is used to perform closed-loop control on the hemispherical resonant gyro. When the hemispherical resonant gyro is in a stable state, E = a 2 , Q=0、q=0,according to the first order of orthogonal error

[0031] Derivative In the form of, the orthogonal control force required for the resonator is obtained as:

[0032] f qs =aωΔωsin4(θ-θ ω )(9)

[0033] In formula (9), ωΔω is defined as formula (10);

[0034]

[0035] Among them, ω x is the vibration frequency on the principal axis x, ω y is the vibration frequency on the principal coordinate axis y.

[0036] The other steps and parameters are the same as those in the first or second embodiment.

[0037] Specific implementation method 4: This implementation method is different from any one of the specific implementation methods 1 to 3 in that, at each data sampling moment during the vibration process of the hemispherical resonant gyroscope, the combined signal S, the combined signal R, the main axis amplitude control amount and the orthogonal control amount are collected, and the standing wave angle at each sampling moment is calculated based on S and R; the specific process is:

[0038] Step 21: Initialize sampling time i=0;

[0039] Step 22: Set the i-th sampling time t i The standing wave rotation angle is denoted as θ(i), and t i The speed at the time is recorded as ω(i), then the i+1th sampling time t i+1 The desired standing wave rotation angle is θ(i+1)=θ(i)+Δθ, where Δθ is the preset angle increment. The i+1th sampling time t i+1 =t i +Δθ / ω(i);

[0040] Step 23: At the i+1th sampling time t i+1 For S(i+1), R(i+1), spindle amplitude control value d a (i+1) and the orthogonal control quantity d q (i+1) collects data;

[0041] Calculate the standing wave rotation angle θ(i+1) at the i+1th sampling moment according to S(i+1) and R(i+1);

[0042] Step 24: Determine whether i<N-1 is satisfied:

[0043] If not satisfied, set i=i+1 and return to step 22;

[0044] If satisfied, execute step 25;

[0045] Step 25: Store the data collected at each sampling moment in the form of an array matrix in formula (13) in the FPGA;

[0046]

[0047] Among them, D a is an array matrix composed of the spindle amplitude control data at each sampling moment, D q It is an array matrix composed of orthogonal control quantity data at each sampling moment, and H is a coefficient array matrix related to the standing wave angle data at each sampling moment.

[0048] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0049] Specific implementation method 5: This implementation method is different from the specific implementation methods 1 to 4 in that the specific process of step 3 is as follows:

[0050] Step 31: Orthogonal control force applied to the resonator and spindle control force They are:

[0051]

[0052] in, is the phase delay error, f a is the spindle amplitude control force, f q is the orthogonal control force;

[0053] In the orthogonal control loop With f qs If they are equal, the actual orthogonal control force output by the orthogonal control loop is:

[0054]

[0055] Transform equation (15) into the orthogonal control quantity equation (16):

[0056]

[0057] Among them, G cis the gain from the hemispherical resonant gyroscope driving signal to the electrostatic excitation force of the resonator;

[0058] Then the orthogonal control quantity equation of formula (16) is rewritten as:

[0059] d q =b 1 +p 1 sin4θ+q 1 cos4θ (17)

[0060] Among them, b 1 、p 1 and q 1 are the coefficients of the orthogonal control quantity equations;

[0061] Let X 1 =[p 1 q 1 b 1 ] T , D q =HX 1 ,but

[0062] X 1 =(H T H) -1 H T D q (18)

[0063] The coefficient of the sinusoidal component in equation (16) is obtained as:

[0064]

[0065] And let the parameter G be:

[0066]

[0067] Step 32: In the amplitude control loop With f ac If they are equal,

[0068]

[0069] Will Substitute into equation (7) under the steady state:

[0070]

[0071] According to the identification parameter G, the terms in equation (22) are transposed and transformed into equation (23):

[0072]

[0073] The control quantity equation of formula (23) is rewritten as:

[0074]

[0075] Let X 2 =[p 2 q 2 b 2 ] T , but

[0076]

[0077] According to formula (23), the damping distribution characteristic parameters and X 2 The relationship between the parameters in

[0078] Step 3: According to the data collected in step 2, 1 =[p 1 q 1 b 1 ] T and X 2 =[p 2 q 2 b 2 ] T Perform joint identification to obtain the damping distribution characteristic parameters and θ τ .

[0079] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0080] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the damping distribution characteristic parameter is 2 The relationship between the parameters is as follows:

[0081] Damping parameters:

[0082]

[0083] Damping Mismatch:

[0084]

[0085] Resonator damping axis angle:

[0086]

[0087] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0088] Specific implementation method 7: This implementation method is different from specific implementation methods 1 to 6 in that the specific process of step 33 is as follows:

[0089] Step 331: Initialize the forgetting factor λ ∈ (0, 1), and initialize the parameter vector estimate X 1 (0) =

[000] T 、X 2 (0) = [0 0 0] T Initialize the covariance matrix P 1 (0) = λI, P 2 (0) = λI, where I is the identity matrix;

[0090] And initialize the iteration number j = 1;

[0091] Step 332: Calculate the gain vector k 1 (j):

[0092] k 1 (j) = P 1 (j - 1)·h T (j) / (λ + h(j)·P 1 (j - 1)·h T (j)) (29)

[0093] where h(j) is the j-th row data of the data matrix H, and h T (j) is the transpose of h(j);

[0094] Update the parameter estimate X 1 (j), X 1 = [p 1 (j)q 1 (j)b 1 (j)] T :

[0095] X 1 (j) = X 1 (j - 1) + k 1 (j)·(d q (j) - h(j)·X 1 (j - 1)) (30)

[0096] where d q (j) is the j-th row data of the data matrix D q ;

[0097] Step 333: Update the covariance matrix P 1 (j):

[0098] P 1 (j) = (P 1 (j - 1) - k 1 (j)·h(j)·P 1 (j - 1)) / λ (31)

[0099] Step 334: Update parameter G(j):

[0100]

[0101] Step 335: Calculate the gain vector k 2 (j):

[0102] k 2 (j) = P 2 (j-1) h T (j) / (λ+h(j)·P 2 (j-1) h T (j)) (33)

[0103] Update parameter estimates X 2 (j), X 2 (j) = [p 2 (j)q 2 (j)b 2 (j)] T :

[0104] X 2 (j) = X 2 (j-1)+k 2 (j)·(G(j)·d a (j)-h(j)·X 2 (j-1)) (34)

[0105] Among them, d a (j) is the data matrix D a The j-th row of data;

[0106] Step 336: Update the covariance matrix P 2 (j):

[0107] P 2 (j)=(P 2 (j-1)-k 2 (j)·h(j)·P 2 (j-1)) / λ (35)

[0108] Step 337: Update damping parameters Damping Mismatch and the resonator damping axis rotation angle θ τ [j]:

[0109]

[0110] Step 338: Whether j=N is satisfied:

[0111] If j=N, execute step 339;

[0112] If j=N is not satisfied, set j=j+1 and return to step 332;

[0113] Step 339: Based on the calculated damping parameters at each sampling time Damping Mismatch and the resonator damping axis rotation angle θ τ [j] Obtain the damping distribution characteristic parameter identification results.

[0114] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0115] Specific implementation eight: This implementation differs from any one of specific implementations one to seven in that, in the damping distribution characteristic parameter, the identification result of the damping parameter is: The mean of .

[0116] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0117] Specific implementation method 9: This implementation method is different from any one of specific implementation methods 1 to 8 in that, in the damping distribution characteristic parameter, the identification result of the damping mismatch is: The mean of .

[0118] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0119] Specific implementation method 10: This implementation method is different from any one of specific implementation methods 1 to 9 in that, in the damping distribution characteristic parameter, the identification result of the resonator damping axis rotation angle is: the resonator damping axis rotation angle θ τ [j], the mean of j=1,2,…,N.

[0120] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0121] Example

[0122] The following is a further description of an online identification method of a hemispherical resonator damping characteristic using active standing wave precession according to the present invention in conjunction with the accompanying drawings:

[0123] The origin of the resonator modal excitation and detection coordinate system is defined as the center of the flat electrode. The x-axis of the resonator modal excitation and detection coordinate system coincides with the symmetry axis of the first electrode in the counterclockwise direction. The positive direction of the x-axis points from the origin to the first electrode. The y-axis is perpendicular to the x-axis, and the positive direction of the y-axis points from the origin to the third electrode in the counterclockwise direction. Figure 1 As shown in the figure, the coordinate axes of the harmonic oscillator modal excitation and detection coordinate system are the main coordinate axes x and y, and the non-ideal harmonic oscillator eigenfrequency axis x ω The angle between the x-axis of the harmonic oscillator modal excitation and detection coordinate system is θ ω, the intrinsic damping axis x of the non-ideal oscillator τ The angle between the x-axis of the harmonic oscillator modal excitation and detection coordinate system is θ τ , and θ ω is called the rigid axis angle of the resonator, and θ τ It is called the oscillator damping axis angle.

[0124] like Figure 2 The figure shows the control loop in the hemispherical resonant gyro control system composed of FPGA. After startup, the output signal of the hemispherical resonant gyro detection channel is demodulated using the sine and cosine reference signals generated by the phase-locked loop circuit, specifically:

[0125] The electrical signal sensed by the detection electrode is:

[0126]

[0127] Among them, a is the main standing wave amplitude, q is the orthogonal wave amplitude, ω is the resonator vibration frequency, θ is the standing wave rotation angle, is the phase delay error.

[0128] The sine and cosine reference signals obtained using the phase-locked loop are:

[0129]

[0130] In the demodulation circuit, the reference signal is multiplied by the electrical signal of formula (1), and then the multiplication result is passed through a low-pass filter to extract the low-frequency signal of the gyroscope to obtain the demodulated low-frequency variable C x , C y , S x , S y :

[0131]

[0132] Wherein, LPF represents low-pass filter.

[0133] By combining and calculating the low-frequency variables in equation (3), we can obtain the typical parameters of the working state of the hemispherical resonant gyroscope, namely, the standing wave energy E and the orthogonal error Q:

[0134]

[0135] According to the parameters in formula (4), the phase delay error can be calculated respectively: And the standing wave angle θ:

[0136]

[0137] The online identification method of the damping characteristics of a hemispherical resonator is described in detail according to the following steps.

[0138] Step 1: By using the Lynch mean method, the first-order derivative of the standing wave energy of the hemispherical resonant gyroscope can be obtained: and the first derivative of the quadrature error The equation form is as follows:

[0139]

[0140] Where E represents the standing wave energy, Q represents the orthogonality error, represents the damping mismatch, represents the damping parameter, θ represents the standing wave rotation angle, a represents the main standing wave amplitude, q represents the orthogonal wave amplitude, ω represents the resonator vibration frequency, and f represents the standing wave amplitude. qs represents the orthogonal control force, f ac represents the amplitude control force, θ τ represents the resonator damping axis angle, θ ω represents the rigid axis angle of the resonator, Δω is the frequency splitting;

[0141] Through the amplitude and orthogonal control loop, the rate integral control mode is used to perform closed-loop control on the hemispherical resonant gyro. When the hemispherical resonant gyro is in a stable state, E = a 2 , Q=0、q=0,according to the first order of orthogonal error

[0142] Derivative In the form of, that is, considering the orthogonal error of formula (8), the orthogonal control force required by the resonator is obtained as:

[0143] f qs =aωΔωsin4(θ-θ ω ) (9)

[0144] In formula (9), ωΔω is defined as formula (10), which represents the frequency difference between the main axis of vibration of the resonator and the orthogonal axis, that is, the frequency splitting;

[0145]

[0146] Among them, ω x is the vibration frequency on the principal axis x, ω y is the vibration frequency on the principal axis y;

[0147] Step 2: Based on amplitude and orthogonal control, introduce a certain active driving force f v Drive the standing wave to precess, using the spindle amplitude control quantity d a , orthogonal control quantity d q And the external drive control quantity d v The driving electrode of the hemispherical resonant gyroscope is superimposed to drive the hemispherical resonant gyroscope to vibrate;

[0148] The control loop control force and the external active driving force act on the resonator through the excitation electrode. The gain of the hemispherical resonant gyroscope drive signal to the electrostatic excitation force of the resonator is G c , spindle amplitude control value d a , orthogonal control quantity d q And the external drive control quantity d v The relationship with the control quantity is:

[0149]

[0150] Among them, f a is the spindle amplitude control force, f q is the orthogonal control force, f v It is an active driving force;

[0151] according to Figure 1 The coordinate system shown will control the spindle control force f output by the control loop a , orthogonal control force f q , Active driving force f v Decomposed into the hemispherical resonator gyroscope base coordinate system XOY, that is, according to formula (12), it is converted into the force of the hemispherical resonator on the two coordinate axes:

[0152]

[0153] Among them, the phase delay error It is generated during circuit signal processing and can be regarded as a constant.

[0154] At each data sampling moment during the vibration of the hemispherical resonant gyro, the combined signal S, the combined signal R, the main axis amplitude control quantity and the orthogonal control quantity are collected, and the standing wave angle at each sampling moment is calculated based on S and R;

[0155] Since the spatial position of the resonator standing wave changes with time, in order to ensure that data acquisition covers as much standing wave angle information as possible, the following method is used to determine the sampling time and position:

[0156] The sampling start time is t 0 , t 0 The standing wave rotation angle at time is denoted as θ(0), and t 0 The rotation speed at the moment is recorded as ω(0), and the standing wave angle expected to be collected at the next sampling moment is θ(1)=θ(0)+Δθ, where Δθ is the preset angle increment. Then the next sampling moment t 1 =t 0 +Δθ / ω(0), and so on. That is, the sampling time and speed of the previous sample point are used to dynamically determine the sampling time of the next sample point. The above sampling process is specifically expressed as:

[0157] Step 21: Initialize sampling time i=0;

[0158] Step 22: Set the i-th sampling time t i The standing wave rotation angle is denoted as θ(i), and t i The speed at the time is recorded as ω(i), then the i+1th sampling time t i+1 The desired standing wave rotation angle is θ(i+1)=θ(i)+Δθ, where Δθ is the preset angle increment. The i+1th sampling time t i+1 =t i +Δθ / ω(i);

[0159] Step 23: At the i+1th sampling time t i+1 For S(i+1), R(i+1), spindle amplitude control value d a (i+1) and the orthogonal control quantity d q (i+1) collects data;

[0160] Calculate the standing wave rotation angle θ(i+1) at the i+1th sampling moment according to S(i+1) and R(i+1);

[0161] Step 24: Determine whether i<N-1 is satisfied:

[0162] If not satisfied, set i=i+1 and return to step 22;

[0163] If satisfied, execute step 25;

[0164] Step 25: Store the data collected at each sampling moment in the form of an array matrix in formula (13) in the FPGA;

[0165] Due to the need to solve the sinusoidal form equation in the subsequent steps, the standing wave angle θ is processed, and sin4θ(i) and cos4θ(i) are calculated and then stored;

[0166] The obtained data points are written in the form of a matrix,

[0167]

[0168] Among them, D a is an array matrix composed of the spindle amplitude control data at each sampling moment, D q is an array matrix composed of orthogonal control quantity data at each sampling moment, and H is a coefficient array matrix related to the standing wave angle data at each sampling moment. The array matrix D shown in formula (13) q and H will be used in the subsequent least squares solution of the damping distribution characteristics.

[0169] The data collection algorithm is shown in Table 1:

[0170] Table 1 Data dynamic sampling algorithm

[0171]

[0172] Step 3: Identify the damping distribution characteristics of the hemispherical resonant gyroscope based on the data collected in step 2.

[0173] The specific process of step three is:

[0174] Step 31: Based on the closed-loop control, use a certain external driving control force f v Drive the standing wave to precess, and obtain the orthogonal control force applied to the resonator at this time and spindle control force They are:

[0175]

[0176] in, is the phase delay error, f a is the spindle amplitude control force, f q is the orthogonal control force;

[0177] In the orthogonal control loop With f qs If they are equal, the actual orthogonal control force output by the orthogonal control loop is:

[0178]

[0179] Transform equation (15) into the orthogonal control quantity equation (16):

[0180]

[0181] Among them, G c is the gain from the hemispherical resonant gyroscope driving signal to the electrostatic excitation force of the resonator;

[0182] This equation is equivalent to the superposition of a DC component and a sinusoidal component in form, so the orthogonal control quantity equation of equation (16) is rewritten as:

[0183] d q =b 1 +p 1 sin4θ+q 1 cos4θ (17)

[0184] Among them, b 1 、p 1 and q 1 are the coefficients of the orthogonal control quantity equations;

[0185] Let X 1 =[p1 q 1 b 1 ] T , apply the least squares formula, D q =HX 1 ,but

[0186] X 1 =(H T H) -1 H T D q (18)

[0187] Returning to the original control quantity equation, the coefficient of the sinusoidal component in equation (16) is obtained as:

[0188]

[0189] Since the frequency splitting Δω depends only on the processing technology of the resonator, it basically does not change in the hemispherical resonant gyroscope system and can be measured and eliminated when the hemispherical resonant gyroscope leaves the factory. Let the parameter G be:

[0190]

[0191] Step 32: In the amplitude control loop With f ac If they are equal,

[0192]

[0193] Will Substitute into the stable state equation (7):

[0194]

[0195] According to the identification parameter G, equation (22) is shifted and transformed into equation (23), and the control quantity expression of the spindle control is obtained:

[0196]

[0197] This equation is also a superposition of a DC component and a sinusoidal component in form, so the control quantity equation of equation (23) can be rewritten as:

[0198]

[0199] Let X 2 =[p 2 q 2 b 2 ] T , using the least squares formula, we get but

[0200]

[0201] According to formula (23), the damping distribution characteristic parameters and X 2 The relationship between the parameters is as follows:

[0202] Damping parameters:

[0203]

[0204] Damping Mismatch:

[0205]

[0206] Resonator damping axis angle:

[0207]

[0208] Step 3: According to the data collected in step 2, 1 =[p 1 q 1 b 1 ] T and X 2 =[p 2 q 2 b 2 ] T Perform joint identification to obtain the damping distribution characteristic parameters and θ τ .

[0209] Step 3: Initialize the forgetting factor λ∈(0,1) and the parameter vector estimate X 1 (0) =

[000] T , X 2 (0) =

[000] T , initialize the covariance matrix P 1 (0) = λI, P 2 (0) = λI, where I is the identity matrix;

[0210] And initialize the number of iterations j = 1;

[0211] Step 332: Calculate the gain vector k 1 (j):

[0212] k 1 (j) = P 1 (j-1) h T (j) / (λ+h(j)·P 1 (j-1) h T (j)) (29)

[0213] Among them, h(j) is the j-th row data of the data matrix H, hT (j) is the transpose of h(j);

[0214] Update parameter estimates X 1 (j), X 1 =[p 1 (j)q 1 (j)b 1 (j)] T :

[0215] X 1 (j) = X 1 (j-1)+k 1 (j)·(d q (j)-h(j)·X 1 (j-1)) (30)

[0216] Among them, d q (j) is the data matrix D q The j-th row of data;

[0217] Step 333: Update the covariance matrix P 1 (j):

[0218] P 1 (j)=(P 1 (j-1)-k 1 (j)·h(j)·P 1 (j-1)) / λ (31)

[0219] Step 334: Update parameter G(j):

[0220]

[0221] Step 335: Calculate the gain vector k 2 (j):

[0222] k 2 (j) = P 2 (j-1) h T (j) / (λ+h(j)·P 2 (j-1) h T (j)) (33)

[0223] Update parameter estimates X 2 (j), X 2 (j) = [p 2 (j) q 2 (j) b 2 (j)] T :

[0224] X 2 (j) = X 2(j-1)+k 2 (j)·(G(j)·d a (j)-h(j)·X 2 (j-1)) (34)

[0225] Among them, d a (j) is the data matrix D a The j-th row of data;

[0226] Step 336: Update the covariance matrix P 2 (j):

[0227] P 2 (j)=(P 2 (j-1)-k 2 (j)·h(j)·P 2 (j-1)) / λ (35)

[0228] Step 337: Update damping parameters Damping Mismatch and the resonator damping axis rotation angle θ τ [j]:

[0229]

[0230] Step 338: Whether j=N is satisfied:

[0231] If j=N, execute step 339;

[0232] If j=N is not satisfied, set j=j+1 and return to step 332;

[0233] Step 339: Calculate damping parameters The mean of the calculated mean is used as the damping parameter identified at the current moment; the damping mismatch is calculated The calculated mean is used as the damping mismatch identified at the current moment; the resonator damping axis rotation angle θ is calculated τ [j], j = the mean of 1, 2, ..., N, and the calculated mean is used as the resonator damping axis angle identified at the current moment.

[0234] The identification method of this embodiment can solve the problem of large amount of calculation of the traditional least square method. The identification process is shown in Table 2.

[0235] Table 2

[0236]

[0237]

[0238] The method of the present invention can identify the damping characteristic parameters online, solving the problem that the traditional method can only identify before starting work and cannot monitor the damping characteristic parameters during the work process. It should be noted that the method of the present invention does not require real-time identification of the damping characteristic parameters. In actual application, the damping characteristic parameters can be identified once at fixed intervals. When identification is required, N groups of data can be collected in advance.

[0239] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. An online identification method of the damping characteristics of a hemispherical resonator using active standing wave precession, characterized in that: The method specifically comprises the following steps: Step 1: Establish equations for the first-order derivative of the standing wave energy and the first-order derivative of the orthogonal error of the hemispherical resonant gyroscope, use the rate integral control mode to perform closed-loop control on the hemispherical resonant gyroscope, and obtain the orthogonal control force required by the resonator when the hemispherical resonant gyroscope is in a stable state; Step 2: Introduce active driving force f v Drive the standing wave to precess, using the spindle amplitude control quantity d a , orthogonal control quantity d q And the external drive control quantity d v The driving electrode of the hemispherical resonant gyroscope is superimposed to drive the hemispherical resonant gyroscope to vibrate; At each data sampling moment during the vibration of the hemispherical resonant gyro, the combined signal S, the combined signal R, the main axis amplitude control quantity and the orthogonal control quantity are collected, and the standing wave angle at each sampling moment is calculated based on S and R; Step three: Identify the damping distribution characteristics of the hemispherical resonant gyroscope based on the orthogonal control force required by the resonator in step one and the data collected in step two.

2. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 1, characterized in that: The equations of the first-order derivative of the standing wave energy and the first-order derivative of the orthogonal error of the hemispherical resonant gyroscope are respectively: Where E represents the standing wave energy, Q represents the orthogonality error, represents the damping mismatch, represents the damping parameter, θ represents the standing wave rotation angle, a represents the main standing wave amplitude, q represents the orthogonal wave amplitude, ω represents the resonator vibration frequency, and f represents the standing wave amplitude. qs represents the orthogonal control force, f ac represents the amplitude control force, θ τ represents the resonator damping axis angle, θ ω represents the rigid axis angle of the resonator, and Δω is the frequency decomposition.

3. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 2, characterized in that: The closed-loop control of the hemispherical resonant gyroscope is performed by adopting the rate integral control mode to obtain the orthogonal control force required by the resonator in the stable state of the hemispherical resonant gyroscope; specifically: Through the amplitude and orthogonal control loop, the rate integral control mode is used to perform closed-loop control on the hemispherical resonant gyro. When the hemispherical resonant gyro is in a stable state, E = a 2 , Q=0、q=0,according to the first-order derivative of the orthogonal error In the form of, the orthogonal control force required for the resonator is obtained as: f qs =aωΔωsin4(θ-θ ω ) (9) In formula (9), ωΔω is defined as formula (10); Among them, ω x is the vibration frequency on the principal axis x, ω y is the vibration frequency on the principal coordinate axis y.

4. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 3 is characterized in that: At each data sampling moment during the vibration of the hemispherical resonant gyroscope, the combined signal S, the combined signal R, the main axis amplitude control amount and the orthogonal control amount are collected, and the standing wave angle at each sampling moment is calculated according to S and R; the specific process is: Step 21: Initialize sampling time i=0; Step 22: Set the i-th sampling time t i The standing wave rotation angle is denoted as θ(i), and t i The speed at the time is recorded as ω(i), then the i+1th sampling time t i+1 The desired standing wave rotation angle is θ(i+1)=θ(i)+Δθ, where Δθ is the preset angle increment. The i+1th sampling time t i+1 =t i +Δθ / ω(i); Step 23: At the i+1th sampling time t i+1 For S(i+1), R(i+1), spindle amplitude control value d a (i+1) and the orthogonal control quantity d q (i+1) collects data; Calculate the standing wave rotation angle θ(i+1) at the i+1th sampling moment according to S(i+1) and R(i+1); Step 24: Determine whether i<N-1 is satisfied: If not satisfied, set i=i+1 and return to step 22; If satisfied, execute step 25; Step 25: Store the data collected at each sampling moment in the form of an array matrix in formula (13) in the FPGA; Among them, D a is an array matrix composed of the spindle amplitude control data at each sampling moment, D q It is an array matrix composed of orthogonal control quantity data at each sampling moment, and H is a coefficient array matrix related to the standing wave angle data at each sampling moment.

5. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 4 is characterized in that: The specific process of step three is: Step 31: Orthogonal control force applied to the resonator and spindle control force They are: in, is the phase delay error, f a is the spindle amplitude control force, f q is the orthogonal control force; In the orthogonal control loop With f qs If they are equal, the actual orthogonal control force output by the orthogonal control loop is: Transform equation (15) into the orthogonal control quantity equation (16): Among them, G c is the gain from the hemispherical resonant gyroscope driving signal to the electrostatic excitation force of the resonator; Then the orthogonal control quantity equation of formula (16) is rewritten as: d q =b1+p1sin4θ+q1cos4θ (17) Among them, b1, p1 and q1 are the coefficients of the orthogonal control quantity equation; Let X1 = [p1 q1 b1] T , D q =HX1, then X1=(H T H) -1 H T D q (18) The coefficient of the sinusoidal component in equation (16) is obtained as: And let the parameter G be: Step 32: In the amplitude control loop With f ac If they are equal, Will Substitute into the stable state equation (7): According to the identification parameter G, the terms in equation (22) are transposed and transformed into equation (23): The control quantity equation of formula (23) is rewritten as: Let X2 = [p2 q2 b2] T , but According to formula (23), the relationship between the damping distribution characteristic parameters and the parameters in X2 is obtained; Step 3: According to the data collected in step 2, X1 = [p1 q1 b1] T and X2 = [p2 q2 b2] T Perform joint identification to obtain the damping distribution characteristic parameters and θ τ .

6. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 5, characterized in that: The relationship between the damping distribution characteristic parameters and the parameters in X2 is as follows: Damping parameters: Damping Mismatch: Resonator damping axis angle:

7. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 6, characterized in that: The specific process of step 33 is as follows: Step 331: Initialize the forgetting factor λ∈(0,1) and the parameter vector estimate X1(0)=[000] T 、X2(0)=[000] T , initialize the covariance matrix P1(0)=λI, P2(0)=λI, where I is the identity matrix; And initialize the number of iterations j = 1; Step 332: Calculate the gain vector k1(j): k1(j)=P1(j-1)·h T (j) / (λ+h(j)·P1(j-1)·h T (j)) (29) Among them, h(j) is the j-th row data of the data matrix H, h T (j) is the transpose of h(j); Update parameter estimates X1(j), X1 = [p1(j)q1(j)b1(j)] T : X1(j)=X1(j-1)+k1(j)·(d q (j)-h(j)·X1(j-1)) (30) Among them, d q (j) is the data matrix D q The j-th row of data; Step 333, update the covariance matrix P1(j): P1(j)=(P1(j-1)-k1(j)·h(j)·P1(j-1)) / λ (31) Step 334: Update parameter G(j): Step 335: Calculate the gain vector k2(j): k2(j)=P2(j-1)·h T (j) / (λ+h(j)·P2(j-1)·h T (j)) (33) Update parameter estimate X2(j), X2(j) = [p2(j)q2(j)b2(j)] T : X2(j)=X2(j-1)+k2(j)·(G(j)·d a (j)-h(j)·X2(j-1)) (34) Among them, d a (j) is the data matrix D a The j-th row of data; Step 336: Update the covariance matrix P2(j): P2(j)=(P2(j-1)-k2(j)·h(j)·P2(j-1)) / λ (35) Step 337: Update damping parameters Damping Mismatch and the resonator damping axis rotation angle θ τ [j]: Step 338: Whether j=N is satisfied: If j=N, execute step 339; If j=N is not satisfied, set j=j+1 and return to step 332; Step 339: Based on the calculated damping parameters at each sampling time Damping Mismatch and the resonator damping axis rotation angle θ τ [j] Obtain the damping distribution characteristic parameter identification results.

8. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 7, characterized in that: In the damping distribution characteristic parameters, the identification result of the damping parameter is: The mean of j=1,2,…,N.

9. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 7, characterized in that: In the damping distribution characteristic parameters, the identification result of damping mismatch is: The mean of j=1,2,…,N.

10. The method for online identification of damping characteristics of a hemispherical resonator using active standing wave precession according to claim 7, characterized in that: Among the damping distribution characteristic parameters, the identification result of the resonator damping axis rotation angle is: τ [j], the mean of j=1,2,…,N.

Images