Satellite on-orbit anomaly detection system based on AI dynamic behavior analysis

Abstract

The present application relates to the technical field of fault prediction, in particular to a satellite on-orbit anomaly detection system based on AI dynamic behavior analysis, which comprises: an attitude signal mode decomposition module, used for collecting attitude angular velocity feedback signals of a satellite attitude control system and flywheel control instruction sequences, extracting high-frequency flutter components of flexible accessories, screening the high-frequency flutter components of the flexible accessories according to signal energy proportion, retaining main vibration modes, and generating a key flexible mode eigenfunction sequence.In the present application, a multi-dimensional on-orbit abnormal state characteristic vector is constructed in combination with an entropy value drop slope feature, and is projected into a fault mode topological space for matching, realizing a leap from single threshold judgment to multi-dimensional feature fusion diagnosis, which can not only distinguish structural resonance from control loop transmission delay, but also provide a classification result containing delay state and stability features in the early stage of fault, thereby improving the accuracy of on-orbit anomaly detection.

Description

Technical Field

[0001] This invention relates to the field of fault prediction technology, and in particular to a satellite on-orbit anomaly detection system based on AI dynamic behavior analysis. Background Technology

[0002] Current technologies for monitoring satellite on-orbit operation typically rely solely on static threshold determinations of telemetry parameters or simple statistical variance analysis, neglecting the complex dynamic coupling effects between the satellite's flexible appendages and attitude control system. When faced with high-frequency, low-amplitude flutter signals, the lack of effective signal mode separation means that critical structural vibration characteristics are easily masked by large-amplitude rigid body attitude maneuvering signals, making it difficult to identify fatigue damage risks in flexible structures at an early stage. Furthermore, alarms are only triggered when the vibration amplitude diverges to a destructive level, resulting in delayed fault detection. Therefore, improvements are needed. Summary of the Invention

[0003] The purpose of this invention is to address the problem that key structural vibration characteristics are easily overwhelmed by large-amplitude rigid body attitude maneuvering signals, making it difficult to identify fatigue damage risks in flexible structures in the early stages. The proposed invention is a satellite on-orbit anomaly detection system based on AI dynamic behavior analysis.

[0004] To achieve the above objectives, the present invention adopts the following technical solution: A satellite on-orbit anomaly detection system based on AI dynamic behavior analysis includes:

[0005] The attitude signal mode decomposition module is used to acquire the attitude angular velocity feedback signal and flywheel control command sequence of the satellite attitude control system, extract the high-frequency flutter component of the flexible accessory, filter the high-frequency flutter component of the flexible accessory according to the signal energy ratio, retain the main vibration mode, and generate the key flexible mode eigenfunction sequence.

[0006] The vibration entropy trend analysis module is used to divide the key flexible mode eigenfunction sequence into continuously sliding time slices, calculate the permutation entropy value of the signal amplitude sequence in each time slice, obtain the vibration complexity time-varying entropy value sequence, identify the time interval in which the permutation entropy value shows a continuous downward trend based on the vibration complexity time-varying entropy value sequence, and generate a structural coupling instability early warning time window.

[0007] The command response cross-correlation module is used to extract the flywheel control command sequence and attitude angular velocity feedback signal segment corresponding to the structural coupling instability early warning time window, calculate the cross-correlation function value, search for the time delay corresponding to the maximum point of cross-correlation function amplitude, generate command response peak offset, and calculate the control loop response synchronization lag fraction based on the command response peak offset.

[0008] The fault vector mapping determination module is used to compare the control loop response synchronization lag score with a preset processor bus blocking threshold, combine the descent slope characteristics within the structural coupling instability early warning time window, construct an on-orbit abnormal state feature vector, and generate satellite on-orbit abnormality detection classification results based on the on-orbit abnormal state feature vector.

[0009] Preferably, the steps for obtaining the key flexible mode eigenfunction sequence are as follows:

[0010] Based on the attitude angular velocity feedback signal of the satellite attitude control system, the attitude angular velocity feedback signal is frequency domain partitioned according to the preset bandwidth limit condition. The center frequency values ​​of different frequency bands in the attitude angular velocity feedback signal are iteratively updated until the center frequency values ​​converge to a stable state. The attitude angular velocity feedback signal is decomposed into multiple intrinsic mode function components of different frequency bands. The signal part with a frequency value greater than the preset rigid body cutoff frequency value in the intrinsic mode function components is extracted. The low frequency components of rigid body motion with a frequency value less than or equal to the rigid body cutoff frequency value are removed to generate high frequency flutter components of flexible attachments.

[0011] The instantaneous amplitude square term of each modal component in the high-frequency flutter component of the flexible attachment is calculated. The instantaneous amplitude square term is integrated and accumulated in the time domain to obtain the absolute energy value of each modal component. The absolute energy value is divided by the sum of the absolute energy values ​​of all modal components to calculate the weight ratio of each modal component in the overall vibration and generate the vibration mode energy ratio value.

[0012] The vibration mode energy ratio is compared with a preset primary mode screening threshold. Modal components whose vibration mode energy ratio is greater than the primary mode screening threshold are identified. These modal components are marked as primary vibration modes. Secondary modes whose vibration mode energy ratio is less than or equal to the primary mode screening threshold are removed. The time-domain waveform data corresponding to the retained primary vibration modes are recombined to generate a sequence of key flexible mode eigenfunctions.

[0013] Preferably, the steps for obtaining the time-varying entropy value sequence of vibration complexity are as follows:

[0014] Set a fixed length of time window width and movement step size, continuously extract the sequence of eigenfunctions of the key flexible modes along the time axis, obtain multiple local time slices, reconstruct the phase space of the signal amplitude data in each time slice and count the occurrence probability of different symbol modes, calculate the quantitative index characterizing the disorder of the system based on Shannon information theory, and generate a time-varying entropy value sequence of vibration complexity.

[0015] Calculate the first-order difference value based on the time-varying entropy value sequence of the vibration complexity.

[0016] Preferably, the step of obtaining the structural coupling instability early warning time window is as follows:

[0017] The first-order difference value is compared with a preset negative change threshold to filter out anomalies where the first-order difference value is less than the negative change threshold. The start and end positions of the continuous occurrence of anomalies on the time axis are identified, the duration of the continuous anomalies is calculated and it is determined whether it exceeds the minimum instability confirmation period. The time period covered by the start and end positions that meet the conditions is marked as the risk interval, and a structural coupling instability early warning time window is generated.

[0018] Preferably, the step of obtaining the peak offset of the instruction response is as follows:

[0019] Using the structural coupling instability early warning time window as the interception benchmark, discrete data segments corresponding to the time period are extracted from the flywheel control command sequence and the attitude angular velocity feedback signal, respectively. Cross-correlation convolution operation is performed on the two discrete data segments to generate a cross-correlation function numerical sequence.

[0020] Based on the numerical sequence of the cross-correlation function, the point with the largest amplitude in the sequence is located by the peak search algorithm. The distance between the time coordinate corresponding to the maximum value point and the zero time is determined, and the peak offset of the command response is generated. At the same time, the span of the peak waveform on the time axis is calculated with half of the maximum amplitude as the boundary, and the full width at half maximum (FWHM) value of the cross-correlation peak point is generated.

[0021] Preferably, the step of obtaining the synchronization lag fraction of the control loop response is as follows:

[0022] The control loop response lag score is calculated based on the peak offset of the command response and the full width at half maximum (FWHM) value of the cross-correlation peak point.

[0023] Preferably, the step of obtaining the on-orbit abnormal state feature vector is as follows:

[0024] The control loop response lag score is compared with a preset processor bus blocking threshold. If the control loop response lag score is greater than the processor bus blocking threshold, it is determined to be a signal blocking delay state. If the control loop response lag score is less than or equal to the processor bus blocking threshold, it is determined to be a signal transmission smooth state, generating a signal transmission delay state for the control loop. At the same time, the structural coupling instability early warning time window is retrieved, and all vibration complexity time-varying entropy values ​​within the structural coupling instability early warning time window are extracted. A linear regression algorithm is used to fit the trend line of the vibration complexity time-varying entropy value, and the slope value of the trend line is calculated to generate an entropy value decrease slope feature.

[0025] Based on the signal transmission delay state of the control loop and the entropy decrease slope feature, a multi-dimensional feature mapping coordinate system is established. The signal transmission delay state of the control loop is quantized into coordinate values ​​of the first dimension, and the entropy decrease slope feature is quantized into coordinate values ​​of the second dimension. The coordinate values ​​of the first dimension and the coordinate values ​​of the second dimension are combined in sequence to generate an on-orbit anomaly state feature vector.

[0026] Preferably, the steps for obtaining the satellite's on-orbit anomaly detection and classification results are as follows:

[0027] The on-orbit anomaly feature vector is projected into the pre-constructed fault mode topology space. The Euclidean distance between the on-orbit anomaly feature vector and each standard fault node in the fault mode topology space is calculated. The standard fault node with the smallest Euclidean distance is searched as the best matching node. The fault type description information and fault level label bound to the best matching node are read to generate the satellite on-orbit anomaly detection and classification results.

[0028] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0029] In this invention, the attitude angular velocity feedback signal and flywheel control command sequence of the satellite attitude control system are acquired, and a variational mode decomposition operation with adaptive optimization of the center frequency is performed. This can separate the low-frequency components of rigid body motion from the high-frequency flutter components of flexible attachments, eliminate the masking effect of rigid body attitude maneuvering on weak flutter signals, and retain the main vibration modes through a screening operation based on the signal energy ratio, ensuring that subsequent analysis focuses on the key structural features with destructive potential. The eigenfunction sequence of key flexible modes is divided into continuously sliding time slices, and the permutation entropy value is calculated, which can quantify the complexity change of the signal on the time axis and capture the continuous decrease of entropy value as the system transforms from disordered random vibration to ordered self-excited oscillation. The system generates a structural coupling instability early warning time window, enabling early location of flutter instability precursors. Within this time window, signal segments are extracted and cross-correlation functions and command response peak offsets are calculated. This decouples the temporal correlation between control commands and physical responses, quantifies the lag fraction of control loop response synchronization, and constructs a multi-dimensional on-orbit anomaly feature vector by combining entropy decrease slope characteristics. This vector is then projected onto the fault mode topology space for matching, achieving a leap from single threshold judgment to multi-dimensional feature fusion diagnosis. It can not only distinguish between structural resonance and control loop transmission delay, but also provide classification results including delay state and stability features in the early stages of a fault, improving the accuracy of on-orbit anomaly detection. Attached Figure Description

[0030] Figure 1 This is a schematic diagram of the principle of the present invention. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0032] Please see Figure 1 This invention provides a technical solution: a satellite on-orbit anomaly detection system based on AI dynamic behavior analysis, comprising:

[0033] The attitude signal mode decomposition module is used to acquire the attitude angular velocity feedback signal and flywheel control command sequence of the satellite attitude control system, extract the high-frequency flutter components of the flexible appendages, filter the high-frequency flutter components of the flexible appendages according to the proportion of signal energy, retain the main vibration modes, and generate the key flexible mode eigenfunction sequence.

[0034] The vibration entropy trend analysis module is used to divide the sequence of eigenfunctions of key flexible modes into continuously sliding time slices, calculate the permutation entropy of the signal amplitude sequence in each time slice, obtain the time-varying entropy sequence of vibration complexity, identify the time interval in which the permutation entropy value shows a continuous downward trend based on the time-varying entropy sequence of vibration complexity, and generate a structural coupling instability early warning time window.

[0035] The command response cross-correlation module is used to extract the flywheel control command sequence and attitude angular velocity feedback signal segment corresponding to the structural coupling instability early warning time window, calculate the cross-correlation function value, search for the time delay corresponding to the maximum point of cross-correlation function amplitude, generate command response peak offset, and calculate the control loop response synchronization lag fraction based on the command response peak offset.

[0036] The fault vector mapping and determination module is used to compare the control loop response lag score with the preset processor bus blocking threshold, combine the descent slope characteristics within the structural coupling instability early warning time window, construct the on-orbit abnormal state feature vector, and generate the satellite on-orbit abnormality detection classification result based on the on-orbit abnormal state feature vector.

[0037] The steps for obtaining the sequence of eigenfunctions of key flexible modes are as follows:

[0038] Based on the attitude angular velocity feedback signal of the satellite attitude control system, the attitude angular velocity feedback signal is frequency domain partitioned according to the preset bandwidth limit condition. The center frequency values ​​of different frequency bands in the attitude angular velocity feedback signal are iteratively updated until the center frequency values ​​converge to a stable state. The attitude angular velocity feedback signal is decomposed into multiple intrinsic mode function components of different frequency bands. The signal part with a frequency value greater than the preset rigid body cutoff frequency value in the intrinsic mode function components is extracted. The low frequency components of rigid body motion with a frequency value less than or equal to the rigid body cutoff frequency value are removed to generate high frequency flutter components of flexible attachments.

[0039] Calculate the instantaneous amplitude squared term of each modal component in the high-frequency flutter component of the flexible attachment, perform an integral and cumulative operation on the instantaneous amplitude squared term in the time domain to obtain the absolute energy value of each modal component, divide the absolute energy value by the sum of the absolute energy values ​​of all modal components, calculate the weight ratio of each modal component in the overall vibration, and generate the vibration mode energy ratio value.

[0040] The vibration mode energy ratio is compared with the preset main mode screening threshold. The modal components whose vibration mode energy ratio is greater than the main mode screening threshold are identified. The modal components whose vibration mode energy ratio is greater than the main mode screening threshold are marked as main vibration modes. The secondary modes whose vibration mode energy ratio is less than or equal to the main mode screening threshold are removed. The time domain waveform data corresponding to the retained main vibration modes are recombined to generate the key flexible mode eigenfunction sequence.

[0041] Specifically, based on the attitude angular velocity feedback signal from the satellite attitude control system, the sampling frequency for signal acquisition is first set to 1000Hz. This sampling frequency setting is based on the Nyquist sampling theorem and can effectively cover the frequency domain range from 0Hz to 500Hz. The acquired attitude angular velocity feedback signal is then input into the variational mode decomposition algorithm to initialize the number of mode decompositions. The value is set to 8, which is based on the theoretical number of modes of the satellite and its accessories, and serves as a penalty factor. To ensure good bandwidth compactness of the decomposed modes, a variational constraint problem is constructed, aiming to minimize the sum of the estimated bandwidths of each mode, by introducing Lagrange multipliers. The constrained variational problem is transformed into an unconstrained problem, and the problem is solved iteratively using the alternating direction multiplier method, updating the modal components in each iteration. and center frequency The formula for updating the center frequency is: ,in Representing the In the nth iteration The center frequency of each mode For frequency variables, For the first After the first iteration update The Fourier transform results of each mode are calculated, the Lagrange multipliers are updated, the convergence condition is checked, the sum of the square norms of the differences between the modal components in the current iteration and the previous iteration is calculated, and it is determined whether the sum is less than the preset convergence tolerance. ,in Set as If the condition is met, the iteration stops, and the result is obtained. Each intrinsic mode function component is identified. The rigid body cutoff frequency is then set to 0.5 Hz, determined based on the upper limit of the closed-loop bandwidth of the satellite attitude control loop. All decomposed intrinsic mode function components are iterated over to extract the center frequency. Signal frequencies greater than 0.5Hz are identified as flexible vibrations, and the center frequency is discarded. Low-frequency components less than or equal to 0.5Hz, which usually correspond to large-angle maneuvers or orbital motions of the satellite body, are used to generate high-frequency flutter components of flexible attachments by time-domain superposition of all selected high-frequency modal components or by retaining their independent channels.

[0042] Calculate the instantaneous amplitude squared term of each modal component in the high-frequency flutter component of the flexible attachment, and apply this to each high-frequency flutter modal component obtained in the previous step. Perform Hilbert transform to construct analytic signal The calculation formula is: ,in It is an analytical signal. These are the original modal components. It is the imaginary unit. The Hilbert transform operation is defined as the interaction between the original signal and... The convolution is used to calculate the instantaneous amplitude based on the analytic signal. The calculation formula is: After obtaining the instantaneous amplitude at each moment, square it to obtain... This represents the instantaneous power of the mode at each time step. Then, the squared instantaneous amplitude term is integrated and accumulated in the time domain, with the integration time interval set as... , The value is set to the length of the current analysis window, for example, 10 seconds. The absolute energy value of each modal component is calculated using the following formula: ,in Representing the The absolute energy of each modal component, Using time as the variable, count the number of all retained high-frequency modal components. The total energy is obtained by summing the absolute energy values ​​of all modal components. The weight ratio of each modal component in the overall vibration is calculated using division. The formula is as follows: ,in For the first The energy percentage of each mode reflects the contribution of the vibration in that frequency band to the overall flexible flutter. The sum of the weights of all modes is 1, and the energy percentage of the vibration modes is generated.

[0043] The vibration mode energy ratio is compared with the preset main mode screening threshold, and the main mode screening threshold is set accordingly. The threshold is set at 0.10. This threshold is derived from statistical analysis of historical on-orbit micro-vibration data. Modes with an energy percentage below 10% are usually sensor noise or non-structural transient interference. The threshold is calculated by iterating through the vibration mode energy percentage of each modal component. Execute conditional judgment logic, if If the modal component contains major structural coupling information, it is marked as the main vibration mode and its corresponding time-domain waveform data is retained. If a mode is identified as a minor mode or background noise, it is removed to reduce computational redundancy in subsequent entropy analysis. After traversing and filtering all modes, the number of retained principal vibration modes is checked. If no mode meets the threshold condition, the mode with the largest energy percentage is forcibly retained as the principal vibration mode to prevent complete signal loss. Time-domain waveform data of all marked principal vibration modes are then extracted. These waveform data are linearly superimposed and recombined; the calculation formula is as follows: ,in The recombined signal sequence, The total number of main vibration modes to be retained. For the first The principal vibration modes are in The amplitude at each moment is used to reconstruct the vibration signal containing only the main flexible features through this superposition process, generating a sequence of key flexible mode eigenfunctions.

[0044] The steps for obtaining the time-varying entropy value sequence of vibration complexity are as follows:

[0045] Set a fixed-length time sliding window width and a moving step size, continuously extract the sequence of key flexible mode eigenfunctions along the time axis, obtain multiple local time slices, reconstruct the phase space of the signal amplitude data in each time slice and count the occurrence probability of different symbol modes, calculate the quantitative index characterizing the disorder of the system based on Shannon information theory, and generate a time-varying entropy value sequence of vibration complexity.

[0046] Based on the time-varying entropy sequence of vibration complexity, the first-order difference value is calculated using the following formula:

[0047] ;

[0048] in, Let be the first-order difference value at time k. Let the time-varying entropy value be the vibration complexity at time k. Let the time-varying entropy value be the vibration complexity at time k-1. The preset downward trend sensitivity index, This is a minimal positive compensation term used to prevent singularity in the denominator.

[0049] Specifically, a fixed-length time window width and movement step size are set. Based on the 1000Hz sampling frequency and 0.5Hz cutoff frequency set in the previous steps, to ensure that each time slice contains at least two complete lowest-frequency vibration cycles, the time window width is set to 4000 sampling points (4 seconds) and the movement step size is set to 1000 sampling points (1 second). The key flexible mode eigenfunction sequence is continuously truncated along the time axis to obtain multiple local time slices. The signal amplitude data in each time slice is reconstructed in phase space. The embedding dimension is set to 4 to adapt to the high-dimensional flexible vibration characteristics, and the delay time is set to 25 sampling points to eliminate sequence autocorrelation. The probability of occurrence of different symbol modes is statistically analyzed. By sorting and encoding the numerical relationship of the elements inside the reconstructed vector, the proportion of each permutation mode in the total number of modes is calculated. A quantitative index characterizing the disorder of the system is calculated based on Shannon information theory. When the vibration is in a random disordered state, the index approaches 1. When coupling instability symptoms appear, causing the vibration to show regular divergence, the index value decreases significantly, generating a time-varying entropy value sequence of vibration complexity.

[0050] In the first-order difference numerical calculation formula, a nonlinear exponential weighting term based on the ratio of past to current entropy values ​​is introduced. When the entropy value shows a downward trend, that is, when the system tends to become orderly unstable, the difference signal is amplified, thereby improving the identification of early weak instability characteristics.

[0051] The steps for obtaining the parameters are as follows: read the current value from the generated time-varying entropy value sequence of vibration complexity. The value at each monitoring moment, such as the entropy value of 0.8520 during the stable phase after the satellite maneuver and orbit change, represents the complexity level of the structural vibration at the current moment.

[0052] The steps to obtain the parameters are as follows: extract the first parameter from the sequence. The entropy value data at the previous moment, i.e., the previous second, was read as 0.8650 and used to establish a differential comparison benchmark.

[0053] The steps to obtain the parameter are as follows, and the value is... ;

[0054] The steps for obtaining the parameter are as follows: This parameter is the decreasing trend sensitivity index, which is determined based on the frequency ratio of the structural fundamental frequency to the control bandwidth and the environmental noise level. The calculation formula is as follows: ,in The first-order fundamental frequency of the satellite's flexible accessory was measured to be 2.5Hz through ground-based modal testing. The upper limit of the attitude control bandwidth is set to 0.5Hz;

[0055] The steps for obtaining the parameter are as follows: This parameter is the environmental noise correction coefficient, used to eliminate the influence of background noise on sensitivity. The calculation formula is... ,in The collected background noise amplitude sequence, For the sample size, This represents the average value of the collected background noise amplitude sequence. Using the reference signal amplitude, and taking the background noise variance calculation result of 0.04, and the reference amplitude of 0.1, the calculation is as follows: ;

[0056] Substitute the above parameters The calculation formula yields: ;

[0057] Calculations based on parameters:

[0058] First, calculate the linear difference component of the entropy value:

[0059] ;

[0060] Calculate the base of the weighted ratio:

[0061] ;

[0062] Calculate the index weighting term:

[0063] ;

[0064] Calculate the final first-order difference value:

[0065] ;

[0066] The results indicate that the entropy value at the current moment shows a downward trend, and after nonlinear weighting, the magnitude of this decrease is quantified as -0.0134. The negative value and the increase in absolute value indicate that the system is accelerating towards a structural coupling instability state, increasing the risk of triggering the early warning mechanism.

[0067] The steps for obtaining the early warning time window for structural coupling instability are as follows:

[0068] The first-order difference value is compared with a preset negative change threshold to filter out anomalies where the first-order difference value is less than the negative change threshold. The start and end positions of the continuous occurrence of anomalies on the time axis are identified, the duration of the continuous anomalies is calculated and it is determined whether it exceeds the minimum instability confirmation period. The time period covered by the start and end positions that meet the conditions is marked as the risk interval, and a structural coupling instability early warning time window is generated.

[0069] Specifically, the first-order difference value is compared with a preset negative change threshold, which is set based on the statistical distribution of historical stable operation data. The mean and standard deviation of the first-order difference of historical data are calculated, and the mean minus three times the standard deviation is taken as the threshold. For example, if the calculated value is -0.0085, outliers with first-order difference values ​​less than -0.0085 are selected. For example, if the currently calculated value is -0.0134, which is less than the threshold, it is judged as an anomaly. The start and end positions of the continuous occurrence of anomalies on the time axis are identified, and the start and end timestamps of the continuous occurrence of anomalies are recorded. The duration of the continuous anomalies is calculated and it is determined whether it exceeds the minimum instability confirmation period. This period is set to 8 seconds based on the response delay of the control system. If the duration of the continuous anomaly is detected to reach 12 seconds, which exceeds the 8-second confirmation threshold, the possibility of occasional interference is excluded. The time period covered by the start and end positions that meet the conditions is marked as the risk interval, and a structural coupling instability early warning time window is generated.

[0070] The steps for obtaining the peak offset of the instruction response are as follows:

[0071] Using the structural coupling instability early warning time window as the interception benchmark, discrete data segments corresponding to the time period are extracted from the flywheel control command sequence and the attitude angular velocity feedback signal, respectively. Cross-correlation convolution operation is performed on the two discrete data segments to generate a numerical sequence of cross-correlation function.

[0072] Based on the numerical sequence of cross-correlation functions, the point with the largest amplitude in the sequence is located by the peak search algorithm. The distance between the time coordinate corresponding to the maximum value point and the zero time is determined, and the peak offset of the command response is generated. At the same time, the span of the peak waveform on the time axis is calculated with half of the maximum amplitude as the boundary, and the full width at half maximum (FWHM) value of the cross-correlation peak point is generated.

[0073] Specifically, using the structural coupling instability early warning time window as the cutoff benchmark, and based on the risk interval determined in the previous steps, the corresponding flywheel control command sequence and attitude angular velocity feedback signal are extracted. To eliminate the phase error caused by asynchronous sampling, cubic spline interpolation is used to resample the two discrete sequences onto a unified high-frequency time axis. The two sequences are then de-DC processed, i.e., their respective time averages are subtracted to eliminate the influence of static bias on correlation analysis. The scanning range of the cross-correlation operation is set. Considering the physical transmission characteristics of the control loop, the maximum lag scanning time is set to 0.5 seconds, corresponding to 500 sampling points at a sampling frequency of 1000Hz. Discrete cross-correlation convolution operation is performed. Keeping the control command sequence stationary, the attitude angular velocity feedback signal is shifted point by point on the time axis. The product integral of the two signal amplitudes at each shift point is calculated to generate a numerical sequence of cross-correlation functions that reflects the change of signal similarity with delay time.

[0074] Based on the numerical sequence of cross-correlation functions, a sliding window extreme value search algorithm is used to traverse the cross-correlation function sequence to find the global maximum point of the correlation coefficient magnitude. The time index position of this maximum point in the sequence is locked. The index value is multiplied by the sampling interval of 0.001 seconds to calculate the time difference between the maximum correlation point and the zero lag time. For example, if the calculated lag time is 0.045 seconds, the peak offset of the command response is generated. At the same time, 50% of the maximum peak amplitude, i.e., half, is used as the cutoff threshold. The moment when the signal strength first decays to this threshold is searched on the left and right sides of the peak point, respectively. The time difference between the two half-power points is calculated. This difference represents the degree of diffusion of the response signal energy on the time axis. For example, if the waveform width is statistically 0.080 seconds, the full width at half maximum (FWHM) value of the cross-correlation peak point is generated.

[0075] The steps for obtaining the control loop response lag fraction are as follows:

[0076] Based on the peak offset of the command response and the full width at half maximum (FWHM) values ​​of the cross-correlation peak points, the synchronization lag fraction of the control loop response is calculated using the following formula:

[0077] ;

[0078] in, The control loop response lag fraction, This represents the peak offset of the command response. This is the preset system nominal delay value. The full width at half maximum (FWHM) value of the cross-correlation peak point. This is a preset waveform dispersion penalty factor.

[0079] Specifically, in the formula for calculating the lag fraction of the control loop response synchronization, the simple time transmission delay and the waveform distortion caused by flexible vibration are unified into an equivalent time lag. The second term utilizes the ratio of the cube of the waveform width to the transmission time to construct a diffuse energy term with the dimension of the square of time. Adjust its weight in the overall failure assessment;

[0080] The steps to obtain the parameters are as follows: read the peak offset of the command response generated by the previous steps. This value represents the actual monitored signal transmission loop delay time, for example, the current calculated value is 0.045 seconds.

[0081] The steps for obtaining the parameter are as follows: the parameter is a preset nominal delay value of the system, which is set according to the ground physical simulation test results of the satellite attitude control system, that is, the theoretical loop delay under the rigid body without interference, which is set to 0.020 seconds.

[0082] The steps to obtain the parameters are as follows: read the full width at half maximum (FWHM) value of the cross-correlation peak point calculated in the previous steps. This value reflects the tailing effect of the signal caused by the flexible coupling oscillation. For example, the current measurement value is 0.080 seconds.

[0083] The steps for obtaining the parameter are as follows: the parameter is a preset waveform dispersion penalty factor, which is set using the critical fault balance method, and the maximum allowable delay error threshold value of the system is set. (Right now ), and the maximum permissible waveform diffusion width threshold. This ensures that the contributions of the two terms in the formula are equal under the two critical states.

[0084] Calculate the critical value of the delay term: ,

[0085] Calculate the critical baseline value of the diffusion term: ,

[0086] Calculate the balance factor: ,

[0087] Rounding settings .

[0088] Calculations based on parameters:

[0089] First, calculate the absolute value of the squared difference of the time delay terms:

[0090] ;

[0091] Calculate the numerator (width cube) of the waveform diffusion term:

[0092] ;

[0093] Calculate the denominator (sum of time delays) of the waveform dispersion term:

[0094] ;

[0095] Calculate the weighted diffusion term value:

[0096] ;

[0097] Calculate the total lag energy within the square root:

[0098] ;

[0099] Calculate the final lagged fraction (square root):

[0100] ;

[0101] The results show that the overall response lag of the current control loop is 0.1318 seconds, which is greater than the nominal delay of 0.020 seconds. Moreover, the contribution of the waveform dispersion term (0.015754) is much greater than that of the pure delay term (0.001625), indicating that the decrease in system synchronicity at this time is mainly caused by signal distortion and trailing caused by flexible vibration, rather than simple transmission blockage.

[0102] The steps for obtaining the feature vector of on-orbit anomaly state are as follows:

[0103] The control loop response synchronicity lag score is compared with a preset processor bus blocking threshold. If the control loop response synchronicity lag score is greater than the processor bus blocking threshold, it is determined to be a signal blocking delay state. If the control loop response synchronicity lag score is less than or equal to the processor bus blocking threshold, it is determined to be a signal transmission smooth state, and the signal transmission delay state of the control loop is generated. At the same time, the structural coupling instability early warning time window is retrieved, and all vibration complexity time-varying entropy values ​​within the structural coupling instability early warning time window are extracted. A linear regression algorithm is used to fit the trend line of the vibration complexity time-varying entropy value, and the slope value of the trend line is calculated to generate the entropy value decrease slope feature.

[0104] Based on the signal transmission delay state and entropy decrease slope characteristics of the control loop, a multi-dimensional feature mapping coordinate system is established. The signal transmission delay state of the control loop is quantized into coordinate values ​​of the first dimension, and the entropy decrease slope characteristics are quantized into coordinate values ​​of the second dimension. The coordinate values ​​of the first dimension and the coordinate values ​​of the second dimension are combined in sequence to generate an on-orbit anomaly state feature vector.

[0105] Specifically, the control loop response synchronization lag score is compared with a preset processor bus congestion threshold. First, the control loop response synchronization lag score calculated in the previous steps is read. This score incorporates physical delays and waveform distortions during satellite signal transmission. Then, the stored processor bus congestion threshold is retrieved. This threshold is set based on the real-time scheduling cycle of the satellite's onboard computer main bus. For example, when the bus clock frequency is 10MHz and the number of cycles for a single data frame transmission is 100, the basic transmission time is 0.01 milliseconds. Considering the waiting queue length under multi-task concurrency, the processor bus congestion threshold is set to 5 times the basic transmission time, i.e., 0.050 seconds. If the control loop response synchronization lag score is greater than 0.050 seconds, bus congestion is determined, and the signal congestion delay state is marked. Conversely, if the score is less than or equal to 0.050 seconds, the data exchange logic is determined to be normal, and the signal transmission is marked as smooth. The signal transmission delay state of the control loop is generated, and the structural coupling instability warning time window is retrieved from memory simultaneously, locating the start time covered by this window. With the deadline The system extracts all time-varying entropy values ​​of vibration complexity within the corresponding interval from the time-varying entropy value sequence of vibration complexity. Since the entropy value reflects the evolution of the orderliness of the signal, the least squares method is used to fit the extracted entropy value sequence to construct a linear regression equation. The slope value of the trend line is calculated by minimizing the sum of squared residuals, and the entropy value decrease slope feature is generated.

[0106] Based on the signal transmission delay state and entropy decrease slope characteristics of the control loop, a multi-dimensional feature mapping coordinate system is established, constructing a two-dimensional Cartesian coordinate plane. The first dimension is set as the time delay dimension, and the second dimension is set as the stability degradation dimension. The signal transmission delay state of the control loop is numerically quantized, and the quantization rule is set as follows: if the state is a signal blocking delay state, the coordinate is assigned as 1; if the state is a signal transmission smooth state, the coordinate is assigned as 0. This transforms the qualitative transmission description into coordinate values ​​of the first dimension. Next, the entropy decrease slope characteristics are processed, and the absolute value is taken and normalized. The maximum instability slope is set as a reference of -0.5. The ratio of the current slope value to this value is used as the coordinate value of the second dimension. For example, when the slope is -0.25, the quantized second dimension value is 0.5. The coordinate values ​​of the first dimension and the second dimension are combined in sequence, and a feature vector containing two elements is constructed by array concatenation, generating the on-orbit anomaly state feature vector.

[0107] The steps for obtaining the satellite's on-orbit anomaly detection and classification results are as follows:

[0108] The on-orbit anomaly feature vector is projected into the pre-constructed fault mode topology space. The Euclidean distance between the on-orbit anomaly feature vector and each standard fault node in the fault mode topology space is calculated. The standard fault node with the smallest Euclidean distance is searched as the best matching node. The fault type description information and fault level label bound to the best matching node are read to generate the satellite on-orbit anomaly detection and classification results.

[0109] Specifically, the on-orbit anomaly state feature vector is projected into a pre-constructed fault mode topology space. This space consists of multiple standard fault nodes, each representing a typical on-orbit fault mode. For example, the coordinates of node one are... Corresponding to the micro-oscillation driven by the solar panel, the coordinates of node two are: For attitude control anomalies caused by scheduling conflicts with the onboard computer, calculate the Euclidean distance between the on-orbit anomaly state feature vector and each standard fault node in the fault mode topology space. The calculation formula is as follows: in, Represents the characteristic vector of on-orbit anomaly state and the first Euclidean distance between standard fault nodes The first dimension coordinates of the on-orbit anomaly state feature vector are given. The coordinates of the second dimension of the on-orbit anomaly state feature vector are given. For the first The x-coordinate of each standard fault node, For the first The system uses the ordinate of a standard fault node, iterates through all calculated distance values, executes a minimum distance search strategy, identifies the node with the smallest value as the best matching node, accesses the pre-associated database fields of that node, retrieves the corresponding text description (e.g., excessive control algorithm loop delay) and fault level label (e.g., Level 2 warning), and generates a satellite on-orbit anomaly detection classification result. This classification covers multi-dimensional scenarios ranging from minor oscillations to complete loss of control. Typical results and warning levels are as follows: Normal operation Level 0 warning: Smooth command response, stable system operation. Solar panel-driven micro-oscillation Level 1 warning: Minor local anomaly, indicating initial fatigue of flexible attachments. Onboard computer scheduling conflict Level 2 warning: Excessive control algorithm loop delay, posing a risk of attitude control loss of synchronization. Flexible attachment fracture or deformation Level 3 warning: Rapid drop in flutter entropy, structurally facing severe coupling instability. Satellite attitude roll loss of control Level 4 warning: Complete disconnect between command response and feedback, facing a crisis of system-wide paralysis. Through this precise classification, the system can effectively cover various on-orbit risks, providing intuitive decision support for subsequent fault isolation and attitude recovery strategies.

Claims

1. A satellite on-orbit anomaly detection system based on AI dynamic behavior analysis, characterized in that, The system includes: The attitude signal mode decomposition module is used to acquire the attitude angular velocity feedback signal and flywheel control command sequence of the satellite attitude control system, extract the high-frequency flutter component of the flexible accessory, and filter the high-frequency flutter component of the flexible accessory according to the signal energy ratio to generate a key flexible mode eigenfunction sequence. The vibration entropy trend analysis module is used to divide the key flexible mode eigenfunction sequence into continuously sliding time slices, calculate the permutation entropy value of the signal amplitude sequence in each time slice, obtain the vibration complexity time-varying entropy value sequence, identify the time interval in which the permutation entropy value shows a continuous downward trend based on the vibration complexity time-varying entropy value sequence, and generate a structural coupling instability early warning time window. The command response cross-correlation module is used to extract the flywheel control command sequence and attitude angular velocity feedback signal segment corresponding to the structural coupling instability early warning time window, calculate the cross-correlation function value, search for the time delay corresponding to the maximum point of cross-correlation function amplitude, generate command response peak offset, and calculate the control loop response synchronization lag fraction based on the command response peak offset. The fault vector mapping determination module is used to compare the control loop response lag score with a preset processor bus blocking threshold to construct an on-orbit abnormal state feature vector, and to generate satellite on-orbit abnormality detection classification results based on the on-orbit abnormal state feature vector.

2. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 1, characterized in that, The steps for obtaining the key flexible mode eigenfunction sequence are as follows: Based on the attitude angular velocity feedback signal of the satellite attitude control system, the attitude angular velocity feedback signal is frequency domain partitioned according to the preset bandwidth limit condition. The center frequency values ​​of different frequency bands in the attitude angular velocity feedback signal are iteratively updated until the center frequency values ​​converge to a stable state. The attitude angular velocity feedback signal is decomposed into multiple intrinsic mode function components of different frequency bands. The signal part with a frequency value greater than the preset rigid body cutoff frequency value in the intrinsic mode function components is extracted. The low frequency components of rigid body motion with a frequency value less than or equal to the rigid body cutoff frequency value are removed to generate high frequency flutter components of flexible attachments. The instantaneous amplitude square term of each modal component in the high-frequency flutter component of the flexible attachment is calculated. The instantaneous amplitude square term is integrated and accumulated in the time domain to obtain the absolute energy value of each modal component. The absolute energy value is divided by the sum of the absolute energy values ​​of all modal components to calculate the weight ratio of each modal component in the overall vibration and generate the vibration mode energy ratio value. The vibration mode energy ratio is compared with a preset primary mode screening threshold. Modal components whose vibration mode energy ratio is greater than the primary mode screening threshold are identified. These modal components are marked as primary vibration modes. Secondary modes whose vibration mode energy ratio is less than or equal to the primary mode screening threshold are removed. The time-domain waveform data corresponding to the retained primary vibration modes are recombined to generate a sequence of key flexible mode eigenfunctions.

3. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 1, characterized in that, The steps for obtaining the time-varying entropy value sequence of vibration complexity are as follows: Set a fixed length of time window width and movement step size, continuously extract the sequence of eigenfunctions of the key flexible modes along the time axis, obtain multiple local time slices, reconstruct the phase space of the signal amplitude data in each time slice and count the occurrence probability of different symbol modes, calculate the quantitative index characterizing the disorder of the system based on Shannon information theory, and generate a time-varying entropy value sequence of vibration complexity. Calculate the first-order difference value based on the time-varying entropy value sequence of the vibration complexity.

4. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 3, characterized in that, The steps for obtaining the structural coupling instability early warning time window are as follows: The first-order difference value is compared with a preset negative change threshold to filter out anomalies where the first-order difference value is less than the negative change threshold. The start and end positions of the continuous occurrence of anomalies on the time axis are identified, the duration of the continuous anomalies is calculated and it is determined whether it exceeds the minimum instability confirmation period. The time period covered by the start and end positions that meet the conditions is marked as the risk interval, and a structural coupling instability early warning time window is generated.

5. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 1, characterized in that, The steps for obtaining the peak offset of the instruction response are as follows: Using the structural coupling instability early warning time window as the interception benchmark, discrete data segments corresponding to the time period are extracted from the flywheel control command sequence and the attitude angular velocity feedback signal, respectively. Cross-correlation convolution operation is performed on the two discrete data segments to generate a cross-correlation function numerical sequence. Based on the numerical sequence of the cross-correlation function, the point with the largest amplitude in the sequence is located by the peak search algorithm. The distance between the time coordinate corresponding to the maximum value point and the zero time is determined, and the peak offset of the command response is generated. At the same time, the span of the peak waveform on the time axis is calculated with half of the maximum amplitude as the boundary, and the full width at half maximum (FWHM) value of the cross-correlation peak point is generated.

6. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 5, characterized in that, The steps for obtaining the synchronization lag fraction of the control loop response are as follows: The control loop response lag score is calculated based on the peak offset of the command response and the full width at half maximum (FWHM) value of the cross-correlation peak point.

7. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 1, characterized in that, The steps for obtaining the on-orbit anomaly state feature vector are as follows: The control loop response lag score is compared with a preset processor bus blocking threshold. If the control loop response lag score is greater than the processor bus blocking threshold, it is determined to be a signal blocking delay state. If the control loop response lag score is less than or equal to the processor bus blocking threshold, it is determined to be a signal transmission smooth state, generating a signal transmission delay state for the control loop. At the same time, the structural coupling instability early warning time window is retrieved, and all vibration complexity time-varying entropy values ​​within the structural coupling instability early warning time window are extracted. A linear regression algorithm is used to fit the trend line of the vibration complexity time-varying entropy value, and the slope value of the trend line is calculated to generate an entropy value decrease slope feature. Based on the signal transmission delay state of the control loop and the entropy decrease slope feature, a multi-dimensional feature mapping coordinate system is established. The signal transmission delay state of the control loop is quantized into coordinate values ​​of the first dimension, and the entropy decrease slope feature is quantized into coordinate values ​​of the second dimension. The coordinate values ​​of the first dimension and the coordinate values ​​of the second dimension are combined in sequence to generate an on-orbit anomaly state feature vector.

8. The satellite on-orbit anomaly detection system based on AI dynamic behavior analysis according to claim 1, characterized in that, The steps for obtaining the satellite's on-orbit anomaly detection and classification results are as follows: The on-orbit anomaly feature vector is projected into the pre-constructed fault mode topology space. The Euclidean distance between the on-orbit anomaly feature vector and each standard fault node in the fault mode topology space is calculated. The standard fault node with the smallest Euclidean distance is searched as the best matching node. The fault type description information and fault level label bound to the best matching node are read to generate the satellite on-orbit anomaly detection classification result.

Images