Aircraft positioning method and system

Abstract

This application relates to the field of aircraft positioning technology and discloses an aircraft positioning method and system. The present invention effectively solves the technical problems of positioning drift and control divergence of inertial navigation systems under high vibration environment, and realizes the accurate removal of vibration rectification error and phase lag compensation. On the one hand, by using the nonlinear model of physical constraints and the energy conservation mechanism, the DC bias generated by the sensor due to high frequency vibration is predicted and eliminated, thus eliminating false maneuvering signals. On the other hand, by deterministic group delay calibration and restricted Taylor feedforward reconstruction, the real-time performance of the signal is restored without amplifying noise. This spatiotemporal dual correction significantly improves the attitude calculation accuracy and flight control stability of the aircraft under complex working conditions.

Description

Technical Field

[0001] This application relates to the field of aircraft positioning technology, and in particular to an aircraft positioning method and system. Background Technology

[0002] With the booming development of the low-altitude economy and precision agriculture, industrial-grade drones are increasingly widely used in heavy-duty logistics transportation, large-area farmland plant protection, and long-endurance inspection. To meet the mission requirements of heavy payload and long endurance, these aircraft typically use high-power brushless motors or hybrid electric systems as propulsion devices. In particular, vertical takeoff and landing fixed-wing aircraft using internal combustion engines as power sources have more complex mechanical structures and stronger power output. However, these high-power power systems generate high-frequency and large-amplitude mechanical vibrations during operation, placing the airborne flight control system in an extremely harsh mechanical environment. As the core sensing device of the flight control system, the microelectromechanical system (MEMS) inertial measurement unit is responsible for collecting the angular velocity and acceleration information of the aircraft in real time, which is crucial for ensuring the attitude stability and navigation accuracy of the aircraft. Because the internal components of MEMS sensors contain tiny movable structures such as cantilever beams or mass blocks, these microstructures are extremely sensitive to stress changes and vibrations in the external environment. When the frequency of external mechanical vibration approaches or overlaps with the resonant frequency of the internal microstructures of the sensor, it will cause a severe resonance phenomenon, leading to a sharp decline in the quality of the sensor output signal.

[0003] In existing technical solutions, a combination of physical damping and software filtering is typically used to suppress the interference of aircraft vibration on sensors. Physical damping often employs rubber damping balls or damping plates to construct suspension devices. However, when faced with the wide-bandwidth and high-intensity vibrations of industrial-grade UAVs, rubber materials are prone to changes in damping characteristics due to aging and temperature variations, leading to damping failure or even secondary resonance. A more critical technical challenge lies in the fact that when a microelectromechanical system (MEMS) gyroscope is subjected to high-frequency sinusoidal vibration input, a phenomenon known as vibration rectification error occurs due to the nonlinearity of the motion of its internal sensing mass and the non-ideal characteristics of the electronic demodulation circuit. This error rectifies the original high-frequency vibration noise into a low-frequency or even DC drift signal, which is directly superimposed on the actual angular velocity signal, causing the flight control system to misjudge that the aircraft is continuously rotating. Although existing digital signal processing technologies generally use low-pass filters to filter out high-frequency noise, this introduces an irreconcilable contradiction: to effectively filter out high-frequency interference related to vibration rectification error, the cutoff frequency of the low-pass filter must be significantly reduced, and an extremely low cutoff frequency inevitably introduces a large phase delay into the control loop. This phase delay severely weakens the phase margin of the flight control system, making the control response sluggish when dealing with sudden gusts or performing high-maneuver flight. This can easily lead to system divergent oscillations, and in severe cases, even cause motor overheating and burnout or aircraft crash. Therefore, how to accurately isolate high-frequency vibration noise and compensate for vibration rectification errors without introducing excessive phase delay is a core technical problem that urgently needs to be solved in the field of high-dynamic aircraft navigation and control. Summary of the Invention

[0004] This application proposes an aircraft positioning method and system to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, this application adopts the following technical solution: an aircraft positioning method, comprising the following steps:

[0006] Step S1: Obtain the aircraft's body dynamic parameters, calculate the physical response cutoff frequency based on the body dynamic parameters, and use the physical response cutoff frequency to perform multi-resolution signal decomposition processing on the original inertial signals collected by the sensors to obtain the low-frequency inertial component characterizing the true motion trend of the aircraft and the high-frequency vibration component characterizing non-aircraft motion disturbances.

[0007] Step S2: Perform mean square value calculation on the high-frequency vibration components obtained in step S1 to obtain the vibration energy density value. Simultaneously acquire the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. Perform time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor.

[0008] Step S3: Input the multidimensional physical feature tensor constructed in step S2 into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained in step S2 to obtain the dynamic zero bias prediction value.

[0009] Step S4: In the time domain, use the dynamic zero-bias prediction value obtained in step S3 to perform a subtraction operation on the low-frequency inertial component obtained in step S1 to obtain a pure motion signal. Determine the algorithm's intrinsic group delay time generated by the multi-resolution signal decomposition processing in step S1. Use the algorithm's intrinsic group delay time to perform first-order derivative feedforward compensation processing on the pure motion signal to output a reconstructed positioning signal.

[0010] Furthermore, in step S1, the specific operation of calculating the physical response cutoff frequency based on the body's dynamic parameters is as follows:

[0011] A10 reads pre-stored airframe dynamics parameters from the non-volatile memory configured in the flight control computer. The airframe dynamics parameters include the maximum resultant moment modulus generated by the aircraft's propulsion system and the airframe rotational inertia tensor of the aircraft's airframe structure.

[0012] A11, perform matrix eigenvalue decomposition on the body moment of inertia tensor, and select the eigenvalue with the largest value from the result of the eigenvalue decomposition as the maximum principal moment of inertia;

[0013] A12, obtain the control loop dead zone threshold of the flight control system, take the maximum combined torque modulus as the dividend, take the product of the maximum principal moment of inertia and the control loop dead zone threshold as the divisor, and perform division to obtain the torque-inertia ratio.

[0014] A13, perform square root operation on the torque-moment of inertia ratio to obtain the angular velocity response limit value, divide the angular velocity response limit value by twice the value of pi to obtain the physical response cutoff frequency, and set the physical response cutoff frequency as the frequency domain boundary value for distinguishing the actual motion trend of the body from non-body motion interference during multi-resolution signal decomposition processing.

[0015] Furthermore, in step S1, the specific operation of performing multi-resolution signal decomposition processing on the raw inertial signal acquired by the sensor using the physical response cutoff frequency is as follows:

[0016] A20, obtain the sensor sampling rate value when the sensor collects the original inertial signal, and divide the sensor sampling rate value by twice the physical response cutoff frequency to obtain the frequency ratio value.

[0017] A21 performs a base-2 logarithmic operation on the frequency ratio value and rounds down the result to obtain the optimal decomposition depth for multi-resolution signal decomposition. A first-in-first-out sliding buffer is constructed with a correlation between the buffer depth and the optimal decomposition depth. The original inertial signal is pushed into the first-in-first-out sliding buffer in the order of acquisition time. A definite algorithm-inherent group delay time is formed during the process of pushing the original inertial signal into the first-in-first-out sliding buffer.

[0018] A22 employs an integer-domain lifting wavelet transform algorithm within a first-in-first-out sliding buffer to perform splitting, prediction, and update processing on the original inertial signal. In the prediction processing, even-numbered position samples are used to perform linear prediction operations on odd-numbered position samples and calculate the prediction error to generate high-frequency detail coefficients. In the update processing, high-frequency detail coefficients are used to perform update operations on even-numbered position samples to generate low-frequency approximation coefficients.

[0019] A23 repeatedly performs splitting, prediction, and update processing according to the optimal decomposition depth, outputs the low-frequency approximation coefficients of the last level as low-frequency inertial components, and performs inverse reconstruction and superposition processing on the high-frequency detail coefficients of all levels to obtain high-frequency vibration components.

[0020] Furthermore, in step S2, the specific operation of performing mean square value calculation on the high-frequency vibration components obtained in step S1 to obtain the vibration energy density value is as follows:

[0021] B10 reads pre-stored sensor sensitivity scaling factors from the non-volatile memory configured in the flight control computer. Sensor sensitivity scaling factors are conversion coefficients used to restore digitally quantized signals to physical acceleration values.

[0022] B11, determine the length of the energy integration window used to calculate the vibration energy, and within the time range defined by the length of the energy integration window, use the high-frequency vibration component obtained in step S1 as the multiplicand and the sensor sensitivity scaling factor as the multiplier to perform point-to-point multiplication to obtain the physical acceleration signal.

[0023] B12 performs a squaring operation on the physical acceleration signal to obtain a squared signal. Within the length of the energy integration window, the squared signal is summed to obtain the total vibration energy. The number of sampling points contained within the length of the energy integration window is obtained. The total vibration energy is used as the dividend, and the number of sampling points is used as the divisor to perform a division operation to obtain the vibration energy density value.

[0024] Furthermore, in step S2, the excitation source state parameters of the aircraft's propulsion system and the environmental physical parameters of the sensors are acquired simultaneously. The vibration energy density value, excitation source state parameters, and environmental physical parameters are then subjected to time alignment and spatial stitching processing to construct a multidimensional physical feature tensor. The specific operations are as follows:

[0025] B20 obtains motor speed data as excitation source status parameters through the communication bus of the aircraft's electronic speed controller;

[0026] B21 uses the temperature sensor built into the inertial measurement unit to obtain the core temperature data of the sensor as environmental physical parameters. For asynchronous sampling conditions where the update frequency of the excitation source state parameters is lower than the update frequency of the vibration energy density value and the update frequency of the environmental physical parameters is lower than the update frequency of the vibration energy density value, a zero-order hold strategy is adopted to perform time alignment processing on the excitation source state parameters and the environmental physical parameters.

[0027] The B22's zero-order hold strategy is to use the most recently updated excitation source state parameters and environmental physical parameters at the current moment, and read the sensor's nominal resonant frequency, resonant frequency temperature coefficient, and damping bandwidth half-width from the non-volatile memory configured in the flight control computer.

[0028] B23 uses time-aligned environmental physical parameters to perform a correction operation on the nominal resonant frequency to obtain the actual natural frequency at the current temperature;

[0029] B24. Calculate the multi-order harmonic frequencies using the time-aligned excitation source state parameters. Combine the actual natural frequency, multi-order harmonic frequencies, and damping bandwidth half-width to perform the calculation and obtain the sum of Lorentz linear resonance factors.

[0030] B25 involves splicing the vibration energy density value, the time-aligned excitation source state parameters, the sum of the Lorentz linear resonance factor, and the time-aligned environmental physical parameters in a predetermined dimensional order to construct a multidimensional physical feature tensor.

[0031] Furthermore, in step S3, the multidimensional physical feature tensor constructed in step S2 is input into the pre-configured nonlinear error prediction model, and the specific operation of outputting the dimensionless rectifier coupling coefficient through the nonlinear error prediction model is as follows:

[0032] C10 reads a pre-stored table of system physical limit parameters from the non-volatile memory configured in the flight control computer;

[0033] C11, the system physical limit parameter table includes the maximum theoretical energy density of the sensor, the maximum physical speed of the power system and the maximum range of ambient temperature variation. The physical boundary scaling diagonal matrix and the physical reference center vector are constructed using the maximum theoretical energy density of the sensor, the maximum physical speed of the power system and the maximum range of ambient temperature variation.

[0034] C12, the physical reference center vector performs centering processing on the environmental physical parameters, and the physical reference center vector keeps the absolute zero value of the sum of the vibration energy density value, the excitation source state parameters and the Lorentz linear resonance factor unchanged;

[0035] C13, the diagonal elements of the physical boundary scaling diagonal matrix are the reciprocal of the sensor's maximum theoretical energy density, the reciprocal of the power system's maximum physical speed, the unit value, and the reciprocal of the maximum range of ambient temperature variation;

[0036] C14 performs a physical limit-based linear manifold mapping operation on the multidimensional physical feature tensor using a physical boundary scaling diagonal matrix and a physical reference center vector to generate a dimensionless normalized input tensor. This dimensionless normalized input tensor is then fed into a nonlinear error prediction model, which employs a gated recurrent unit network structure with memory cells.

[0037] C15 uses the hyperbolic tangent activation function configured in the nonlinear error prediction model to perform amplitude limiting on the calculation results of the nonlinear error prediction model, and obtains the output value of the hyperbolic tangent activation function.

[0038] C16 reads the pre-stored maximum rectified gain limit from the non-volatile memory configured in the flight control computer, performs a multiplication operation on the output value of the hyperbolic tangent activation function and the maximum rectified gain limit, and obtains the rectified coupling coefficient.

[0039] Furthermore, the corresponding axial vibration rectification error sensitivity reference is read from the non-volatile memory configured in the flight control computer;

[0040] C20 corresponds to the axial vibration rectification error sensitivity benchmark, which is the unit energy error conversion rate of the target axial direction at the most sensitive frequency point, measured by frequency sweep bench test during the factory calibration stage.

[0041] C21, corresponding to the axial vibration rectification error sensitivity benchmark, has the physical dimension balancing function of restoring the square energy term to the first acceleration term;

[0042] C22, use the dimensionless rectifier coupling coefficient output in step S3 as the first multiplier;

[0043] C23, using the corresponding axial vibration rectification error sensitivity benchmark as the second multiplier;

[0044] C24, use the vibration energy density value obtained in step S2 as the third multiplier;

[0045] C25 performs a three-term multiplication operation on the first, second, and third multipliers to obtain the dynamic zero-bias prediction value. The three-term multiplication operation forces the dynamic zero-bias prediction value to be zero when the vibration energy density is zero.

[0046] Furthermore, in step S4, the specific operation of subtracting the low-frequency inertial component obtained in step S1 from the dynamic zero-bias prediction value obtained in step S3 in the time domain to obtain the pure motion signal is as follows:

[0047] D10, confirm that the low-frequency inertial component obtained in step S1 and the dynamic zero-bias prediction value obtained in step S3 are both constrained by the first-in-first-out sliding buffer mechanism in step S1. Under the premise that the low-frequency inertial component and the dynamic zero-bias prediction value have the same phase lag property, take the low-frequency inertial component as the minuend and the dynamic zero-bias prediction value as the subtrahend, and perform a subtraction operation on the minuend and the subtrahend to obtain a pure motion signal.

[0048] D11, the specific operation to determine the inherent group delay time of the algorithm generated by the multi-resolution signal decomposition processing in step S1 is as follows: obtain the optimal decomposition depth and sensor sampling rate value determined in step S1, perform a power operation on the optimal decomposition depth to obtain the causal buffer cost value introduced by the first-in-first-out sliding buffer, and obtain the filter tap length value of the lifting wavelet filter.

[0049] D12, subtract one from the filter tap length value to get the filter length difference value, divide the filter length difference by two to get the filter phase delay value;

[0050] D13 adds the causal buffer cost value to the filter phase delay value to obtain the total delay sampling points. Dividing the total delay sampling points by the sensor sampling rate value yields the algorithm's intrinsic group delay time.

[0051] Furthermore, in step S4, the specific operation of performing first-order derivative feedforward compensation processing on the pure motion signal using the algorithm's inherent group delay time to output the reconstructed positioning signal is as follows:

[0052] D20, construct a restricted Taylor reconstruction operator. The restricted Taylor reconstruction operator uses the three-point Lagrange backward difference algorithm to perform differential operations on the pure motion signal to obtain the rate of change value of the pure motion signal.

[0053] D21 reads the pre-stored maximum physical acceleration value of the aircraft from the non-volatile memory configured in the flight control computer, and performs kinematic truncation processing on the rate of change value using the maximum physical acceleration value of the aircraft.

[0054] D22, the specific steps of the kinematic truncation process are as follows: determine whether the absolute value of the rate of change is greater than the maximum physical acceleration value of the aircraft. If the absolute value of the rate of change is greater than the maximum physical acceleration value of the aircraft, then set the rate of change to the maximum physical acceleration value of the aircraft and output it as the safe acceleration value. If the absolute value of the rate of change is not greater than the maximum physical acceleration value of the aircraft, then keep the rate of change unchanged and output it as the safe acceleration value.

[0055] D23, using the safety acceleration value as the multiplicand and the algorithm's intrinsic group delay time as the multiplier, performs a multiplication operation on the safety acceleration value and the algorithm's intrinsic group delay time to obtain the time-domain feedforward compensation amount;

[0056] D24 performs an addition operation on the pure motion signal and the time-domain feedforward compensation, and outputs the reconstructed positioning signal.

[0057] An aircraft positioning system includes: a signal decomposition module, a feature construction module, an error prediction module, and a signal reconstruction module, wherein;

[0058] The signal decomposition module is used to acquire the aircraft's body dynamic parameters, calculate the physical response cutoff frequency based on the body dynamic parameters, and perform multi-resolution signal decomposition processing on the raw inertial signals collected by the sensors using the physical response cutoff frequency to obtain the low-frequency inertial component that characterizes the true motion trend of the aircraft and the high-frequency vibration component that characterizes non-aircraft motion interference.

[0059] The feature construction module is used to perform mean square value calculation on the high-frequency vibration components obtained by the signal decomposition module to obtain the vibration energy density value. Simultaneously, it acquires the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. It performs time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor.

[0060] The error prediction module is used to input the multidimensional physical feature tensor constructed by the feature construction module into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained by the feature construction module to obtain the dynamic zero bias prediction value.

[0061] The signal reconstruction module is used to perform a subtraction operation on the low-frequency inertial component obtained by the signal decomposition module in the time domain using the dynamic zero-bias prediction value obtained by the error prediction module to obtain a clean motion signal. It determines the inherent group delay time of the algorithm generated by the multi-resolution signal decomposition processing in the signal decomposition module, performs first-order derivative feedforward compensation processing on the clean motion signal using the inherent group delay time of the algorithm, and outputs the reconstructed positioning signal.

[0062] The beneficial effects of this invention are as follows:

[0063] This invention effectively solves the technical problems of positioning drift and control divergence in inertial navigation systems under high vibration environments. It achieves precise removal of vibration rectification errors and compensation for phase lag. On the one hand, by using a nonlinear model of physical constraints and energy conservation mechanism, it predicts and eliminates the DC bias generated by the sensor due to high-frequency vibration, thus eliminating false maneuvering signals. On the other hand, through deterministic group delay calibration and restricted Taylor feedforward reconstruction, it restores the real-time performance of the signal without amplifying noise. This spatiotemporal dual correction significantly improves the attitude calculation accuracy and flight control stability of the aircraft under complex operating conditions. Attached Figure Description

[0064] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort:

[0065] Figure 1 This is a flowchart of the method of the present invention;

[0066] Figure 2 This is a system framework diagram of the present invention. Detailed Implementation

[0067] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0068] Example 1

[0069] like Figure 1 As shown, this invention discloses a method for locating an aircraft, comprising the following steps:

[0070] Step S1: Obtain the aircraft's dynamic parameters, calculate the physical response cutoff frequency based on the dynamic parameters, and use the physical response cutoff frequency to perform multi-resolution signal decomposition processing on the raw inertial signals collected by the sensors to obtain low frequencies that characterize the true motion trend of the aircraft.

[0071] In this embodiment, the specific operation of calculating the physical response cutoff frequency based on the body dynamics parameters in step S1 is as follows: read the pre-stored body dynamics parameters from the non-volatile memory configured in the flight control computer. The body dynamics parameters include the maximum resultant torque modulus generated by the aircraft power system and the body rotational inertia tensor of the aircraft body structure.

[0072] In the specific execution process, the method first accesses the electrically erasable programmable read-only memory (EEPROM) or flash memory of the flight control computer to read the airframe dynamics configuration table written during the factory calibration stage.

[0073] In this embodiment, the maximum resultant moment modulus is selected as a fixed value between 30 and 50 Newton-meters. The selection of the maximum resultant moment modulus is based on the bench limit test data of the aircraft's power system, which represents the physical upper limit of energy input that all actuators of the aircraft can provide under full throttle conditions. The airframe moment of inertia tensor is selected as a 3x3 second-order tensor matrix obtained by simulation calculation through computer-aided design software. The airframe moment of inertia tensor accurately characterizes the inertial properties of the aircraft's mass distribution in three-dimensional space.

[0074] Perform matrix eigenvalue decomposition on the moment of inertia tensor of the machine body, and select the eigenvalue with the largest value from the results of the eigenvalue decomposition as the maximum principal moment of inertia.

[0075] During the specific execution process, the method calls the linear algebra library to perform Jacobi iteration or QR decomposition on the three-by-three body rotational inertia tensor to calculate three positive real eigenvalues. The method then executes numerical comparison logic to lock the maximum value among the three eigenvalues ​​as the maximum principal rotational inertia.

[0076] In this embodiment, the technical basis for selecting the maximum principal moment of inertia is to follow the principle of the most unfavorable operating condition, that is, to take the rotation axis (usually the yaw axis) of the aircraft with the slowest response as the bottleneck of physical response capability, to ensure that the cutoff frequency calculated subsequently covers the physical hysteresis characteristics of all axes, and to prevent the filter bandwidth from being set too high due to underestimation of inertia, thereby introducing high-frequency noise from non-body motion.

[0077] The dead zone threshold of the control loop of the flight control system is obtained. The maximum resultant torque modulus is used as the dividend, and the product of the maximum principal moment of inertia and the dead zone threshold of the control loop is used as the divisor. The division operation is performed to obtain the torque-inertia ratio.

[0078] In the specific implementation process, the control loop dead zone threshold is preferably set to 0.017 radians (approximately equal to one degree) in this embodiment. The selection of the control loop dead zone threshold is based on the response characteristics of the flight control algorithm (such as a proportional-integral-derivative controller): for small attitude jitters with an amplitude less than 0.017 radians, the flight control algorithm is physically set to not respond to avoid servo jitter. If the control loop dead zone threshold is not introduced, when the amplitude approaches zero, the theoretical physical response frequency will tend to infinity, causing the mathematical model to collapse. Therefore, the introduction of the control loop dead zone threshold couples and anchors the physical limit with the control requirements. The method performs a division operation to calculate the quotient of the maximum combined torque modulus divided by the product of the maximum principal moment of inertia and the control loop dead zone threshold, and obtains the torque-inertia ratio with the square dimension of the angular acceleration frequency.

[0079] The angular velocity response limit is obtained by taking the square root of the torque-moment of inertia ratio. The angular velocity response limit is then divided by twice the value of pi to obtain the physical response cutoff frequency. This physical response cutoff frequency is set as the frequency domain boundary value used to distinguish between the actual motion trend of the machine body and non-machine body motion interference during multi-resolution signal decomposition processing.

[0080] In practice, the method calculates the arithmetic square root of the torque-inertia ratio and divides the result by 6.28 (i.e., twice pi), thereby converting the angular frequency into the physical response cutoff frequency in Hertz.

[0081] In the heavy-duty UAV scenario of this embodiment, the calculated physical response cutoff frequency is usually distributed in the closed interval of 5 Hz to 40 Hz. Compared with the prior art scheme of setting a fixed 20 Hz low-pass filter, the physical response cutoff frequency calculated in this embodiment is dynamic and has physical reference: any signal energy with a frequency higher than the physical response cutoff frequency, according to Newton's second law, cannot be physically generated by a motor with limited power driving a large inertia body, and therefore must be judged as non-body motion interference.

[0082] This first-principles-based parameter definition method effectively solves the technical problems of phase delay or vibration aliasing caused by fixed parameters in traditional filters.

[0083] In this embodiment, the specific operation of performing multi-resolution signal decomposition processing on the original inertial signal collected by the sensor using the physical response cutoff frequency in step S1 is as follows: obtain the sensor sampling rate value when the sensor collects the original inertial signal, and divide the sensor sampling rate value by twice the physical response cutoff frequency to obtain the frequency ratio value.

[0084] In the specific execution process, the method reads the sensor sampling rate value from the analog-to-digital converter configuration register. In this embodiment, the sensor sampling rate value is preferably set to 1 kHz to meet the Nyquist sampling theorem's requirement for capturing high-frequency vibration details. The method calculates the ratio of the sensor sampling rate value to twice the physical response cutoff frequency. This frequency ratio value quantifies the multiple relationship between the physical limit frequency and the digital sampling frequency.

[0085] A base-2 logarithmic operation is performed on the frequency ratio value, and the result of the logarithmic operation is rounded down to obtain the optimal decomposition depth for multi-resolution signal decomposition processing. A first-in-first-out sliding buffer with a correlation between the buffer depth and the optimal decomposition depth is constructed. The original inertial signal is pushed into the first-in-first-out sliding buffer in the order of acquisition time. A definite algorithm-inherent group delay time is formed during the process of pushing the original inertial signal into the first-in-first-out sliding buffer.

[0086] In practice, the calculation of the optimal decomposition depth ensures that the upper limit of the frequency band of the Nth decomposition is close to the physical response cutoff frequency. In order to solve the problem of non-causality that the prediction step in the wavelet transform needs to use future data, the method allocates a first-in-first-out sliding buffer in the random access memory.

[0087] In this embodiment, the buffer depth of the first-in-first-out sliding buffer is set to the optimal decomposition depth plus a power of one sample point. The residence process of the original inertial signal in the first-in-first-out sliding buffer physically introduces a definite algorithmic intrinsic group delay time. This algorithmic intrinsic group delay time is pre-calculated and constant, aiming to exchange a definite time lag for the accuracy of frequency domain processing, thereby transforming the mathematically non-causal filtering algorithm into an engineering-implementable causal processing logic.

[0088] Within a first-in-first-out sliding buffer, an integer-domain lifting wavelet transform algorithm is used to perform splitting, prediction, and update processing on the original inertial signal. In the prediction processing, even-numbered position samples are used to perform linear prediction operations on odd-numbered position samples and calculate the prediction error to generate high-frequency detail coefficients. In the update processing, high-frequency detail coefficients are used to perform update operations on even-numbered position samples to generate low-frequency approximation coefficients.

[0089] In specific implementation, this embodiment preferably adopts The method uses a wavelet basis and restricts all operations to the integer domain. It replaces division with bit shifting operations (e.g., right shift by one bit instead of division by two), thus avoiding the accumulation of truncation errors caused by floating-point operations.

[0090] In the predictive processing, the method calculates the average of the current even-numbered position sample and the next even-numbered position sample in the buffer. The odd-numbered position samples are then subtracted from this average to obtain the high-frequency detail coefficients, which capture transient components in the signal whose rate of change exceeds the physical inertial constraints.

[0091] In the update process, the method uses high-frequency detail coefficients to smooth and correct even-numbered position samples to obtain low-frequency approximation coefficients, so as to maintain the energy conservation characteristics of the signal.

[0092] Repeatedly perform splitting, prediction, and update processing according to the optimal decomposition depth, output the low-frequency approximation coefficients of the last level as low-frequency inertial components, and perform inverse reconstruction and superposition processing on the high-frequency detail coefficients of all levels to obtain high-frequency vibration components.

[0093] In the specific execution process, the method recursively calls the above logic until the optimal decomposition depth is reached. At this point, the low-frequency approximation coefficients output are the low-frequency inertial components. This component physically represents the pure body motion signal that lags behind the inherent group delay time of the algorithm. At the same time, the method reconstructs all the high-frequency detail coefficients generated at each level back to the time domain through the inverse lifting operator and performs linear superposition to generate high-frequency vibration components.

[0094] In summary, this embodiment solves the technical problem of the overlap and difficulty in separating body motion and environmental vibration in the frequency domain by introducing cutoff frequency calculation logic based on the physical properties of the body and causal wavelet transform logic based on the buffer. By using the physical limit as the hard boundary of signal decomposition, this embodiment can completely separate non-body motion interference (i.e., high-frequency vibration components) as independent and orthogonal data sources without losing real maneuver information. This provides a clean, reliable and physically meaningful data foundation for constructing a high-precision physical feature tensor in the subsequent step S2 and for error prediction based on energy conservation in step S3.

[0095] Step S2: Perform mean square value calculation on the high-frequency vibration components obtained in step S1 to obtain the vibration energy density value. Simultaneously acquire the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. Perform time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor.

[0096] In step S2, the mean square value calculation is performed on the high-frequency vibration component obtained in step S1 to obtain the vibration energy density value. The specific operation is as follows: read the pre-stored sensor sensitivity scaling factor from the non-volatile memory configured in the flight control computer. The sensor sensitivity scaling factor is a conversion coefficient used to restore the digital quantization format signal to the physical acceleration value.

[0097] In this embodiment, for a microelectromechanical system accelerometer with a range configuration of ±16 times the gravitational acceleration (±16g) and a 16-bit analog-to-digital converter, the sensor sensitivity scaling factor is set to 0.000488 (i.e., 16 divided by 32,768). The technical basis for selecting the sensor sensitivity scaling factor is to establish a linear mapping relationship between digital sampling values ​​and the real physical environment. In some existing technologies, some solutions directly use digital quantization values ​​for feature extraction, ignoring the impact of changes in sensor range configuration on energy calculation. By introducing the sensor sensitivity scaling factor, this method ensures that subsequent calculations are based on the real physical dimension of meters per second squared, making the vibration energy calculation results under different range configurations physically consistent.

[0098] Determine the length of the energy integration window used to calculate the vibration energy. Within the time range defined by the energy integration window length, use the high-frequency vibration component obtained in step S1 as the multiplicand and the sensor sensitivity scaling factor as the multiplier to perform point-to-point multiplication to obtain the physical acceleration signal.

[0099] In this embodiment, the energy integration window length is set to fifty sampling points (corresponding to a sampling time of fifty milliseconds). The selection of the energy integration window length parameter follows the matching principle of the Nyquist sampling theorem and the low-frequency cutoff frequency: the time length of fifty milliseconds completely covers the time scale of the lowest frequency band (about forty milliseconds) in the multi-resolution signal decomposition processing in step S1 in the time domain, ensuring that the energy integration window length contains the complete vibration cycle.

[0100] Comparative experimental data show that when the energy integration window length is less than 20 milliseconds, the variance of the calculated vibration energy density value increases by 30%, exhibiting significant random jitter; while when the energy integration window length is set to 50 milliseconds, the smoothness and response speed of the vibration energy density value reach the best balance, which can accurately reflect the changing trend of instantaneous vibration power.

[0101] The physical acceleration signal is squared to obtain a squared signal. The squared signal is then summed within the energy integration window to obtain the total vibration energy. The number of sampling points included in the energy integration window is obtained. The total vibration energy is used as the dividend, and the number of sampling points is used as the divisor to perform a division operation to obtain the vibration energy density value.

[0102] In this embodiment, the method strictly performs square operation rather than root mean square operation. The parameter selection is based on the second-order nonlinear rectification principle in the failure physics of microelectromechanical systems: the magnitude of the vibration rectification error is linearly proportional to the square of the input vibration acceleration (i.e., the energy term), rather than the amplitude.

[0103] Existing technologies typically use root mean square (RMS) values ​​as feature inputs, forcing subsequent models to refit square relationships, increasing computational complexity and reducing fitting accuracy. Through comparative testing, under the same vibration conditions, using vibration energy density values ​​calculated from square signals as input features resulted in a 15.4% reduction in prediction residuals for subsequent nonlinear error prediction models compared to using RMS values. This demonstrates the beneficial effect of feature construction based on first principles. The vibration energy density value accurately quantifies the average vibration kinetic energy borne by a unit mass of sensor microstructure at the current moment, providing physical energy data that directly characterizes the causes of errors in subsequent steps.

[0104] In step S2, the excitation source state parameters of the aircraft's power system and the environmental physical parameters of the sensors are acquired synchronously. The vibration energy density value, excitation source state parameters, and environmental physical parameters are then subjected to time alignment and spatial stitching to construct a multidimensional physical feature tensor. The specific operations are as follows: the motor speed data is acquired through the communication bus of the aircraft's electronic speed controller as the excitation source state parameters; the core temperature data of the sensor is acquired through the temperature sensor built into the inertial measurement unit as the environmental physical parameters. For asynchronous sampling conditions where the update frequency of the excitation source state parameters is lower than that of the vibration energy density value and the update frequency of the environmental physical parameters is lower than that of the vibration energy density value, a zero-order hold strategy is used to perform time alignment processing on the excitation source state parameters and the environmental physical parameters.

[0105] In this embodiment, the update frequency of the vibration energy density value is 1 kHz, while the update frequency of the motor speed data is 50 Hz, and the update frequency of the sensor core temperature data is 20 Hz. Existing technologies often use linear interpolation algorithms for alignment, but linear interpolation algorithms require data points from future moments, which violates the causality of real-time control systems. The zero-order hold strategy adopted in this method only needs to read the latest written value in the register in physical implementation, without waiting for the next frame of data, eliminating the interpolation waiting delay of about 20 to 50 milliseconds, and meeting the stringent requirements of flight control systems for microsecond-level real-time performance.

[0106] The zero-order hold strategy uses the most recently updated excitation source state parameters and environmental physical parameters at the current moment, and reads the sensor's nominal resonant frequency, resonant frequency temperature coefficient, and damping bandwidth half-width from the non-volatile memory configured in the flight control computer.

[0107] In this embodiment, the nominal resonant frequency is set to 27 kHz, the temperature coefficient of the resonant frequency is set to -30 parts per million per degree Celsius, and the damping bandwidth half-width is set to 50 Hz. The selection of the above parameters is based on the wafer-level test report of the microelectromechanical system sensor.

[0108] The nominal resonant frequency is corrected using time-aligned environmental physical parameters to obtain the actual natural frequency at the current temperature. The multi-harmonic frequencies are calculated using time-aligned excitation source state parameters. The actual natural frequency, multi-harmonic frequencies, and damping bandwidth half-width are combined to obtain the sum of Lorentz linear resonance factors.

[0109] In this embodiment, the method first uses the current temperature to correct the nominal resonant frequency. The calculation formula is the nominal resonant frequency multiplied by (one plus the temperature coefficient of the resonant frequency multiplied by the temperature change) to obtain the actual natural frequency.

[0110] Theoretical analysis shows that the Young's modulus of the micromechanical structure decreases with increasing temperature, causing the natural frequency to drift. Without temperature correction, the drift in natural frequency can reach 40 Hz when the temperature difference reaches 50 degrees Celsius, which is enough to cause the resonance peak prediction to fail.

[0111] Subsequently, the method calculates the fundamental frequency, second harmonic frequency, and third harmonic frequency of the motor speed as the multiple harmonic frequencies, and uses the Lorentz function to calculate the overlap between each harmonic and the actual natural frequency, and accumulates them to obtain the cumulative sum of the Lorentz linear resonance factors.

[0112] Comparative experiments show that, during the process of the motor accelerating through the resonance point, the scheme that introduces the sum of Lorentz linear resonance factors improves the ability to capture resonance mutation errors by 40% compared to the scheme that only uses speed characteristics, effectively solving the technical problem that a single speed characteristic cannot characterize nonlinear resonance peaks.

[0113] The vibration energy density value, the time-aligned excitation source state parameters, the sum of the Lorentz linear resonance factor, and the time-aligned environmental physical parameters are spliced ​​together in a predetermined dimensional order to construct a multidimensional physical feature tensor.

[0114] In this embodiment, the constructed multidimensional physical feature tensor is a four-dimensional column vector. This multidimensional physical feature tensor achieves strict synchronization of heterogeneous data in the time dimension and integrates the four core physical elements that cause vibration rectification error in the spatial dimension: energy source (vibration energy density value), excitation source (excitation source state parameters), coupling source (Lorentz linear resonance factor accumulation) and environmental source (environmental physical parameters).

[0115] The method in this embodiment transforms the original sensor data, which is discrete, asynchronous, and has a single physical meaning, into a set of features that are spatiotemporally aligned, have clear physical causality, and contain high-order nonlinear prior knowledge. This provides complete and pure input data for the nonlinear error prediction model in the subsequent step S3, significantly reducing the training difficulty of the model and improving its generalization ability.

[0116] Step S3: Input the multidimensional physical feature tensor constructed in step S2 into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained in step S2 to obtain the dynamic zero bias prediction value.

[0117] In step S3, the multidimensional physical feature tensor constructed in step S2 is input into the pre-configured nonlinear error prediction model, and the specific operation of outputting the dimensionless rectification coupling coefficient through the nonlinear error prediction model is as follows: read the pre-stored system physical limit parameter table from the non-volatile memory configured in the flight control computer.

[0118] In this embodiment, the system physical limit parameter table stored in the non-volatile memory is based on the absolute boundary values ​​determined during the design phase of the aircraft hardware physical properties. The technical purpose of introducing the system physical limit parameter table is to solve the problem that neural networks have difficulty processing raw physical data with huge dimensional spans. For example, the vibration energy density is on the order of 10 to the power of negative cube, while the rotational speed of the propulsion system is on the order of 10 to the power of cube. If physical normalization is not performed, it will lead to an imbalance in weight updates during the gradient descent process.

[0119] The system physical limit parameter table includes the sensor's maximum theoretical energy density, the power system's maximum physical rotation speed, and the maximum range of ambient temperature variation. The physical boundary scaling diagonal matrix and the physical reference center vector are constructed using the sensor's maximum theoretical energy density, the power system's maximum physical rotation speed, and the maximum range of ambient temperature variation.

[0120] In this embodiment, the maximum theoretical energy density of the sensor is set to 24,601.6 meters per second squared. This value is calculated based on the square of the inertial sensor's range (e.g., ±16 times the gravitational acceleration) and represents the upper limit of energy that the hardware can sense. Inputs exceeding this value are considered sensor saturation. The maximum physical rotational speed of the power system is set to 100 Hz, calculated based on the motor's back electromotive force constant and the battery voltage. The maximum range of ambient temperature variation is set to 80 degrees Celsius (i.e., -40 degrees Celsius to +40 degrees Celsius). These parameters define the hard boundaries of the physical world.

[0121] The physical reference center vector performs centering processing on the environmental physical parameters, and the physical reference center vector keeps the absolute zero value of the sum of the vibration energy density value, the excitation source state parameters and the Lorentz linear resonance factor unchanged.

[0122] In this embodiment, the fourth dimension component (corresponding to temperature) of the physical reference center vector is set to the calibration temperature (e.g., 25 degrees Celsius), while the first three dimension components (corresponding to energy, frequency, and resonance factor) are set to zero.

[0123] The basis for this asymmetric centering process is that the effect of temperature on micromechanical structures manifests as a bidirectional drift relative to the calibration point, thus requiring mean removal. Energy and frequency, on the other hand, have an absolute zero-point physical meaning (i.e., a static state), and translation operations are strictly prohibited, otherwise the sparsity of the data will be destroyed, causing the nonlinear error prediction model to receive non-zero input noise in a static state, thereby generating false drift predictions.

[0124] The diagonal elements of the physical boundary scaling diagonal matrix are the reciprocal of the sensor's maximum theoretical energy density, the reciprocal of the power system's maximum physical rotational speed, the unit value, and the reciprocal of the maximum range of ambient temperature variation.

[0125] In this embodiment, the processor constructs a 4x4 diagonal matrix and maps each physical component to a dimensionless interval from negative one to positive one through matrix multiplication.

[0126] A linear manifold mapping operation based on physical limits is performed on the multidimensional physical feature tensor using the physical boundary scaling diagonal matrix and the physical reference center vector to generate a dimensionless normalized input tensor. The dimensionless normalized input tensor is then input into a nonlinear error prediction model, which employs a gated recurrent unit network structure with memory units.

[0127] In this embodiment, the processor strictly executes linear mapping operations and refuses to use logarithmic transformations. Comparative experimental data shows that since the vibration rectification error and the vibration energy density value follow the principle of linear superposition, if nonlinear transformations (such as logarithmic transformations) are used to process the energy term, the gradient in the high-energy range will be compressed, resulting in a 25% increase in the prediction residual of the model under severe vibration conditions. The physical basis for selecting the gated recurrent unit network structure is that the change of Young's modulus of the micromechanical structure has thermal hysteresis characteristics, that is, the stress release paths of the heating process and the cooling process are different. The gated recurrent unit can capture this thermodynamic history dependency relationship across time steps through its internal state memory unit, thereby significantly improving the prediction accuracy under variable temperature conditions.

[0128] The hyperbolic tangent activation function configured by the nonlinear error prediction model is used to perform amplitude limiting on the calculation results of the nonlinear error prediction model, and the output value of the hyperbolic tangent activation function is obtained.

[0129] In this embodiment, considering that the response gain of the physical system is finite, the processor forces the use of the hyperbolic tangent activation function in the output layer of the nonlinear error prediction model, strictly limiting the output range to the interval between negative one and positive one, thus eliminating the risk of numerical explosion from the algorithm structure.

[0130] The maximum rectified gain limit is read from the non-volatile memory configured in the flight control computer. The output value of the hyperbolic tangent activation function is multiplied by the maximum rectified gain limit to obtain the rectified coupling coefficient.

[0131] In this embodiment, the maximum rectified gain limit is a dimensionless constant, preferably set to 1.2. The value is selected based on the ratio of the error sensitivity to the nominal sensitivity under the worst operating condition in the device datasheet. The rectified coupling coefficient characterizes the conversion efficiency and polarity of the microstructure in converting unit vibration energy into DC bias under the current frequency and temperature combination.

[0132] Introducing a maximum rectified gain limit as a physical constraint prevents the nonlinear error prediction model from outputting coefficients that exceed physical possibilities when dealing with extreme data outside the training set.

[0133] In step S3, the specific operation of multiplying the rectification coupling coefficient with the vibration energy density value obtained in step S2 to obtain the dynamic zero bias prediction value is as follows: read the pre-stored corresponding axial vibration rectification error sensitivity reference from the non-volatile memory configured in the flight control computer.

[0134] In this embodiment, the processor reads the corresponding axial vibration rectification error sensitivity reference based on the currently processed sensitive axis (e.g., the Z-axis).

[0135] The corresponding axial vibration rectification error sensitivity benchmark is the unit energy error conversion rate of the target axial direction at the most sensitive frequency point, which is measured by frequency sweep bench testing during the factory calibration stage.

[0136] In this embodiment, the sensitivity benchmark for the axial vibration rectification error is obtained by applying a full-band scan to the sensor on a vibration table, capturing the ratio of the drift amount that causes the maximum zero-bias drift to the vibration energy. This value is stored in a non-volatile memory and represents the potential pathological degree of the sensor under the worst-case scenario.

[0137] The corresponding axial vibration rectification error sensitivity benchmark has the physical dimension balancing function of restoring the square energy term to the first-order acceleration term.

[0138] In this embodiment, the dimension of the vibration energy density value is meters per second squared, while the dimension of the target dynamic zero bias prediction value is meters per second squared. Therefore, the physical dimension of the corresponding axial vibration rectification error sensitivity benchmark is seconds squared per meter. Introducing the corresponding axial vibration rectification error sensitivity benchmark solves the dimensional paradox that the dimensionless coefficient multiplied by energy is not equal to the error, ensuring the dimensional homogeneity of both sides of the physical equation.

[0139] The dimensionless rectification coupling coefficient output in step S3 is used as the first multiplier; the corresponding axial vibration rectification error sensitivity benchmark is used as the second multiplier; the vibration energy density value obtained in step S2 is used as the third multiplier; the three-term multiplication operation is performed on the first multiplier, the second multiplier and the third multiplier to obtain the dynamic zero bias prediction value. The three-term multiplication operation forces the dynamic zero bias prediction value to be zero when the vibration energy density value is zero.

[0140] In this embodiment, the processor performs a three-term multiplication operation to construct an embedded zero-input zero-output physical constraint equation. In this equation, the rectified coupling coefficient output by the nonlinear error prediction model essentially acts as an attenuator.

[0141] For example, when at a non-resonant frequency, the network outputs a small coefficient (e.g., 0.1), indicating that the current error conversion efficiency is only 10% of the worst-case scenario. Crucially, when there is no external vibration energy input (i.e., the vibration energy density as the third multiplier is zero), regardless of how the internal state of the nonlinear error prediction model fluctuates or what rectifier coupling coefficients are output, the final dynamic zero-bias prediction value is strictly zero. This mechanism fundamentally eliminates the risk of static drift illusion that a purely data-driven model might produce when at rest.

[0142] Comparative experiments show that after introducing this energy conservation hard constraint, the miscompensation rate of the positioning method in the static phase is reduced from five percent to zero, which significantly improves the robustness of the system.

[0143] Through the feature manifold mapping and physical constraint synthesis in this embodiment, the present invention successfully combines the nonlinear fitting capability of deep learning with the rigid constraints of physical laws, outputting a high-confidence, physically interpretable dynamic zero-bias prediction value, providing a reliable numerical basis for the accurate subtraction correction of the original inertial signal in step S4.

[0144] Step S4: In the time domain, use the dynamic zero-bias prediction value obtained in step S3 to perform a subtraction operation on the low-frequency inertial component obtained in step S1 to obtain a pure motion signal. Determine the algorithm's intrinsic group delay time generated by the multi-resolution signal decomposition processing in step S1. Use the algorithm's intrinsic group delay time to perform first-order derivative feedforward compensation processing on the pure motion signal to output a reconstructed positioning signal.

[0145] In step S4, the specific operation of subtracting the low-frequency inertial component obtained in step S1 from the dynamic zero-bias prediction value obtained in step S3 in the time domain to obtain a pure motion signal is as follows: confirm that both the low-frequency inertial component obtained in step S1 and the dynamic zero-bias prediction value obtained in step S3 are constrained by the first-in-first-out sliding buffer mechanism in step S1.

[0146] In this embodiment, the low-frequency inertial component generated by multi-resolution signal decomposition and the dynamic zero-bias prediction value calculated based on multi-dimensional physical feature tensor both pass through a first-in-first-out sliding buffer with a length of sixty-four sampling points on the data flow path. Physically and logically, this means that the two signals have exactly the same physical time lag relative to the actual motion of the aircraft.

[0147] Assuming that the low-frequency inertial component and the dynamic zero-bias prediction value have the same phase lag property, the low-frequency inertial component is used as the minuend and the dynamic zero-bias prediction value is used as the subtrahend. Subtraction is performed on the minuend and the subtrahend to obtain a pure motion signal.

[0148] In this embodiment, the processor performs vector subtraction at each sampling moment. Based on the inverse operation logic of the linear superposition principle, since the nonlinear vibration rectification error is manifested as an additive DC bias superimposed on the real acceleration signal, the bias can be eliminated in the amplitude domain by performing subtraction.

[0149] The pure motion signal has eliminated the false DC component caused by the nonlinear vibration rectification effect in terms of amplitude, but it still has a physical lag relative to the actual motion of the body in terms of time.

[0150] Comparative experimental data show that performing subtraction directly without prior hysteresis alignment introduces a residual error of about 15% due to phase mismatch; however, by adopting the prior hysteresis alignment strategy of this embodiment, the residual error is reduced to less than 0.5%, verifying the key role of phase consistency in signal cleaning.

[0151] The specific operation for determining the intrinsic group delay time of the algorithm generated by the multi-resolution signal decomposition processing in step S1 is as follows: obtain the optimal decomposition depth and sensor sampling rate value determined in step S1, perform a power operation on the optimal decomposition depth to obtain the causal buffer cost value introduced by the first-in-first-out sliding buffer.

[0152] In this embodiment, the optimal decomposition depth is set to six layers based on the correspondence between the physical response cutoff frequency (e.g., 30 Hz) and the sensor sampling rate (e.g., 1 kHz). The processor calculates the number 2 raised to the power of 6 to obtain a causal buffer cost of 64. The physical meaning of this value is: in order to transform the lifting wavelet transform, which originally depends on future data, into a causal system, a waiting period equal to the power of two decomposition layers and sampling points must be introduced.

[0153] Obtain the filter tap length value of the lifting wavelet filter; subtract one from the filter tap length value to obtain the filter length difference value, and divide the filter length difference value by two to obtain the filter phase delay value. In this embodiment, the following is selected: The biorthogonal wavelet basis has a fixed filter tap length of five.

[0154] The processor calculates the difference between five and one, then divides it by two to obtain the filter phase delay value of two. This value precisely quantifies the inherent phase shift introduced by the digital filter during processing.

[0155] The causal buffer cost value and the filter phase delay value are added together to obtain the total number of delay sampling points. The total number of delay sampling points is then divided by the sensor sampling rate value to obtain the algorithm's intrinsic group delay time.

[0156] In this embodiment, the processor adds 64 to 2 to obtain a total number of delay sampling points of 66. Then, the processor divides 66 by 1000 (i.e., the sensor sampling rate value) to calculate the algorithm's intrinsic group delay time as 0.066 seconds. This intrinsic group delay time is a deterministic physical quantity determined by the algorithm structure, rather than an empirical estimate.

[0157] Compared to existing technologies that use fixed empirical values ​​(such as 0.05 seconds) for compensation, this embodiment eliminates compensation deviations caused by changes in sampling rate or decomposition level by using the intrinsic group delay time of the algorithm obtained through analytical calculation, thus ensuring the accuracy of the physical reference for subsequent time compensation.

[0158] In step S4, the first-order derivative feedforward compensation processing of the pure motion signal is performed using the inherent group delay time of the algorithm, and the specific operation of outputting the reconstructed positioning signal is as follows: a restricted Taylor reconstruction operator is constructed, and the restricted Taylor reconstruction operator performs differential operation on the pure motion signal using the three-point Lagrange backward difference algorithm to obtain the rate of change value of the pure motion signal.

[0159] In this embodiment, in order to extract smooth derivatives (i.e. acceleration) from discrete digital signals, the two-point difference algorithm is strictly prohibited to avoid amplifying high-frequency quantization noise.

[0160] The three-point Lagrange backward difference algorithm uses sample points from the current time, the previous time, and the two time points before that to calculate the slope of the current time through second-order polynomial fitting. The calculation formula is: multiply the current time value by three, subtract the previous time value multiplied by four, add the two time points before that, and finally divide by twice the sampling interval. This algorithm reduces the variance of differential noise by about 40% while ensuring that it does not use future data.

[0161] The maximum physical acceleration value of the aircraft is read from the non-volatile memory configured in the flight control computer, and the kinematic truncation process is performed on the rate of change value using the maximum physical acceleration value of the aircraft.

[0162] In this embodiment, the maximum physical acceleration of the aircraft is preferably set to 500 meters per cubic second. The selection of this parameter is based on bench testing of the aircraft's power system: due to the physical limitations of the motor torque response speed and rotational inertia, the rate of change of acceleration of a macroscopic rigid body cannot be infinite. Any rate of change exceeding this physical limit must be a non-physical artifact caused by sensor quantization noise or transient impact.

[0163] The specific steps of the kinematic truncation process are as follows: determine whether the absolute value of the rate of change is greater than the maximum physical acceleration of the aircraft. If the absolute value of the rate of change is greater than the maximum physical acceleration of the aircraft, set the rate of change to the maximum physical acceleration of the aircraft and output it as the safe acceleration value. If the absolute value of the rate of change is not greater than the maximum physical acceleration of the aircraft, keep the rate of change unchanged and output it as the safe acceleration value.

[0164] In this embodiment, the processor executes numerical clamping logic. Comparative experiments show that without kinematic truncation processing, differential compensation leads to a decrease of 3 dB in the signal-to-noise ratio. However, after introducing kinematic truncation processing, not only are high-frequency glitches eliminated, but the signal-to-noise ratio is also improved by 6 dB, demonstrating the necessity of physical constraints in signal reconstruction.

[0165] The safe acceleration value is used as the multiplicand, and the algorithm's intrinsic group delay time is used as the multiplier. The safe acceleration value and the algorithm's intrinsic group delay time are multiplied to obtain the time-domain feedforward compensation amount.

[0166] In this embodiment, the time-domain feedforward compensation amount physically represents the predicted value of the acceleration increment generated based on the current motion trend within the inherent group delay time of the delayed algorithm.

[0167] The pure motion signal is added to the time-domain feedforward compensation value to output a reconstructed positioning signal. In this embodiment, based on the first-order Taylor series expansion principle, the value at the lag time is added to the change within the lag time to approximate the true value at the current time.

[0168] The final reconstructed positioning signal not only eliminates vibration rectification error in amplitude but also compensates for the inherent phase lag of the algorithm in time, achieving real-time output with zero phase delay.

[0169] Experimental verification shows that the reconstructed positioning signal improves the phase margin of the flight control system by 15 degrees, and in high-dynamic flight tests, the position estimation error is reduced by 72% compared with the traditional low-pass filtering scheme.

[0170] Through amplitude cleaning and spatiotemporal reconstruction in this embodiment, the present invention completes the final closed loop of returning from mathematical calculation space to physical reality space. The output reconstructed positioning signal has high purity and can be directly and seamlessly called by the attitude calculation module of the flight controller. It effectively solves the positioning drift and control divergence problems caused by data distortion in inertial navigation systems under high vibration environment, and provides a solid data foundation for the safe flight of aircraft under complex conditions.

[0171] Example 2

[0172] like Figure 2 As shown, the present invention also discloses an aircraft positioning system, comprising: a signal decomposition module, a feature construction module, an error prediction module, and a signal reconstruction module, wherein;

[0173] The signal decomposition module is used to acquire the aircraft's body dynamic parameters, calculate the physical response cutoff frequency based on the body dynamic parameters, and perform multi-resolution signal decomposition processing on the raw inertial signals collected by the sensors using the physical response cutoff frequency to obtain the low-frequency inertial component that characterizes the true motion trend of the aircraft and the high-frequency vibration component that characterizes non-aircraft motion interference.

[0174] The feature construction module is used to perform mean square value calculation on the high-frequency vibration components obtained by the signal decomposition module to obtain the vibration energy density value. Simultaneously, it acquires the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. It performs time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor.

[0175] The error prediction module is used to input the multidimensional physical feature tensor constructed by the feature construction module into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained by the feature construction module to obtain the dynamic zero bias prediction value.

[0176] The signal reconstruction module is used to perform a subtraction operation on the low-frequency inertial component obtained by the signal decomposition module in the time domain using the dynamic zero-bias prediction value obtained by the error prediction module to obtain a clean motion signal. It determines the inherent group delay time of the algorithm generated by the multi-resolution signal decomposition processing in the signal decomposition module, performs first-order derivative feedforward compensation processing on the clean motion signal using the inherent group delay time of the algorithm, and outputs the reconstructed positioning signal.

[0177] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for locating an aircraft, characterized in that, Includes the following steps: Step S1: Obtain the aircraft's body dynamic parameters, calculate the physical response cutoff frequency based on the body dynamic parameters, and use the physical response cutoff frequency to perform multi-resolution signal decomposition processing on the original inertial signals collected by the sensors to obtain the low-frequency inertial component characterizing the true motion trend of the aircraft and the high-frequency vibration component characterizing non-aircraft motion disturbances. Step S2: Perform mean square value calculation on the high-frequency vibration components obtained in step S1 to obtain the vibration energy density value. Simultaneously acquire the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. Perform time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor. Step S3: Input the multidimensional physical feature tensor constructed in step S2 into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained in step S2 to obtain the dynamic zero bias prediction value. Step S4: In the time domain, use the dynamic zero-bias prediction value obtained in step S3 to perform a subtraction operation on the low-frequency inertial component obtained in step S1 to obtain a pure motion signal. Determine the algorithm's intrinsic group delay time generated by the multi-resolution signal decomposition processing in step S1. Use the algorithm's intrinsic group delay time to perform first-order derivative feedforward compensation processing on the pure motion signal to output a reconstructed positioning signal.

2. The aircraft positioning method according to claim 1, characterized in that, In step S1, the specific operation of calculating the physical response cutoff frequency based on the body's dynamic parameters is as follows: A10 reads pre-stored airframe dynamics parameters from the non-volatile memory configured in the flight control computer. The airframe dynamics parameters include the maximum resultant moment modulus generated by the aircraft's propulsion system and the airframe rotational inertia tensor of the aircraft's airframe structure. A11, perform matrix eigenvalue decomposition on the body moment of inertia tensor, and select the eigenvalue with the largest value from the result of the eigenvalue decomposition as the maximum principal moment of inertia; A12, obtain the control loop dead zone threshold of the flight control system, take the maximum combined torque modulus as the dividend, take the product of the maximum principal moment of inertia and the control loop dead zone threshold as the divisor, and perform division to obtain the torque-inertia ratio. A13, perform square root operation on the torque-moment of inertia ratio to obtain the angular velocity response limit value, divide the angular velocity response limit value by twice the value of pi to obtain the physical response cutoff frequency, and set the physical response cutoff frequency as the frequency domain boundary value for distinguishing the actual motion trend of the body from non-body motion interference during multi-resolution signal decomposition processing.

3. The aircraft positioning method according to claim 2, characterized in that, In step S1, the specific operation of performing multi-resolution signal decomposition processing on the raw inertial signal acquired by the sensor using the physical response cutoff frequency is as follows: A20, obtain the sensor sampling rate value when the sensor collects the original inertial signal, and divide the sensor sampling rate value by twice the physical response cutoff frequency to obtain the frequency ratio value. A21 performs a base-2 logarithmic operation on the frequency ratio value and rounds down the result to obtain the optimal decomposition depth for multi-resolution signal decomposition. A first-in-first-out sliding buffer is constructed with a correlation between the buffer depth and the optimal decomposition depth. The original inertial signal is pushed into the first-in-first-out sliding buffer in the order of acquisition time. A definite algorithm-inherent group delay time is formed during the process of pushing the original inertial signal into the first-in-first-out sliding buffer. A22 employs an integer-domain lifting wavelet transform algorithm within a first-in-first-out sliding buffer to perform splitting, prediction, and update processing on the original inertial signal. In the prediction processing, even-numbered position samples are used to perform linear prediction operations on odd-numbered position samples and calculate the prediction error to generate high-frequency detail coefficients. In the update processing, high-frequency detail coefficients are used to perform update operations on even-numbered position samples to generate low-frequency approximation coefficients. A23 repeatedly performs splitting, prediction, and update processing according to the optimal decomposition depth, outputs the low-frequency approximation coefficients of the last level as low-frequency inertial components, and performs inverse reconstruction and superposition processing on the high-frequency detail coefficients of all levels to obtain high-frequency vibration components.

4. The aircraft positioning method according to claim 1, characterized in that, In step S2, the mean square value calculation is performed on the high-frequency vibration components obtained in step S1 to obtain the vibration energy density value. The specific operation is as follows: B10 reads pre-stored sensor sensitivity scaling factors from the non-volatile memory configured in the flight control computer. Sensor sensitivity scaling factors are conversion coefficients used to restore digitally quantized signals to physical acceleration values. B11, determine the length of the energy integration window used to calculate the vibration energy, and within the time range defined by the length of the energy integration window, use the high-frequency vibration component obtained in step S1 as the multiplicand and the sensor sensitivity scaling factor as the multiplier to perform point-to-point multiplication to obtain the physical acceleration signal. B12 performs a squaring operation on the physical acceleration signal to obtain a squared signal. Within the length of the energy integration window, the squared signal is summed to obtain the total vibration energy. The number of sampling points contained within the length of the energy integration window is obtained. The total vibration energy is used as the dividend, and the number of sampling points is used as the divisor to perform a division operation to obtain the vibration energy density value.

5. The aircraft positioning method according to claim 4, characterized in that, In step S2, the excitation source state parameters of the aircraft's propulsion system and the environmental physical parameters of the sensors are acquired synchronously. The vibration energy density value, excitation source state parameters, and environmental physical parameters are then subjected to time alignment and spatial stitching to construct a multidimensional physical feature tensor. The specific operations are as follows: B20 obtains motor speed data as excitation source status parameters through the communication bus of the aircraft's electronic speed controller; B21 uses the temperature sensor built into the inertial measurement unit to obtain the core temperature data of the sensor as environmental physical parameters. For asynchronous sampling conditions where the update frequency of the excitation source state parameters is lower than the update frequency of the vibration energy density value and the update frequency of the environmental physical parameters is lower than the update frequency of the vibration energy density value, a zero-order hold strategy is adopted to perform time alignment processing on the excitation source state parameters and the environmental physical parameters. The B22's zero-order hold strategy is to use the most recently updated excitation source state parameters and environmental physical parameters at the current moment, and read the sensor's nominal resonant frequency, resonant frequency temperature coefficient, and damping bandwidth half-width from the non-volatile memory configured in the flight control computer. B23 uses time-aligned environmental physical parameters to perform a correction operation on the nominal resonant frequency to obtain the actual natural frequency at the current temperature; B24. Calculate the multi-order harmonic frequencies using the time-aligned excitation source state parameters. Combine the actual natural frequency, multi-order harmonic frequencies, and damping bandwidth half-width to perform the calculation and obtain the sum of Lorentz linear resonance factors. B25 involves splicing the vibration energy density value, the time-aligned excitation source state parameters, the sum of the Lorentz linear resonance factor, and the time-aligned environmental physical parameters in a predetermined dimensional order to construct a multidimensional physical feature tensor.

6. The aircraft positioning method according to claim 1, characterized in that, In step S3, the multidimensional physical feature tensor constructed in step S2 is input into the pre-configured nonlinear error prediction model, and the specific operation of outputting the dimensionless rectifier coupling coefficient through the nonlinear error prediction model is as follows: C10 reads a pre-stored table of system physical limit parameters from the non-volatile memory configured in the flight control computer; C11, the system physical limit parameter table includes the maximum theoretical energy density of the sensor, the maximum physical speed of the power system and the maximum range of ambient temperature variation. The physical boundary scaling diagonal matrix and the physical reference center vector are constructed using the maximum theoretical energy density of the sensor, the maximum physical speed of the power system and the maximum range of ambient temperature variation. C12, the physical reference center vector performs centering processing on the environmental physical parameters, and the physical reference center vector keeps the absolute zero value of the sum of the vibration energy density value, the excitation source state parameters and the Lorentz linear resonance factor unchanged; C13, the diagonal elements of the physical boundary scaling diagonal matrix are the reciprocal of the sensor's maximum theoretical energy density, the reciprocal of the power system's maximum physical speed, the unit value, and the reciprocal of the maximum range of ambient temperature variation; C14 performs a linear manifold mapping operation based on physical limits on the multidimensional physical feature tensor using the physical boundary scaling diagonal matrix and the physical reference center vector to generate a dimensionless normalized input tensor; the dimensionless normalized input tensor is then input into a nonlinear error prediction model, which adopts a gated recurrent unit network structure with memory units. C15 uses the hyperbolic tangent activation function configured in the nonlinear error prediction model to perform amplitude limiting on the calculation results of the nonlinear error prediction model, and obtains the output value of the hyperbolic tangent activation function. C16 reads the pre-stored maximum rectified gain limit from the non-volatile memory configured in the flight control computer, performs a multiplication operation on the output value of the hyperbolic tangent activation function and the maximum rectified gain limit, and obtains the rectified coupling coefficient.

7. The aircraft positioning method according to claim 6, characterized in that, Read the pre-stored corresponding axial vibration rectification error sensitivity reference from the non-volatile memory configured in the flight control computer; C20 corresponds to the axial vibration rectification error sensitivity benchmark, which is the unit energy error conversion rate of the target axial direction at the most sensitive frequency point, measured by frequency sweep bench test during the factory calibration stage. C21, corresponding to the axial vibration rectification error sensitivity benchmark, has the physical dimension balancing function of restoring the square energy term to the first acceleration term; C22, use the dimensionless rectifier coupling coefficient output in step S3 as the first multiplier; C23, using the corresponding axial vibration rectification error sensitivity benchmark as the second multiplier; C24, use the vibration energy density value obtained in step S2 as the third multiplier; C25 performs a three-term multiplication operation on the first, second, and third multipliers to obtain the dynamic zero-bias prediction value. The three-term multiplication operation forces the dynamic zero-bias prediction value to be zero when the vibration energy density is zero.

8. The aircraft positioning method according to claim 1, characterized in that, In step S4, the specific operation of subtracting the low-frequency inertial component obtained in step S1 from the dynamic zero-bias prediction value obtained in step S3 in the time domain to obtain the pure motion signal is as follows: D10, confirm that the low-frequency inertial component obtained in step S1 and the dynamic zero-bias prediction value obtained in step S3 are both constrained by the first-in-first-out sliding buffer mechanism in step S1. Under the premise that the low-frequency inertial component and the dynamic zero-bias prediction value have the same phase lag property, take the low-frequency inertial component as the minuend and the dynamic zero-bias prediction value as the subtrahend, and perform a subtraction operation on the minuend and the subtrahend to obtain a pure motion signal. D11, the specific operation to determine the inherent group delay time of the algorithm generated by the multi-resolution signal decomposition processing in step S1 is as follows: obtain the optimal decomposition depth and sensor sampling rate value determined in step S1, perform a power operation on the optimal decomposition depth to obtain the causal buffer cost value introduced by the first-in-first-out sliding buffer, and obtain the filter tap length value of the lifting wavelet filter. D12, subtract one from the filter tap length value to get the filter length difference value, divide the filter length difference by two to get the filter phase delay value; D13 adds the causal buffer cost value to the filter phase delay value to obtain the total delay sampling points. Dividing the total delay sampling points by the sensor sampling rate value yields the algorithm's intrinsic group delay time.

9. The aircraft positioning method according to claim 8, characterized in that, In step S4, the specific operation of performing first-order derivative feedforward compensation processing on the pure motion signal using the algorithm's intrinsic group delay time to output the reconstructed positioning signal is as follows: D20, construct a restricted Taylor reconstruction operator. The restricted Taylor reconstruction operator uses the three-point Lagrange backward difference algorithm to perform differential operations on the pure motion signal to obtain the rate of change value of the pure motion signal. D21 reads the pre-stored maximum physical acceleration value of the aircraft from the non-volatile memory configured in the flight control computer, and performs kinematic truncation processing on the rate of change value using the maximum physical acceleration value of the aircraft. D22, the specific steps of the kinematic truncation process are as follows: determine whether the absolute value of the rate of change is greater than the maximum physical acceleration value of the aircraft. If the absolute value of the rate of change is greater than the maximum physical acceleration value of the aircraft, then set the rate of change to the maximum physical acceleration value of the aircraft and output it as the safe acceleration value. If the absolute value of the rate of change is not greater than the maximum physical acceleration value of the aircraft, then keep the rate of change unchanged and output it as the safe acceleration value. D23, using the safety acceleration value as the multiplicand and the algorithm's intrinsic group delay time as the multiplier, performs a multiplication operation on the safety acceleration value and the algorithm's intrinsic group delay time to obtain the time-domain feedforward compensation amount; D24 performs an addition operation on the pure motion signal and the time-domain feedforward compensation, and outputs the reconstructed positioning signal.

10. An aircraft positioning system, employing an aircraft positioning method as described in any one of claims 1-9, characterized in that, include: The signal decomposition module, feature construction module, error prediction module, and signal reconstruction module are included. The signal decomposition module is used to acquire the aircraft's body dynamic parameters, calculate the physical response cutoff frequency based on the body dynamic parameters, and perform multi-resolution signal decomposition processing on the raw inertial signals collected by the sensors using the physical response cutoff frequency to obtain the low-frequency inertial component that characterizes the true motion trend of the aircraft and the high-frequency vibration component that characterizes non-aircraft motion interference. The feature construction module is used to perform mean square value calculation on the high-frequency vibration components obtained by the signal decomposition module to obtain the vibration energy density value. Simultaneously, it acquires the excitation source state parameters of the aircraft power system and the environmental physical parameters of the sensor. It performs time alignment and spatial stitching processing on the vibration energy density value, excitation source state parameters and environmental physical parameters to construct a multidimensional physical feature tensor. The error prediction module is used to input the multidimensional physical feature tensor constructed by the feature construction module into the pre-configured nonlinear error prediction model, output the dimensionless rectification coupling coefficient through the nonlinear error prediction model, and perform a product operation on the rectification coupling coefficient and the vibration energy density value obtained by the feature construction module to obtain the dynamic zero bias prediction value. The signal reconstruction module is used to perform a subtraction operation on the low-frequency inertial component obtained by the signal decomposition module in the time domain using the dynamic zero-bias prediction value obtained by the error prediction module to obtain a clean motion signal. It determines the inherent group delay time of the algorithm generated by the multi-resolution signal decomposition processing in the signal decomposition module, performs first-order derivative feedforward compensation processing on the clean motion signal using the inherent group delay time of the algorithm, and outputs the reconstructed positioning signal.

Images