Inertial measurement array fault diagnosis method

Abstract

The invention relates to an inertial measurement array fault diagnosis method. An inertial measurement array is placed on a horizontal vibration isolation table, and signals output by three-axis gyroscopes and three-axis accelerometers of all inertial measurement units in the inertial measurement array are collected. The method comprises the following steps: firstly, performing feature extraction on an acquired signal by using discrete multivariate variational mode decomposition; secondly, constructing a generalized likelihood ratio test fault detection function, and performing fault detection on all gyroscopes and accelerometers; and finally, constructing a generalized likelihood ratio test fault isolation function to remove faults of the gyroscope and the accelerometer one by one. And after fault isolation is completed, outputting serial numbers of the gyroscope and the accelerometer with faults. The invention belongs to the technical field of inertial navigation and can be applied to an inertial measurement array system.

Description

Technical Field

[0001] The invention relates to an inertial measurement array fault diagnosis method, which is suitable for occasions of inertial measurement array fault diagnosis. Background Art

[0002] The inertial measurement array (IMU) concept, originating from the "array gyroscope" concept proposed by Bayard et al. in 2003 (Bayard DS. Combining multiple gyroscope outputs for increased accuracy[R]. Pasadena: NASA's Jet Propulsion Laboratory, 2003), combines multiple chip-level micro-electromechanical (MEMS) inertial measurement units (IMUs) into an array. Through data fusion, this significantly reduces random errors and improves navigation accuracy. MEMS IMUs are technically capable of mass production and have a large application scale in the consumer sector. They can be produced on a large scale and at low cost. Array technology holds promise for building low-cost, high-precision inertial navigation systems.

[0003] Fault diagnosis for inertial measurement arrays (IMAs) is currently underresearch. Therefore, research has drawn on the fault diagnosis of redundant IMAs. This approach typically uses the generalized likelihood ratio test (GLT), which uses the inconsistency between the residual parity vectors in the presence and absence of a fault as the basis for fault detection (Kevin C. Daly, Eliezer Gai, James V. Harrison. Generalized Likelihood Test for FDI in Redundant Sensor Configurations [J]. Guidance and Control, 1978.). During GLT fault diagnosis, the MEMSIMU output contains significant noise, outliers, and other errors, which severely impact IMA fault diagnosis and increase the false alarm rate. However, GLT fault diagnosis requires a relatively accurate dynamic model and error statistics to compensate for the residual parity vectors, which are difficult to obtain in practice.

[0004] Therefore, for the fault diagnosis of inertial measurement arrays, it is of great significance to remove the influence of noise and outliers on fault diagnosis, effectively isolate the faulty sensor, prevent the inertial measurement array from being destroyed during data fusion, and ensure the reliability of the inertial measurement array. Summary of the Invention

[0005] The technical problems solved by the present invention are:

[0006] When the inertial measurement array is used for fault diagnosis, the noise and outliers contained in the outputs of the gyroscope and accelerometer will increase the false alarm rate of the fault diagnosis and reduce the accuracy of the fault diagnosis.

[0007] The technical solution of the present invention:

[0008] Before detecting and isolating the output faults of the gyroscope and accelerometer, discrete multivariate mode decomposition is used to extract features from the output signals. This effectively solves the problem that the noise and outliers contained in the gyroscope and accelerometer outputs will increase the false alarm rate of fault diagnosis, thereby improving the accuracy of fault diagnosis.

[0009] The overall process steps are as follows:

[0010] 1. Place the inertial measurement array on a horizontal vibration isolation table and collect the angular rate signals output by the gyroscopes of the inertial measurement array. The number of IMUs is represented by m, the total number of gyroscopes is M = 3m, and the total number of sampling points is N, satisfying N ≥ 256.

[0011] 2. Extract features from the collected gyroscope angular rate signal. The specific steps are as follows:

[0012] (1) The gyroscope angular rate signal collected at the i-th sampling point is recorded as g(i):

[0013] g(i)=[g1(i),g2(i),g3(i)…g c (i)] (27)

[0014] Where g c (i) represents the time domain signal of the cth gyroscope, c = 1, 2, 3, ..., M, i = 0, 1, 2, 3, ..., N-1.

[0015] (2) For g c (i) Set K multi-modulation signals u k (i)

[0016] u k (i)=[u k,1 (i),u k,2 (i)…u k,c (i)] (28)

[0017] Where u k,c (i) represents the modal component of the kth mode of the cth gyroscope, k = 1, 2, 3, ..., K.

[0018] (3) Gyro signal g c (i) Perform discrete Fourier transform:

[0019]

[0020] Where g c (ω l ) represents the frequency domain signal of the cth gyroscope, ω l Indicates the lth actual frequency, l represents the frequency index number, f s represents the sampling frequency, l = 0, 1, 2, 3, ..., N-1, ∑ represents the summation, j represents the imaginary unit, e represents the exponent, and c = 1, 2, 3, ..., M. Let c = 1 and plot the spectrum of g1(i) after the discrete Fourier transform. K is the number of peaks in the spectrum.

[0021] (4) Initialization parameter u k,c (ω l ),ω k and λ c (ω l )

[0022] u k,c (ω l ) represents the frequency domain modal component of the cth gyroscope and the kth mode, and the initial value ω k Indicates the center frequency, initial value Take the frequency corresponding to the first peak of the spectrum after the discrete Fourier transform of g1(i); λ c (ω l ) represents the Lagrange multiplier, the initial value The alternating direction multiplier method is applied to iteratively update the modes, center frequencies and Lagrange multipliers.

[0023] (5) Modal Update

[0024]

[0025] Where, Indicates u k,c (ω l )’s value at the n+1th iteration, Indicates u k,c (ω l )’s value for the nth iteration, represents the Lagrange multiplier of the nth iteration, represents the center frequency of the nth iteration, α represents the penalty coefficient, usually 10 3 ~10 6 , k=1,2,3,…,K, c=1,2,3,…,M.

[0026] (6) Center frequency update

[0027]

[0028] Where, represents the center frequency of the (n+1)th iteration, |·| represents the absolute value, k=1,2,3,…,K, c=1,2,3,…,M.

[0029] (7) Lagrange multiplier update

[0030]

[0031] Where, represents the Lagrange multiplier of the n+1th iteration, τ represents the step size of the Lagrange multiplier, which is usually 0.5, k=1,2,3,…,K,c=1,2,3,…,M.

[0032] (8) Yes and Calculate the root mean square error:

[0033]

[0034] Where, RMSE k,c Represents root mean square error, RMSE k,c ≤10 -3 When the iterative update stops, otherwise, continue to repeat steps (5) to (8), k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., M. After the iterative update stops, the gyro signal and the modal component satisfy the following relationship:

[0035]

[0036] (9) Obtaining data after iteration

[0037] After the iterative update stops, the frequency domain modal component u k,c (ω l ) performs inverse discrete Fourier transform to obtain the time domain signals of K multivariate modulation signals:

[0038]

[0039] Where, k = 1, 2, 3,…, K, c = 1, 2, 3,…, C, i = 0, 1, 2, 3,…, N-1, l = 0, 1, 2, 3,…, N-1.

[0040] (10) Selecting characteristic modal signals

[0041] To u k,c (i) and the corresponding gyro signal g c (i) Calculate the correlation coefficient r:

[0042]

[0043] Where, Indicates root opening. and Represents u k,c (i) and g c The mean of (i), k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., C, i = 0, 1, 2, 3, ..., N-1.

[0044] Take the multivariate modulation u corresponding to the maximum value in r gbest (i) is the signal extracted from the gyroscope feature.

[0045] 3. Fault detection and isolation. The specific steps are as follows:

[0046] a. Determine the gyroscope installation configuration matrix H based on the specific installation configuration of the inertial measurement array g :

[0047]

[0048] b. Obtain the parity matrix V:

[0049] Calculate the orthogonal projection matrix W:

[0050] W=IH g (H g T H g ) -1 H g T (38)

[0051] Obtain the parity matrix V by Schmidt orthogonalization:

[0052]

[0053] Where, v q is the qth column vector of V, W q is the qth column vector of W, q∈AB, initialize the index number set A={1,2,3,…,M}, initialize the fault set

[0054] c. Calculate the residual parity vector:

[0055] p=Vu gbest (i) (40)

[0056] Where p represents the residual parity vector.

[0057] d. Calculate the fault detection function DF D , as follows:

[0058]

[0059] Where T represents the transpose and σ represents the standard deviation.

[0060] The fault decision function criterion is:

[0061]

[0062] Where, T D represents the fault detection threshold, the query freedom is M-3, and the significance level is 0.01. 2 The value corresponding to the distribution table is T D If there is a fault, go to step e for fault isolation. If there is no fault, update the sampling point index number i and repeat steps c to d of fault diagnosis and isolation.

[0063] e. Calculate the fault isolation function DFI q :

[0064]

[0065] Get DFI q The maximum value, at this time let max=q, the maximum value is DFI max , DFI max The corresponding gyroscope is the faulty gyroscope. Delete v in the parity matrix V corresponding to the faulty gyroscope. max , the remaining v q Reconstruct a parity matrix V, update the parity matrix V; delete the u corresponding to the faulty gyroscope gbest For the data points in (i), update u gbest (i) Add max to the fault set B, update the index number q, and repeat steps c to d of fault diagnosis and isolation until all collected data points are detected.

[0066] f. Output the set B of fault index numbers.

[0067] 4. Place the inertial measurement array on a horizontal vibration isolation table and collect the specific force of all accelerometers in the inertial measurement array. The number of IMUs is represented by m, the total number of accelerometers is M = 3m, and the number of sampling points is N, which should satisfy N ≥ 256. Replace the corresponding data and repeat the overall process steps 2 to 3.

[0068] The angular rate data collected at the i-th sampling point is f(i), f c (i) represents the signal of the cth accelerometer, f c (ω l ) represents the discrete frequency signal of the cth gyroscope and the lth accelerometer. The installation configuration matrix H of the accelerometer is determined by the specific installation configuration of the inertial measurement array. f , feature extraction signal u gbest (i) Replace with u fbest(i).

[0069] At this point, the design of the inertial measurement array fault diagnosis method is completed. This method can be used to diagnose faults in the inertial measurement array and ensure the reliability of the inertial measurement array.

[0070] The inventive principle of the present invention is:

[0071] The signal features of the collected data are extracted using discrete multivariate variational mode decomposition, and a generalized likelihood ratio test fault detection function is constructed to perform fault detection on all gyroscopes and accelerometers. A generalized likelihood ratio test fault isolation function is constructed to troubleshoot the gyroscopes and accelerometers one by one.

[0072] Compared with the existing solutions, the solution of the present invention has the following main advantages:

[0073] Compared with traditional fault diagnosis design methods, the present invention uses discrete multivariate variational modal decomposition to extract signal features, which can effectively solve the problem of high false alarm rate caused by noise and outliers when performing fault diagnosis in traditional algorithms, and achieve more accurate fault diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 Specific implementation plan diagram;

[0075] Figure 2 Comparison of gyroscope feature extraction signal and original signal;

[0076] Figure 3 Comparison of gyroscope fault detection functions;

[0077] Figure 4 Gyroscope raw signal fault isolation function;

[0078] Figure 5 Fault isolation function after feature extraction of gyroscope signal. DETAILED DESCRIPTION

[0079] Specific embodiments of the present invention are as follows Figure 1 As shown in the figure, an inertial measurement array consisting of 16 IMUs is taken as an example.

[0080] 1. Place the inertial measurement array on a horizontal table and collect the angular rate signals output by the inertial measurement array gyroscopes. The number of IMUs (m) is 16, the total number of gyroscopes (M) is 48, and the total number of sampling points is 1000, satisfying N ≥ 256.

[0081] 2. Extract features from the collected gyroscope angular rate signal.

[0082] 2.1. The gyroscope angular rate signal collected at the i-th sampling point is recorded as g(i):

[0083] g(i)=[g1(i),g2(i),g3(i)…g c (i)] (44)

[0084] Where g c (i) represents the time domain signal of the cth gyroscope, c = 1, 2, 3, ..., M, i = 0, 1, 2, 3, ..., N-1.

[0085] 2.2, g c (i) Set K multi-modulation signals u k (i)

[0086] u k (i)=[u k,1 (i),u k,2 (i)…u k,c (i)] (45)

[0087] Where u k,c (i) represents the modal component of the kth mode of the cth gyroscope, k = 1, 2, 3, ..., K.

[0088] 2.3、Gyro signal g c (i) Perform discrete Fourier transform:

[0089]

[0090] Where g c (ω l ) represents the frequency domain signal of the cth gyroscope, ω l Indicates the lth actual frequency, l represents the frequency index number, f s represents the sampling frequency, l = 0, 1, 2, 3, ..., N-1, ∑ represents the summation, j represents the imaginary unit, e represents the exponent, and c = 1, 2, 3, ..., M. Let c = 1 and plot the spectrum of g1(i) after the discrete Fourier transform. K is the number of peaks in the spectrum, and K = 5.

[0091] 2.4. Initialization parameter u k,c (ω l ),ω k and λ c (ω l )

[0092] u k,c (ω l ) represents the frequency domain modal component of the cth gyroscope and the kth mode, and the initial value ω k Indicates the center frequency, initial value Take the frequency corresponding to the first peak of the spectrum graph after the discrete Fourier transform of g1(i), λ c(ω l ) represents the Lagrange multiplier, the initial value The alternating direction multiplier method is applied to iteratively update the modes, center frequencies and Lagrange multipliers.

[0093] 2.5 Modal Update

[0094]

[0095] Where, Indicates u k,c (ω l )’s value at the n+1th iteration, Indicates u k,c (ω l )’s value for the nth iteration, represents the Lagrange multiplier of the nth iteration, represents the center frequency of the nth iteration, α represents the penalty coefficient, α=2000, k=1,2,3,…,K, c=1,2,3,…,M.

[0096] 2.6. Center frequency update

[0097]

[0098] Where, represents the center frequency of the (n+1)th iteration, |·| represents the absolute value, k=1,2,3,…,K, c=1,2,3,…,M.

[0099] 2.7 Lagrange multiplier update

[0100]

[0101] Where, represents the Lagrange multiplier of the n+1th iteration, τ represents the step size of the Lagrange multiplier, τ = 0.5, k = 1, 2, 3, …, K, c = 1, 2, 3, …, M.

[0102] 2.8, Yes and Calculate the root mean square error:

[0103]

[0104] Where, RMSE k,c Represents root mean square error, RMSE k,c ≤10 -3 When the iterative update stops, otherwise, continue to repeat steps 2.5 to 2.7, k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., M. After the iterative update stops, the gyro signal and the modal component satisfy the following relationship:

[0105]

[0106] 2.9. Get the data after iteration

[0107] After the iterative update stops, the frequency domain modal component u k,c (ω l ) performs inverse discrete Fourier transform to obtain the time domain signals of K multivariate modulation signals:

[0108]

[0109] Where, k = 1, 2, 3,…, K, c = 1, 2, 3,…, C, i = 0, 1, 2, 3,…, N-1, l = 0, 1, 2, 3,…, N-1.

[0110] 2.10. Selecting characteristic modal signals

[0111] To u k,c (i) and the corresponding gyro signal g c (i) Calculate the correlation coefficient r:

[0112]

[0113] Where, Indicates root opening. and Represents u k,c (i) and g c The mean of (i), k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., C, i = 0, 1, 2, 3, ..., N-1.

[0114] Take the multivariate modulation u corresponding to the maximum value in r gbest (i) is the signal extracted from the gyroscope feature.

[0115] 3. Fault detection and isolation. The specific steps are as follows:

[0116] 3.1. Determine the gyroscope installation configuration matrix H based on the specific installation configuration of the inertial measurement array g :

[0117]

[0118] 3.2. Obtain the parity matrix V:

[0119] Calculate the orthogonal projection matrix W:

[0120] W=IH g (H g T H g ) -1 Hg T (55)

[0121] Obtain the parity matrix V by Schmidt orthogonalization:

[0122]

[0123]

[0124] Where, v q is the qth column vector of V, W q is the qth column vector of W, q∈AB, the set of initialized index numbers q is A={1,2,3,…,M}, the set of initialized fault index numbers is

[0125] 3.3. Calculate the residual parity vector:

[0126] p=Vu gbest (i) (58)

[0127] Where p represents the residual parity vector.

[0128] 3.4. Calculation of fault detection function DF D , as follows:

[0129]

[0130] Where T represents the transpose and σ represents the standard deviation.

[0131] The fault decision function criterion is:

[0132]

[0133] Where, T D represents the fault detection threshold, and the query freedom is M-3. 2 The value corresponding to the distribution table is T D , T D =63.691. If there is a fault, go to step 3.5 for fault isolation. If there is no fault, update the sampling point index number i and repeat steps 3.3 to 3.4 for fault diagnosis and isolation.

[0134] 3.5. Calculation of fault isolation function DFI q :

[0135]

[0136] Get DFI q The maximum value, at this time let max=q, the maximum value is DFI max , DFI maxThe corresponding gyroscope is the faulty gyroscope. Delete v in the parity matrix V corresponding to the faulty gyroscope. max , the remaining v q Reconstruct a parity matrix V, update the parity matrix V; delete the u corresponding to the faulty gyroscope gbest For the data points in (i), update u gbest (i) Add max to the set B of fault index numbers, update the index number q, and repeat steps 3.3 to 3.4 of fault diagnosis and isolation until all collected data points are detected.

[0137] 3.6. Output the set B of fault index numbers.

[0138] 4. Place the inertial measurement array on a horizontal table and collect the specific force signal output by the accelerometer of the inertial measurement array. The number of IMUs m = 16, the total number of accelerometers M = 48, and the total number of sampling points is 1000, satisfying N ≥ 256.

[0139] 5. Perform feature extraction on the collected accelerometer force signal.

[0140] 5.1. The accelerometer force signal collected at the i-th sampling point is recorded as f(i):

[0141] f(i)=[f1(i),f2(i),f3(i)…f c (i)] (62)

[0142] Where, f c (i) represents the time domain signal of the cth accelerometer, c = 1, 2, 3, ..., M, i = 0, 1, 2, 3, ..., N-1.

[0143] 5.2, f c (i) Set K multi-modulation signals u k (i)

[0144] u k (i)=[u k,1 (i),u k,2 (i)…u k,c (i)] (63)

[0145] Where u k,c (i) represents the modal component of the kth mode of the cth accelerometer, k = 1, 2, 3, …, K.

[0146] 5.3、Accelerometer signal f c (i) Perform discrete Fourier transform:

[0147]

[0148] Where, fc (ω l ) represents the frequency domain signal of the cth accelerometer, ω l Indicates the lth actual frequency, l represents the frequency index number, f s represents the sampling frequency, l = 0, 1, 2, 3, ..., N-1, Σ represents the summation, j represents the imaginary unit, e represents the exponent, and c = 1, 2, 3, ..., M. Let c = 1 and plot the spectrum of f1(i) after the discrete Fourier transform. K is the number of peaks in the spectrum, and K = 5.

[0149] 5.4. Initialization parameter u k,c (ω l ),ω k and λ c (ω l )

[0150] u k,c (ω l ) represents the frequency domain modal component of the kth mode of the cth accelerometer, and the initial value ω k Indicates the center frequency, initial value Take the frequency corresponding to the first peak of the spectrum graph after the discrete Fourier transform of f1(i), λ c (ω l ) represents the Lagrange multiplier, the initial value The alternating direction multiplier method is applied to iteratively update the modes, center frequencies and Lagrange multipliers.

[0151] 5.5 Modal Update

[0152]

[0153] Where, Indicates u k,c (ω l )’s value at the n+1th iteration, Indicates u k,c (ω l )’s value for the nth iteration, represents the Lagrange multiplier of the nth iteration, represents the center frequency of the nth iteration, α represents the penalty coefficient, α=2000, k=1,2,3,…,K, c=1,2,3,…,M.

[0154] 5.6. Center frequency update

[0155]

[0156] Where, represents the center frequency of the (n+1)th iteration, |·| represents the absolute value, k=1,2,3,…,K, c=1,2,3,…,M.

[0157] 5.7. Lagrange multiplier update

[0158]

[0159] Where, represents the Lagrange multiplier of the n+1th iteration, τ represents the step size of the Lagrange multiplier, τ = 0.5, k = 1, 2, 3, …, K, c = 1, 2, 3, …, M.

[0160] 5.8. Yes and Calculate the root mean square error:

[0161]

[0162] Where, RMSE k,c Represents root mean square error, RMSE k,c ≤10 -3 When the iterative update stops, otherwise, continue to repeat steps 5.5 to 5.7, k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., M. After the iterative update stops, the gyro signal and the modal component satisfy the following relationship:

[0163]

[0164] 5.9. Obtaining the signal after iteration

[0165] After the iterative update stops, the frequency domain modal component u k,c (ω l ) performs inverse discrete Fourier transform to obtain the time domain signals of K multivariate modulation signals:

[0166]

[0167] Where, k = 1, 2, 3,…, K, c = 1, 2, 3,…, C, i = 0, 1, 2, 3,…, N-1, l = 0, 1, 2, 3,…, N-1.

[0168] 5.10. Selecting characteristic modal signals

[0169] To u k,c (i) and the corresponding accelerometer signal f c (i) Calculate the correlation coefficient r:

[0170]

[0171] Where, Indicates root opening. and Represents u k,c (i) and f c The mean of (i), k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., C, i = 0, 1, 2, 3, ..., N-1.

[0172] Take the multivariate modulation u corresponding to the maximum value in r fbest (i) is the signal extracted from the accelerometer features. 6. Fault detection and isolation. The specific steps are as follows:

[0173] 6.1. Determine the gyroscope installation configuration matrix H from the specific installation configuration of the inertial measurement array f :

[0174]

[0175] 6.2. Find the parity matrix V:

[0176] Calculate the orthogonal projection matrix W:

[0177] W=IH f (H f T H f ) -1 H f T (73)

[0178] Obtain the parity matrix V by Schmidt orthogonalization:

[0179]

[0180] Where, v q is the qth column vector of V, W q is the qth column vector of W, q∈AB, initializes the index number set A={1,2,3,…,M}, and initializes the set of fault index numbers

[0181] 6.3. Calculate the residual parity vector:

[0182] p=Vu fbest (i) (76)

[0183] Where p represents the residual parity vector.

[0184] 6.4. Calculation of Fault Detection Function DF D , as follows:

[0185]

[0186] Where T represents the transpose and σ represents the standard deviation.

[0187] The fault decision function criterion is:

[0188]

[0189] Where, T D represents the fault detection threshold value, and the query freedom is M-3. 2 The value corresponding to the distribution table is T D , T D =63.691. If there is a fault, go to step 6.5 for fault isolation. If there is no fault, update the sampling point index number i and repeat steps 6.3 to 6.4 for fault diagnosis and isolation.

[0190] 6.5. Calculation of Fault Isolation Function DFI q :

[0191]

[0192] Get DFI q The maximum value, at this time let max=q, the maximum value is DFI max , DFI max The corresponding accelerometer is the faulty accelerometer. Delete v in the parity matrix V corresponding to the faulty accelerometer. max , the remaining v q Reconstruct a parity matrix V, update the parity matrix V; delete the u corresponding to the faulty accelerometer fbest For the data points in (i), update u fbest (i) Add max to the set of fault index numbers B, update the index number q, and repeat steps 6.3 to 6.4 of fault diagnosis and isolation until all collected data points have been detected;

[0193] 6.6. Output the set B of fault index numbers.

[0194] At this point, the specific implementation plan for inertial measurement array fault diagnosis has been designed and this method can be used to diagnose inertial measurement array faults and ensure the reliability of the inertial measurement array.

[0195] In order to verify the effectiveness and superiority of the method of the present invention, an experiment was conducted on an inertial measurement array using a gyroscope as an example. The method of the present invention was used to extract the characteristics of the gyroscope signal, and the signal extracted from the gyroscope characteristics was compared with the original signal. Figure 2 As shown in the figure, the results show that after the feature extraction, the original gyroscope signal containing noise and outliers is converted into a stable gyroscope signal, which effectively removes the noise and outliers. A step fault with an amplitude of 3° / s is added to the gyroscope 1 at the 500th sampling point, and the gyroscope signal fault detection method of the present invention is used. The fault detection comparison is as follows: Figure 3As shown in the figure, the results show that the fault detection method of the present invention effectively eliminates the occurrence of false alarms and false alarms. The fault isolation method of the present invention is used to isolate the fault of the gyroscope. Figure 4 and Figure 5 As shown in the figure, due to the large number of gyroscopes, only four gyroscopes are shown for comparison. The results show that the fault isolation of the method of the present invention is relatively stable, effectively eliminating the situation where the traditional method is not properly isolated, and improving the accuracy of fault isolation.

[0196] Therefore, the method of the present invention can effectively solve the problem that the traditional algorithm will cause a high false alarm rate when performing fault diagnosis due to noise and outliers, and achieve more accurate fault diagnosis.

[0197] The contents not described in detail in this specification belong to the prior art known to those skilled in the art.

Claims

1. A method for diagnosing faults in an inertial measurement array, characterized in that: The method comprises: Collect the angular rate output by the three-axis gyroscope and the specific force output by the accelerometer of all inertial measurement units (IMU) in the inertial measurement array; Extract features from the collected signals; Fault detection and isolation are performed on the extracted characteristic signals.

2. The method according to claim 1, characterized in that The inertial measurement array is constructed by m IMUs in an array manner, and the three axes of each IMU are kept coincident. In a single IMU, three gyroscopes and three accelerometers are installed along three orthogonal axes respectively. The total number of gyroscopes is consistent with the total number of accelerometers, and the total number M=3m.

3. The method according to claim 1, characterized in that The acquisition of the gyroscope output is specifically as follows: the angular rate signal is collected, the total number of sampling points is N, and N ≥ 256, and the gyroscope angular rate signal collected at the i-th sampling point is recorded as g(i): g(i)=[g1(i),g2(i),g3(i)…g c (i)] (1) Where g c (i) represents the time domain signal of the cth gyroscope, c = 1, 2, 3, ..., M, i = 0, 1, 2, 3, ..., N-1.

4. The method according to claim 1, wherein The acquisition of the accelerometer output is specifically as follows: collecting the specific force signal, the total number of sampling points is N, and N ≥ 256, and the accelerometer specific force signal collected at the i-th sampling point is recorded as f(i): f(i)=[f1(i),f2(i),f3(i)…f c (i)] (2) Where, f c (i) represents the time domain signal of the cth accelerometer, c = 1, 2, 3, ..., M, i = 0, 1, 2, 3, ..., N-1.

5. The method according to claim 1, wherein The specific steps of extracting the features of the gyroscope signal are as follows: 5.1, g c (i) Set K multi-modulation signals u k (i) in k (i)=[in k,1 (and),in k,2 (and)…in k,c (and)] (3) Where u k,c (i) represents the modal component of the kth mode of the cth gyroscope, k = 1, 2, 3, ..., K; 5.

2. Gyro signal g c (i) Perform discrete Fourier transform: Where g c (ω l ) represents the frequency domain signal of the cth gyroscope, ω l Indicates the lth actual frequency, l represents the frequency index number, f s represents the sampling frequency, l = 0, 1, 2, 3, ..., N-1, Σ represents the summation, j represents the imaginary unit, e represents the exponent, c = 1, 2, 3, ..., M; let c = 1, plot the spectrum of g1(i) after discrete Fourier transform, and K is the number of peaks in the spectrum; 5.

3. Initialization parameter u k,c (ω l ),ω k and λ c (ω l ) u k,c (ω l ) represents the frequency domain modal component of the cth gyroscope and the kth mode, and the initial value ω k Indicates the center frequency, initial value Take the frequency corresponding to the first peak of the spectrum after the discrete Fourier transform of g1(i); λ c (ω l ) represents the Lagrange multiplier, the initial value Apply the alternating direction multiplier method to iteratively update the modes, center frequencies and Lagrange multipliers; 5.4 Modal Update Where, Indicates u k,c (ω l )’s value at the n+1th iteration, Indicates u k,c (ω l )’s value for the nth iteration, represents the Lagrange multiplier of the nth iteration, represents the center frequency of the nth iteration, α represents the penalty coefficient, usually 10 3 ~10 6 , k=1,2,3,…,K, c=1,2,3,…,M; 5.

5. Center frequency update Where, represents the center frequency of the n+1th iteration, |·| represents the absolute value, k=1,2,3,…,K, c=1,2,3,…,M; 5.6 Lagrange multiplier update Where, represents the Lagrange multiplier of the n+1th iteration, τ represents the Lagrange multiplier step size, usually 0.5, k=1,2,3,…,K,c=1,2,3,…,M; 5.

7. Yes and Calculate the root mean square error: Where, RMSE k,c Represents root mean square error, RMSE k,c ≤10 -3 When the iterative update stops, otherwise, continue to repeat steps 5.4 to 5.7, k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., M. After the iterative update stops, the gyro signal and the modal component satisfy the following relationship: 5.

8. Obtaining the signal after iteration After the iterative update stops, the frequency domain modal component u k,c (ω l ) performs inverse discrete Fourier transform to obtain the time domain signals of K multivariate modulation signals: Where, k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., C, i = 0, 1, 2, 3, ..., N-1, l = 0, 1, 2, 3, ..., N-1; 5.

9. Selecting characteristic modal signals To u k,c (i) and the corresponding gyro signal g c (i) Calculate the correlation coefficient r: Where, Indicates root opening. and Respectively represent u k,c (i) and g c (i) mean, k = 1, 2, 3, ..., K, c = 1, 2, 3, ..., C, i = 0, 1, 2, 3, ..., N-1; Take the multivariate modulation u corresponding to the maximum value in r gbest (i) is the signal extracted from the gyroscope feature.

6. The method according to claim 1, characterized in that The specific steps of extracting features from accelerometer signals are as follows: 6.1、f c (i) Set K multi-modulation signals u k (i), as shown in formula (3); 6.

2. Accelerometer signal f c (i) Perform discrete Fourier transform: Where, f c (ω l ) represents the frequency domain signal of the cth accelerometer; 6.

3. Initialization parameter u k,c (ω l ),ω k and λ c (ω l ), as shown in step 5.3 of claim 5; 6.4 Modal Update 6.

5. Update the center frequency, as shown in formula (6); 6.6 Lagrange multiplier update 6.

7. Yes and Calculate the root mean square error as shown in formula (8); RMSE k,c ≤10 -3 When the iterative update stops, otherwise, continue to repeat steps 6.4 to 6.

7. After the iterative update stops, the accelerometer signal and the modal component satisfy the following relationship: 6.

8. Obtaining the signal after iteration After the iterative update stops, the frequency domain modal component u k,c (ω l ) performs inverse discrete Fourier transform to obtain the time domain signals of K multivariate modulated signals, as shown in formula (10); 6.

9. Selecting characteristic modal signals To u k,c (i) and the corresponding accelerometer signal f c (i) Calculate the correlation coefficient r: Take the multivariate modulation u corresponding to the maximum value in r fbest (i) is the signal extracted from the accelerometer features.

7. The method according to claim 1, characterized in that The specific steps of performing fault detection and isolation on the gyro signal after feature extraction are as follows: 7.

1. Determine the gyroscope installation configuration matrix H from the specific installation configuration of the inertial measurement array g : 7.

2. Find the parity matrix V: Calculate the orthogonal projection matrix W: W=I-H g (H g T H g ) -1 H g T (18) Obtain the parity matrix V by Schmidt orthogonalization: Where, v q is the qth column vector of V, W q is the qth column vector of W, q∈AB, initializes the index number set A={1,2,3,…,M}, and initializes the set of fault index numbers 7.

3. Calculate the residual parity vector: p=Vu gbest (i)(20) Where p represents the residual parity vector; 7.

4. Calculation of Fault Detection Function DF D , as follows: In the formula, T means transpose, σ means standard deviation; The fault decision function criterion is: Where, T D represents the fault detection threshold, the query freedom is M-3, and the significance level is 0.

01. 2 The value corresponding to the distribution table is T D If there is a fault, go to step 7.5 for fault isolation. If there is no fault, update the sampling point index number i and repeat steps 7.3 to 7.4 for fault diagnosis and isolation. 7.

5. Calculation of Fault Isolation Function DFI q : Get DFI q The maximum value, at this time let max=q, the maximum value is DFI max , DFI max The corresponding gyroscope is the faulty gyroscope. Delete v in the parity matrix V corresponding to the faulty gyroscope. max , the remaining v q Reconstruct a parity matrix V, update the parity matrix V; delete the u corresponding to the faulty gyroscope gbest For the data points in (i), update u gbest (i) Add max to the fault set B, update the index number q, and repeat steps 7.3 to 7.4 of fault diagnosis and isolation until all collected data points are detected; 7.

6. Output the set B of fault index numbers.

8. The method according to claim 1, characterized in that The specific steps of performing fault detection and isolation on the accelerometer signal after feature extraction are as follows: 8.

1. Determine the accelerometer installation configuration matrix H from the specific installation configuration of the inertial measurement array f : 8.

2. Find the parity matrix V: Calculate the orthogonal projection matrix W: W=I-H f (H f T H f ) -1 H f T (25) Obtain the parity matrix V through Schmidt orthogonalization, as shown in formula (19); 8.

3. Calculate the residual parity vector: p=Vu fbest (i)(26) 8.

4. Calculation of Fault Detection Function DF D , as shown in formula (21); The fault decision function criterion is shown in formula (22); Query χ with M-3 degrees of freedom 2 The value corresponding to the distribution table is T D If there is a fault, go to step 8.5 for fault isolation. If there is no fault, update the sampling point index number i and repeat steps 8.3 to 8.4 for fault diagnosis and isolation. 8.

5. Calculating the Fault Isolation Function DFI q , as shown in formula (21); Get DFI q The maximum value, at this time let max=q, the maximum value is DFI max , DFI max The corresponding accelerometer is the faulty accelerometer; delete v in the parity matrix V corresponding to the faulty accelerometer. max , the remaining v q Reconstruct a parity matrix V and update the parity matrix V; Delete the faulty accelerometer corresponding to u fbest For the data points in (i), update u fbest (i) Add max to the fault set B, update the index number q, and repeat steps 8.3 to 8.4 of fault diagnosis and isolation until all collected data points are detected; 8.

6. Output the set B of fault index numbers.

Images