A method and system for image stabilization processing in shipborne mobile satellite communication

Abstract

This invention provides an image stabilization processing method and system for shipborne satellite communication in motion, relating to the fields of satellite communication and digital image processing technology. The method includes: performing time-series alignment processing based on sampling clock deviation on synchronously acquired three-axis angular velocity signals and original image frame sequences to obtain a time-synchronized initial rigid body motion sequence of the ship; extracting corresponding image data based on the initial rigid body motion sequence of the ship to obtain a time-matched original image dataset; performing multi-scale fluid-structure interaction (FSI) frequency domain separation operations based on the initial rigid body motion sequence of the ship to obtain a mixed frequency domain spectrum; and filtering out low-frequency rigid body swaying components based on the mixed frequency domain spectrum to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. This invention effectively overcomes the shortcomings of existing technologies, such as insufficient image stabilization accuracy and increased satellite beam pointing error due to neglecting local elastic deformation of the ship.

Description

Technical Field

[0001] This invention relates to the fields of satellite communication and digital image processing technology, and in particular to an image stabilization processing method and system for shipborne mobile satellite communication. Background Technology

[0002] With the development of the marine economy and the growth in demand for ocean shipping, shipborne satellite communication systems have become the core equipment for high-speed data transmission, video surveillance and emergency communication on ships in the open sea. During navigation, the hull is subjected to complex loads such as waves, sea winds and ocean currents, which cause violent swaying, rolling and vibration, resulting in significant changes in the attitude of the antenna and imaging equipment, affecting image quality and beam pointing accuracy.

[0003] Currently, rigid body motion compensation technology based on inertial measurement units is widely used. While it can achieve some image stabilization when sea conditions are good and elastic deformation is small, it has the following drawbacks: Based on the assumption of an ideal rigid body, it ignores the influence of local elastic deformation of the hull. Specifically, the image stabilization accuracy is insufficient, making it difficult to eliminate pixel-level aperiodic phase distortion caused by local elastic deformation. In medium to high sea states, the high-frequency, small-amplitude elastic vibrations generated by the hull cannot be effectively compensated by rigid body motion sensors and servo mechanisms, manifesting as residual jitter and motion blur in imaging, leading to loss of image details and blurred edges. Simultaneously, elastic deformation causes a small but high-frequency deflection of the antenna optical axis, exceeding the response bandwidth of traditional beam tracking systems. This results in the beam not being accurately aligned with the satellite, causing signal attenuation, increased bit error rate, and even communication interruption in severe sea conditions. Since the imaging equipment and antenna in shipborne mobile communication terminals typically use a coaxial rigid mounting structure, elastic deformation, while causing antenna optical axis deflection, also leads to a small high-frequency deflection of the camera optical axis, further exacerbating image distortion.

[0004] The aforementioned problems are particularly prominent on large container ships, ocean-going research vessels, and other vessels with long hulls and high flexibility. Summary of the Invention

[0005] This invention provides an image stabilization processing method and system for shipborne mobile satellite communication, which can simultaneously compensate for the overall rigid body motion and local elastic deformation of the ship, effectively eliminate pixel-level aperiodic phase distortion, reduce satellite beam pointing error, and improve the imaging quality and communication reliability of the shipborne mobile satellite communication system in complex marine environments.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: A first aspect includes an image stabilization processing method for shipborne mobile satellite communication, the method comprising: Step 1: Perform timing alignment processing based on sampling clock deviation on the synchronously acquired three-axis angular velocity signals and the original image frame sequence to obtain a timing-synchronized initial rigid body motion sequence of the hull; extract the corresponding image data according to the initial rigid body motion sequence of the hull to obtain a timing-matched original image dataset. Step 2: Perform multi-scale fluid-structure rigid-flexible frequency domain separation operation based on the initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum; filter out the low-frequency rigid body swaying component based on the mixed frequency domain spectrum to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. Step 3: Construct the optomechanical coupled spatial projection mapping relationship based on the high-frequency flexible deformation component sequence to obtain the deformation projection transformation matrix; Analyze the undersampling aliasing trajectory within the exposure period based on the deformation projection transformation matrix to obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. Step 4: Correct the initial rigid body motion sequence of the hull according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory; use the preset rigid body rotation inertia smoothing constraint to perform time-frequency domain nonlinear fusion on the first-order compensated motion trajectory to obtain a continuous and non-aliased high-precision steady-state motion vector field. Step 5: Perform inverse geometric transformation and motion blur suppression processing on the original image dataset of the time-matched sequence frame by frame according to the high-precision steady-state motion vector field to obtain a set of corrected frame images; perform closed-loop verification of pointing error according to the set of corrected frame images to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

[0007] Secondly, an image stabilization processing system for shipborne satellite communication in motion includes: The timing alignment module is used to perform timing alignment processing based on sampling clock deviation on the synchronously acquired three-axis angular velocity signals and the original image frame sequence to obtain the timing-synchronized initial rigid body motion sequence of the hull, and extract the corresponding image data to obtain the timing-matched original image dataset. The rigid-flexible frequency domain separation module is used to perform multi-scale fluid-structure rigid-flexible frequency domain separation operations on the initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum. Based on the mixed frequency domain spectrum, the low-frequency rigid body swaying component is filtered out to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. The distortion field analysis module is used to construct the optomechanical coupling spatial projection mapping relationship based on the high-frequency flexible deformation component sequence, obtain the deformation projection transformation matrix, analyze the undersampled aliasing trajectory within the exposure period, and obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. The motion trajectory fusion module is used to correct the initial rigid body motion sequence of the hull according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory. The first-order compensated motion trajectory is nonlinearly fused in the time and frequency domain using a preset rigid body rotation inertia smoothing constraint to obtain a continuous and non-aliased high-precision steady-state motion vector field. The image correction and closed-loop verification module is used to perform inverse geometric transformation and motion blur suppression processing on the original image dataset of the time-matched model frame by frame according to the high-precision steady-state motion vector field to obtain a set of corrected frame images, and then perform closed-loop verification of pointing error to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

[0008] The above-described solution of the present invention has at least the following beneficial effects: By employing multi-scale fluid-structure interaction (FSI) rigid-flexible frequency domain separation, low-frequency rigid body sway components are precisely filtered out from the mixed frequency domain spectrum, and high-frequency flexible deformation characterizing the local elastic deformation of the hull is extracted. Based on optomechanical coupling projection, a deformation mapping model is constructed to analyze the pixel-level displacement distortion field formed by aperiodic phase distortion within the exposure cycle. This overcomes the technical difficulties of traditional rigid body compensation methods that ignore hull elastic deformation and cannot eliminate high-frequency residual jitter and motion blur. Furthermore, through reverse phase compensation and rigid body rotational inertia smoothing constraints, the initial compensation trajectory is nonlinearly fused in the time and frequency domains to generate a continuous, non-aliased, high-precision steady-state motion vector field, effectively suppressing image edge blurring and detail loss. Simultaneously, combined with pointing error closed-loop verification, satellite beam pointing error is effectively reduced, ultimately obtaining a high-quality stable video sequence with clear and stable images and reliable communication links. This effectively improves the performance of the on-the-move communication system for large flexible ships in complex sea conditions. Attached Figure Description

[0009] Figure 1 This is a flowchart illustrating an image stabilization processing method for shipborne satellite communication in motion, provided by an embodiment of the present invention. Figure 2 This is a schematic diagram of an image stabilization processing system for shipborne satellite communication in motion, provided by an embodiment of the present invention. Detailed Implementation

[0010] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0011] like Figure 1As shown, an embodiment of the present invention proposes an image stabilization processing method for shipborne satellite communication in motion, the method comprising the following steps: Step 1: Perform timing alignment processing based on sampling clock deviation on the synchronously acquired three-axis angular velocity signals and the original image frame sequence to obtain a timing-synchronized initial rigid body motion sequence of the hull; extract the corresponding image data according to the initial rigid body motion sequence of the hull to obtain a timing-matched original image dataset. Step 2: Perform multi-scale fluid-structure rigid-flexible frequency domain separation operation based on the initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum; filter out the low-frequency rigid body swaying component based on the mixed frequency domain spectrum to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. Step 3: Construct the optomechanical coupled spatial projection mapping relationship based on the high-frequency flexible deformation component sequence to obtain the deformation projection transformation matrix; Analyze the undersampling aliasing trajectory within the exposure period based on the deformation projection transformation matrix to obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. Step 4: Correct the initial rigid body motion sequence of the hull according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory; use the preset rigid body rotation inertia smoothing constraint to perform time-frequency domain nonlinear fusion on the first-order compensated motion trajectory to obtain a continuous and non-aliased high-precision steady-state motion vector field. Step 5: Perform inverse geometric transformation and motion blur suppression processing on the original image dataset of the time-matched sequence frame by frame according to the high-precision steady-state motion vector field to obtain a set of corrected frame images; perform closed-loop verification of pointing error according to the set of corrected frame images to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

[0012] In this embodiment of the invention, precise timing alignment of the three-axis angular velocity signal with the original image frame ensures strict synchronization between motion data and image content. Based on this, a multi-scale fluid-structure interaction (FSI) technique is used to effectively extract the high-frequency flexible deformation component generated by local elastic deformation from the ship's hybrid motion, avoiding interference from low-frequency swaying. Furthermore, an optomechanical coupled projection model is constructed to convert high-frequency deformation into a pixel-level displacement distortion field, accurately characterizing the aperiodic phase distortion within the exposure cycle. The distortion field is used to correct the initial motion trajectory, and rigid body rotational inertia constraints are introduced for time-frequency domain nonlinear fusion, generating a continuous, non-aliased, high-precision steady-state motion vector field. Finally, inverse geometric transformation and adaptive deblurring eliminate the imaging effects of residual jitter, combined with closed-loop verification of pointing error, achieving clear and stable video output that simultaneously meets the requirements of motion blur radius and satellite beam pointing accuracy, effectively improving the communication and imaging reliability of shipborne mobile communications in complex sea conditions.

[0013] In a preferred embodiment of the present invention, step 1 above may include: Step 1.1: Analyze the local time-frequency energy distribution characteristics of the triaxial angular velocity signal, and construct a dynamic convolution constraint matrix based on these characteristics. Specifically, this includes: performing time-frequency transformation processing on the synchronously acquired triaxial angular velocity signal, mapping it from the time domain to the time-frequency joint domain to obtain the time-frequency energy distribution spectrum of the signal; this transformation process uses a sliding window short-time Fourier transform, with a fixed-length analysis window sliding along the time axis point by point. At each sliding step, a Fourier transform is performed on the signal segment within the window to obtain the complex amplitude of the signal at each frequency component within that local time period; the modulus of the complex amplitude is taken and the square is calculated to obtain the energy density value of that time-frequency unit. This is achieved by traversing the entire signal duration... For all window positions, a two-dimensional time-frequency energy distribution map is ultimately generated, with time as the horizontal axis, frequency as the vertical axis, and energy density as the pixel brightness. Based on this, the time-frequency energy distribution map is divided into several frequency slices along the frequency axis, with each slice corresponding to a narrow frequency range. For each frequency slice, the energy value of the current time-frequency unit is calculated point by point along the time axis, and the energy value of the adjacent time-frequency units is calculated to obtain the energy gradient amplitude at each time-frequency position on the slice. The energy gradient amplitudes of all frequency slices at the same time position are summarized, and the maximum value of the energy gradient amplitude at the corresponding time position in all slices is taken as the comprehensive gradient value at that time, forming a comprehensive gradient curve that varies along the time axis.

[0014] The sections of the curve with amplitudes significantly higher than the average level represent non-stationary energy distribution regions, indicating rapid switching of the ship's motion state during these periods. Furthermore, a dynamic curvature constraint matrix is ​​constructed based on these local time-frequency energy distribution characteristics. This matrix is ​​a two-dimensional grid with the same size as the two-dimensional time-frequency energy distribution map. Each element in the grid corresponds to a time-frequency unit on the time-frequency map. The element values ​​are determined by: presetting a baseline constraint value and two adjustment coefficients. In a preferred embodiment, the baseline constraint value is based on the sampling period of the triaxial angular velocity signal and the image frame period. The ratio is determined and is usually set to 0.5 to balance the requirements of time alignment accuracy and path smoothness. The relaxation factor is between 1.5 and 3.0. When the corresponding time position is in a non-stationary change region, the reference constraint value is multiplied by the relaxation factor to obtain a higher element value, indicating that a larger amount of time remapping curvature is allowed at this position. The tightening factor is between 0.2 and 0.5. When the corresponding time position is in a stationary region, the reference constraint value is multiplied by the tightening factor to obtain a lower element value, indicating that the time mapping at this position must maintain a higher rigidity.

[0015] The resulting dynamic curl constraint matrix has each element representing the maximum local deformation cost allowed at its corresponding time-mapping node. This matrix differentiates the bending amount at each step of the path search process, preventing unnecessary fluctuations in the path where the signal energy is stable and allowing the path to be flexibly adjusted to match the image frame sequence where the signal energy changes abruptly. This fully constrains the morphological range of the time-mapping path between the triaxial angular velocity signal and the original image frame sequence.

[0016] Step 1.2: Based on the dynamic curl constraint matrix and the sampling clock deviation, path optimization calculation is performed to obtain an adaptive time remapping trajectory. Specifically, this includes: using the sampling clock deviation as the initial offset for time alignment. This sampling clock deviation is caused by the frequency difference of the hardware clock source between the inertial measurement unit and the camera, manifested as an unknown but constant time difference between the angular velocity sampling time and the image exposure time. This clock deviation value is set as a fixed offset parameter for the initial mapping node in the path search space. Combined with the dynamic curl constraint matrix, path optimization calculation is performed in a two-dimensional discrete path search space spanned by the angular velocity sampling point index and the image frame index. The goal of this calculation is to find a path that minimizes the total cost path of the accumulated difference between the three-axis angular velocity signal sequence and the original image frame sequence, thereby obtaining an adaptive time remapping trajectory. Specifically, the path optimization calculation is achieved by solving a problem minimizing a cumulative cost function, the expression of which is: In the formula, This indicates the first of the three-axis angular velocity signals. The sampling point is mapped to the first sampling point in the original image frame sequence. The total cost accumulated along the path at frame rate; This indicates the arrival at the current mapping point. Previously, the minimum cumulative cost among the possible preceding path points, these three preceding points correspond to diagonal matching on the time axis, jump matching on the angular velocity axis and jump matching on the image axis, respectively; Indicates the first The angular velocity sampling point and the first The local matching distance of feature vectors between image frames, where the feature vectors are composed of the magnitude features of angular velocity sampling points and the motion blur features or gradient features of the image frames; local matching distance The construction of this feature space is crucial for connecting heterogeneous data domains. Since angular velocity signals and image frames belong to different physical domains, distance cannot be directly calculated. It is necessary to map both to an intermediate feature space representing the intensity of motion before measurement. The calculation formula is as follows: ;

[0017] In the formula, From the triaxial angular velocity signal, the first Motion intensity feature values ​​extracted at each sampling point For from the first Motion intensity feature values ​​extracted from frame images To perform absolute value operations, both are pre-normalized to an interval. ; The construction method is based on the first Taking a sampling point as the center, samples are taken forward and backward. Each sampling point forms a local window. The square root of the mean of the sum of the squares of the three-axis angular velocity components within the window is calculated as the local root-mean-square angular velocity of that sampling point, as shown in the following formula: ; in, The width of the local window is half the window width, and its value is the number of sampling points corresponding to the exposure time of a single frame of image that makes the window length cover the window width. This is the offset index of the sampling point within the window; , , These represent the components of the three-axis angular velocity vector along the roll, pitch, and yaw axes, respectively. The construction method is as follows: for the first The negative logarithmic transform of the Laplacian operator response mean of the entire frame's grayscale values ​​is calculated using the following formula: ; In the formula, This represents the total number of valid pixels in the image. For pixel index, For the first The first frame of the image grayscale value of each pixel. Represents the discrete Laplace operator. For taking values The smallest positive number is used to avoid taking the logarithm of 0; when the image is blurred due to motion, the edge sharpness decreases, the mean of the Laplacian response decreases, and after the negative logarithmic transformation... The increase is consistent with the trend of angular velocity increasing during intense motion. By recursively solving this function using dynamic programming, calculating the minimum cumulative cost point-by-point across the entire two-dimensional search space, and backtracking the sequence of nodes traversed by the optimal path, a path that minimizes the total cost can be determined. The minimum optimal path, and the mapping relationship corresponding to this optimal path, is the adaptive time remapping trajectory.

[0018] Step 1.3: Perform variable-step interpolation on the triaxial angular velocity signal according to the adaptive time remapping trajectory to obtain a resampled angular velocity sequence. Specifically, this includes: the adaptive time remapping trajectory defines a non-uniform mapping relationship between the original triaxial angular velocity signal sampling time and the resampled time that should be output; this mapping relationship is represented by a series of discrete mapping pairs, each mapping pair indicating an original sampling time and its corresponding target resampled time. Based on this mapping relationship, perform variable-step interpolation on the original triaxial angular velocity signal to generate a resampled angular velocity sequence that is precisely aligned with the image frame in time; the variable-step interpolation operation is implemented through a time-distance-based weighted average algorithm, the calculation formula of which is: ;

[0019] in, Indicates the time of target resampling The angular velocity value calculated at that location; This indicates that in the original triaxial angular velocity signal, the time position is... The first in the neighborhood angular velocity values ​​of a known sampling point; Indicates the first Known sampling points The contribution weight to the interpolation result, which is based on the known sampling point time. With the target time The time distance between them was calculated. This represents a summation operation performed on all known sampling points within the neighborhood involved in the calculation; in a preferred embodiment, the weights... The Gaussian kernel function is used for definition, and its calculation formula is as follows: In the formula, Represented by natural constant An exponential function with base 0; The bandwidth parameter of the Gaussian kernel controls the rate at which the weights decay with increasing time distance. The value is preset based on the temporal resolution and noise level of the angular velocity signal; the weight is calculated using this Gaussian kernel function, so that the closer the known sampling point is in time, the greater its contribution to the interpolation result, and the farther away the contribution is, the smaller the contribution, thus obtaining a smooth and locally adaptive interpolation effect; by traversing all target resampling times and performing this weighted average operation, the angular velocity value at each target time can be obtained, completing the non-uniform resampling of the signal on the time axis, and obtaining a resampled angular velocity sequence that is precisely synchronized with the original image frame sequence on the time reference.

[0020] Step 1.4: Extract the frame time boundary of the resampled angular velocity sequence and the exposure window of the original image frame sequence; perform phase alignment matching on the frame time boundary and the exposure window to obtain a time-synchronized initial rigid body motion sequence of the hull, specifically including: based on this, after obtaining the resampled angular velocity sequence, parse the timestamp information attached to each sampling point in the sequence. The timestamp accurately records the absolute time corresponding to each angular velocity value after resampling; calculate the nominal time interval between two consecutive images according to the frame rate of the original image frame sequence, and accurately define the time boundary of the start and end of the exposure of each image frame, i.e., the frame time boundary. The start time of the frame time boundary corresponds to the start time of the exposure of the corresponding frame, and the end time corresponds to the end time of the exposure of the corresponding frame. The difference between the two is the exposure duration; at the same time, extract the exposure window of each frame in the original image frame sequence. This window is jointly defined by the exposure start time and exposure duration recorded in the camera hardware trigger signal or image metadata.

[0021] Furthermore, a phase alignment matching operation is performed on each pair of frame time boundaries and exposure windows. Specifically, all angular velocity sampling points in the resampled angular velocity sequence whose timestamps fall completely within the current frame time boundary are formed into a data segment. The weighted time centroid of the angular velocity values ​​within this data segment is calculated, and this weighted time centroid is associated with the center time of the image exposure window. The formula for calculating the weighted time centroid of the angular velocity values ​​within this data segment is as follows: In the formula, This indicates the weighted time centroid of the calculated data segment; Indicates the first [number]th [item] in this data segment Timestamp of each sampling point; Indicates the first [number]th [item] in this data segment The triaxial angular velocity vector of each sampling point consists of three scalar components, namely ; This represents the three-axis angular velocity vector. The modulus is calculated as follows: ; This is used to measure the instantaneous intensity of the ship's motion at the sampling moment; this weighted temporal centroid is weighted by the angular velocity modulus, so that the more intense the motion, the greater the contribution to the centroid position, thus more accurately reflecting the temporal concentration of the ship's motion energy during the exposure of that frame; by... By associating and matching with the center moment of the image exposure window, it is ensured that the hull motion state corresponding to the exposure period of the image frame can be completely and unbiasedly represented by the angular velocity data segment that is precisely aligned with it. After traversing all image frames and completing the above phase alignment matching, a set of time-synchronized initial rigid body motion sequences of the hull is finally constructed, in which each image frame has a unique corresponding rigid body motion data segment.

[0022] Step 1.5: Extract the corresponding original image data based on the time-synchronized initial rigid body motion sequence of the hull, and perform frame-level data serialization and encapsulation to obtain a time-matched original image dataset. Specifically, this includes: using the time-synchronized initial rigid body motion sequence of the hull as an index, extracting the original image data frame that perfectly corresponds in timestamp from the original video stream one by one. The corresponding judgment criterion is: the difference between the timestamp of the original image frame and the timestamp of the data segment in the initial rigid body motion sequence of the hull is less than a preset synchronization tolerance threshold. In a preferred embodiment, the synchronization tolerance threshold is set according to the synchronization hardware characteristics of the shipborne inertial measurement unit and the image acquisition device, and its value range is usually set to the exposure time of the image frame. to For a typical shipborne imaging system, if the exposure time of a single frame is 10ms, the synchronization tolerance threshold is 1 to 3ms; if the exposure time is 30ms, the synchronization tolerance threshold is 3 to 10ms. The reason for this is that if the threshold is too high, it will introduce motion information that is not in the current frame, resulting in inter-frame crosstalk in motion compensation; if the threshold is too low, it may discard too many valid frames due to hardware clock jitter, resulting in a decrease in data utilization.

[0023] Each extracted original image frame is bound to its synchronized rigid body motion data segment to form a data pair containing image pixel data and corresponding three-axis angular velocity data. All data pairs are encapsulated into frame-level data serialization according to a strict and unified time order, so that each encapsulated data unit constitutes a logically indivisible processing unit. All such data units are arranged in chronological order to form a time-matched original image dataset. The correspondence between image frames and motion data in this dataset has undergone strict clock deviation correction and phase alignment, providing a precise and synchronized data foundation for subsequent rigid-flexible frequency domain separation and accurate compensation.

[0024] In a preferred embodiment of the present invention, step 2 above may include: Step 2.1 involves performing a multi-scale fluid-structure-rigid-flexible frequency domain expansion operation on the time-synchronized initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum. Specifically, this includes: using the time-synchronized initial rigid body motion sequence of the hull as input, where each data segment contains discrete sampled values ​​of the three-axis angular velocity vector, precisely aligned with a frame of image, varying over time; the sampling rate is typically 100Hz to 200Hz; and performing a multi-scale fluid-structure-rigid-flexible frequency domain separation operation on the motion sequence along the time axis. This multi-scale fluid-structure-rigid-flexible frequency domain separation operation refers to using multiple sets of bandpass filter banks with different time and frequency resolutions to expand the initial rigid body motion sequence of the hull at multiple time and frequency scales to separate the low-frequency rigid body swaying component excited by wave fluid loads from the high-frequency flexible deformation component of the hull structure's elastic response. The core idea is to construct a set of complex analytic bandpass filters with different quality factors, where the center frequency of each filter is determined according to... The octave intervals are evenly distributed on the logarithmic frequency axis, covering from 0.03Hz to the sampling rate. The complete analysis frequency band; at each scale level, the motion sequence is segmented using an analysis time window length corresponding to the center frequency of that level. The window length is set to 8 to 16 times the period corresponding to the center frequency of that level to ensure sufficient frequency resolution for low-frequency bands and sufficient time resolution for high-frequency bands.

[0025] The process involves sliding the window along the time axis in 50% increments, performing a Fast Fourier Transform on the angular velocity signal segment within each window to obtain the complex spectrum within each local time region at that scale level. For the complex spectra at the same time window position across all scale levels, they are stacked and pieced together along the frequency axis in ascending order of frequency to form a local two-dimensional time-spectrum slice for that moment. These local two-dimensional time-spectrum slices are then arranged sequentially along the time axis to ultimately form a complete hybrid frequency domain spectrum jointly unfolded along the time, frequency, and scale axes. When the length of the lowest frequency band analysis window exceeds the total length of available data... The window data is supplemented by endpoint mirror extension or zero-phase filling to ensure the integrity of multi-scale frequency domain expansion calculation. The hybrid frequency domain spectrum completely preserves all motion information from low frequency to high frequency in the initial rigid body motion sequence of the hull. In the low frequency region of 0.03Hz to 0.3Hz, the overall rigid body swaying motion component of the hull directly excited by wave fluid load is mainly concentrated. In the mid-high frequency region of 0.3Hz to 10Hz, the elastic vibration response of the local structure of the hull and the excitation of the deck machinery are mixed. In the high frequency region above 10Hz, the inherent measurement noise component of the inertial measurement unit is mainly included.

[0026] Step 2.2: Analyze the spectral energy distribution gradient of the mixed frequency domain spectrum, and extract the motion component within the wave fluid load excitation frequency band based on the spectral energy distribution gradient to obtain the low-frequency interference characteristic spectrum. Specifically, this includes: after obtaining the mixed frequency domain spectrum, scanning the energy density distribution of the spectrum layer by layer along the frequency axis. Specifically, the mixed frequency domain spectrum is represented as a two-dimensional matrix, where the row index corresponds to discrete frequency points and the column index corresponds to discrete time frames; let the frequency index be... , ,in Corresponding to the lowest frequency, Corresponding to the highest frequency; time index is , For each fixed time series Starting from the lowest frequency, calculate the index position of each frequency sequentially. Energy density value It is immediately adjacent to the previous frequency index position Energy density value The first-order forward difference, whose value constitutes the difference at that moment. The spectral energy distribution gradient sequence; in a preferred embodiment, the calculation formula for the first-order forward difference is: ; In the formula For time frames Frequency Index The energy gradient along the frequency axis; The frequency index in the mixed frequency domain spectrum is Time index is The energy density value of the time-frequency unit; The frequency index of the adjacent low-frequency side is: The energy density value of the time-frequency unit; by analyzing all Perform the above difference operation on each time series to obtain a complete time-varying gradient matrix. Then, calculate the arithmetic mean of this gradient matrix along the time axis, that is, for each frequency index... The mean of the gradient values ​​across all time frames is calculated to obtain a one-dimensional global gradient curve that varies only with frequency. Its form is as follows: This global gradient curve The average rate of change of spectral energy density along the frequency direction was characterized throughout the entire observation period; based on this, the global gradient curve was calculated. Overall average value across the entire analysis band Traverse from the lowest frequency end towards the higher frequency end. The curve searches for the first curve that satisfies the frequency index within the 0.3Hz to 1Hz frequency band. The local peak position, which marks the abrupt boundary of the transition from rigid body swaying motion energy to elastic vibration energy; the frequency value corresponding to this peak position. The cutoff frequency was determined, and the frequency range below the cutoff frequency up to 0Hz was identified as the wave fluid load excitation frequency band. This excitation frequency band precisely corresponds to the low-frequency rigid body motion response generated by the coupling effect between waves and the ship hull in the actual marine environment. Specifically, under typical medium-high sea states of level four to six, the wave energy is mainly concentrated in the frequency band with a period of 4 to 12 seconds, corresponding to a wave frequency range of 0.08Hz to 0.25Hz. The natural frequencies of the ship's roll, pitch, and heave motions in the waves usually fall between 0.05Hz and 0.3Hz, and the wave encounter frequency will undergo a Doppler shift due to the influence of ship speed and wave angle, so that the wave fluid load excitation energy actually acting on the ship hull is concentrated in a wide frequency band of 0.03Hz to 0.3Hz.

[0027] Therefore, the identified cutoff frequency typically falls around 0.3Hz, which can fully cover the main energy frequency band of the wave; after determining the excitation frequency band, all frequency components in the mixed frequency domain spectrum are set to 0Hz to... The time-frequency components within the interval are extracted separately and recombined according to their original time and frequency positions to form a low-frequency interference characteristic spectrum that is strictly isomorphic to the mixed frequency domain spectrum in the time and frequency dimensions. This characteristic spectrum accurately depicts the intensity fluctuations on the time axis and the energy broadening characteristics on the frequency axis of the three rigid body motion modes of overall roll, pitch and yaw generated by the ship under the action of waves.

[0028] Step 2.3: Construct a frequency domain dynamic isolation mask based on the low-frequency interference feature spectrum. Perform low-frequency energy removal operation on the mixed frequency domain spectrum using the frequency domain dynamic isolation mask to obtain the high-frequency structural response frequency domain components. Specifically, this includes: constructing a frequency domain dynamic isolation mask with the same size as the low-frequency interference feature spectrum extracted in Step 2.2. This frequency domain dynamic isolation mask is a two-dimensional matrix consisting only of values ​​0 and 1, with rows corresponding to frequency indices and columns corresponding to time indices. The construction rule is as follows: traverse each time-frequency unit in the low-frequency interference feature spectrum. If the energy density value of the unit is greater than the preset energy detection threshold, then the element at the same row and column position in the mask matrix is ​​assigned a value of 0, indicating that the time-frequency position needs to be isolated and filtered out; if the energy density value of the unit is less than or equal to the energy detection threshold, then the corresponding element in the mask matrix is ​​assigned a value of 1, indicating that the time-frequency position is allowed to pass.

[0029] The energy detection threshold is set to 20% of the average energy density of all time-frequency units in the low-frequency interference characteristic spectrum to ensure that the core energy region of wave load excitation is completely covered. For situations where sea conditions change over time, causing wave load frequency drift, the 0-value region of the mask dynamically adjusts according to the time evolution of energy distribution in the low-frequency interference characteristic spectrum, allowing the mask to have an adaptive isolation range for different time periods. After construction, the frequency domain dynamic isolation mask is multiplied element-wise with the original mixed frequency domain spectrum: the low-frequency time-frequency components in the mixed frequency domain spectrum located within the wave fluid load excitation frequency band are multiplied by the 0 of the mask and set to 0, while the mid-to-high frequency time-frequency components located above the excitation cutoff frequency are fully preserved after multiplying by the 1 of the mask. Through this low-frequency energy removal operation, the low-frequency rigid body rocking components in the mixed frequency domain spectrum are cleanly filtered out, leaving only the high-frequency structural response frequency domain component containing only the high-frequency dynamic response components of the hull structure to wave load. In this component, the elastic vibration signal in the 0.3Hz to 10Hz frequency band and the measurement noise above 10Hz are clearly distinguishable in the spectrum.

[0030] Step 2.4 involves performing inverse time-domain phase reconstruction based on the frequency domain components of the high-frequency structural response to obtain a sequence of high-frequency flexible deformation components characterizing local elastic deformation. Specifically, this includes: using the obtained frequency domain components of the high-frequency structural response as the reconstruction object, performing an inverse transformation reconstruction from the frequency domain to the time domain; for each local time window in the frequency domain components of the high-frequency structural response, extracting the amplitude spectrum and phase spectrum information of its complex spectrum; the amplitude spectrum directly characterizes the intensity distribution of high-frequency structural vibration at each frequency component within that local time period, while the phase spectrum records the relative time delay relationship between each frequency component in radians; in the high-frequency band above 10Hz, an amplitude threshold is set to determine frequency domain components below this threshold as measurement noise and set them to 0. This amplitude threshold is set to the mean of the high-frequency amplitude spectrum plus three times the standard deviation of the amplitude spectrum, in order to suppress noise to the maximum extent while preserving the elastic vibration signal.

[0031] After processing, the inverse fast Fourier transform algorithm, corresponding to the forward fast Fourier transform in step 2.1, is used to map the corrected complex spectrum back to the time domain window by window, obtaining high-frequency vibration time-domain waveform segments in each local time region. Since the window sliding step size in step 2.1 is 50%, there is an overlapping region between adjacent windows corresponding to different time periods. Within the overlapping region, a cosine square function is used as a transition weight to gradually cross-merge the overlapping portion of the tail of the previous window segment and the head of the next window segment. In a preferred embodiment, the specific expression of the cosine square function is: In the formula, Indicates the first [unit] within the overlapping region The weighting coefficients applied to the previous window segment at each sampling point; This indicates the position index of the sampling point within the overlapping region. ,common There are points, among which The starting point of the corresponding overlapping region, The end point of the corresponding overlapping region; This represents the total number of sampling points within the overlapping area; the corresponding weight coefficient applied to the next window segment is... This ensures that the sum of the two weights is always 1, and the fused waveform maintains dual continuity of amplitude and phase within the overlapping region, without any jumps or discontinuities at the splicing point. After splicing, a complete and phase-continuous high-frequency flexible deformation component sequence is obtained. This sequence is strictly aligned with the original initial rigid body motion sequence of the hull on the time axis, and the sampling rate remains consistent. However, the motion information it contains has been transformed from broadband rigid-flexible hybrid motion into a set of small-amplitude, high-frequency vibrations that only characterize the local elastic deformation of the hull. The value at each sampling moment in the sequence is a three-dimensional vector, and the three components correspond to the local elastic angular deformation rates along the roll axis, pitch axis, and bow axis at the antenna and imaging equipment mounting base, respectively. This provides accurate elastic deformation excitation source data for constructing the optomechanical coupling projection mapping relationship in subsequent steps.

[0032] In a preferred embodiment of the present invention, step 3 above may include: Step 3.1: Analyze the spatial deformation vector of the high-frequency flexible deformation component sequence. Based on the spatial deformation vector, establish a camera optical axis deflection mapping model using the optomechanical coupling spatial projection mapping relationship. The optomechanical coupling spatial projection mapping relationship refers to the rigid connection and transmission relationship between the optical imaging system and the antenna mechanical support structure in the shipborne terminal, coupling and mapping the angular displacement of the mounting base caused by the local elastic deformation of the hull into the deflection attitude of the camera optical axis. Then, through camera projection geometry, the optical axis deflection is mapped into a spatial coordinate transformation relationship of pixel displacement on the image plane. Spatial coordinate mapping calculations are performed using the camera optical axis deflection mapping model to obtain the terminal's local deformation spatial coordinate set, specifically including: The high-frequency flexible deformation component sequence is used as input. Each sampling moment in the sequence contains a three-dimensional vector, whose three components correspond to the local elastic angular deformation rates along the roll axis, pitch axis, and yaw axis at the antenna and imaging device mounting base, respectively. For each sampling moment in the sequence, an integral time microelement is traced back from that moment. The length of this integral time microelement is set to be consistent with the sampling period of the high-frequency flexible deformation component sequence, typically up to 10 ms. Within this microelement period, the angular deformation rates along the three axes are integrated over time to obtain the cumulative angular deformation amounts along the three axes from the start time of the microelement to the current sampling moment. In a preferred embodiment, the formula for calculating the cumulative angular deformation amount is: ; , Indicates at the sampling time Around The cumulative angular deformation of the shaft, in radians; The first in the high-frequency flexible deformation component sequence Each sampling time; It is the length of the integral time element, and its value is equal to the sampling period of the high-frequency flexible deformation component sequence; For variables in continuous time Around Instantaneous angular deformation rate of the shaft, Pick The corresponding component of the roll axis. Pick The corresponding component of the pitch axis. Pick The corresponding component of the bow rocker axis; This represents the definite integral operator, with the lower limit of integration being... The maximum number of points is In practical discrete signal processing, this integration operation is numerically approximated by applying the trapezoidal rule to the angular deformation rate values ​​of two sampling points within the integration interval. Through the above integration calculation, the cumulative angular deformation along the three axes together constitutes the spatial deformation vector at that sampling moment. Arrange the spatial deformation vectors of all sampling moments in the entire high-frequency flexible deformation component sequence in chronological order to form a complete spatial deformation vector sequence. This sequence accurately depicts the evolution of the elastic angular deformation at the mounting base on the time axis.

[0033] Based on this, a camera optical axis deflection mapping model is established using the aforementioned optomechanical coupling spatial projection mapping relationship. Specifically, it is described that there is a definite motion transmission relationship between the optical imaging system and the mechanical support structure in the shipborne mobile communication terminal. The micro-geometric displacement caused by the elastic deformation of the mechanical structure will be coupled to the optical system through the rigid body kinematic chain, thereby changing the orientation of the camera optical axis in the object space, and ultimately reflecting as a systematic offset of the image points on the image plane. This principle unifies the mechanical response of the mechanical structure and the optical projection geometry within the same mathematical framework, establishing an end-to-end mapping channel from structural deformation input to image distortion output. The model constructed based on this principle is based on the following physical facts: the antenna parabolic reflector of the shipborne mobile communication terminal and the optical imaging equipment are coaxially fixedly connected through a rigid support truss structure made of high-strength alloy steel, and the relative geometric relationship between the two remains constant after installation and calibration.

[0034] When the elastic deformation of a local area of ​​the ship's deck is transmitted to the mounting base of the supporting truss, the tiny angular displacement generated at the base is transmitted almost without loss to the camera mounting flange in the form of stress waves through the various members of the truss, causing a corresponding tiny deflection of the camera's optical axis in three-dimensional space. This deflection process follows the rigid body motion transmission law in structural mechanics, that is, within the order of 0.1 degrees of elastic deformation angular displacement, the transmission relationship can be regarded as linear. The model establishes a three-dimensional rectangular coordinate system that moves with the ship with the front node of the camera optical system as the origin. The first axis of this coordinate system points towards the bow along the horizontal plane, the second axis points towards the starboard side along the horizontal plane, and the third axis points towards the zenith along the vertical horizontal plane. The three axes form a right-handed system. In the reference state without any elastic deformation, the camera's optical axis looks forward along the first axis, the row direction of the image plane is parallel to the second axis, and the column direction is parallel to the third axis.

[0035] The model construction and parameter calibration process is as follows: During the shipyard dock mooring test, the ship is in a state of zero speed in still water, and the elastic deformation of the hull deck under known load conditions is negligible; three mutually orthogonal high-precision laser vibrometers are temporarily installed at the support truss mounting base, and a planar calibration target is placed directly in front of the camera's field of view. The target is printed with an array of equally spaced circular markers with a spacing of 50 mm, and the target plane is perpendicular to the camera's optical axis at a distance equal to the nominal working distance; a 5Hz frequency is applied near the support truss mounting base using a micro-torque electric vibrator. A 50Hz, amplitude-controllable sweeping sinusoidal excitation is used to simulate the local elastic deformation of the ship's hull. Simultaneously, the time-domain signals of the angular displacements along three axes at the mounting base, measured by three laser vibrometers, and the target image sequence captured by the camera are recorded. Cross-correlation analysis is performed on the collected mounting base angular displacement data and the pixel offset data of the marked points in the images. The least-squares fitting method is used to solve for the transfer coefficient matrix between the angular displacements along the three axes at the mounting base and the small deflection angles of the camera's optical axis around the second and third axes. In a preferred embodiment, the least-squares fitting expression is: ; Let be the transfer coefficient matrix to be solved, with a size of 2 rows and 3 columns; This represents the matrix that minimizes the summation expression that follows. Operation; It is the total number of sample groups collected for calibration data; ; It is the first The measured deflection angle vectors of the camera optical axis around the second and third axes, obtained from the target image in the sample group, are in the form of a column vector containing two elements. The first element is the pitch deflection angle around the second axis, and the second element is the azimuth deflection angle around the third axis. Indicates the first The three axial angular displacement vectors at the mounting base measured by the laser vibrometer in the sample group are in the form of a column vector containing three elements, which correspond to the angular displacements of the roll axis, pitch axis and yaw axis, respectively. Let represent the square of the Euclidean norm of a vector, which is the sum of the squares of all its elements. The optimal transfer coefficient matrix is ​​obtained by solving the least squares problem using singular value decomposition or normal equation methods. The element values ​​of this transfer coefficient matrix are jointly determined by the geometric dimensions of the supporting truss, the elastic modulus of the material, and the stiffness characteristics of the connection nodes. Through the above calibration process, a deterministic mapping relationship between the spatial deformation vector in the high-frequency flexible deformation component sequence and the instantaneous deflection attitude of the camera optical axis is established, and the camera optical axis deflection mapping model is completed.

[0036] During routine online operation, at each sampling moment in the spatial deformation vector sequence, the three components of the spatial deformation vector at that moment are substituted into the calibrated transfer coefficient matrix for multiplication to obtain the pitch deflection angle of the camera optical axis about the second axis and the azimuth deflection angle about the third axis at that moment. These two deflection angles, together with the initial optical axis orientation, determine the instantaneous deflection attitude of the camera optical axis at that moment. Based on this instantaneous deflection attitude, a three-dimensional rotation matrix describing the rotation of the camera coordinate system relative to the reference state is generated. The imaging terminal surface refers to the end face of the feed horn at the focal point of the satellite antenna parabolic reflector, on which 9 rows and 9 columns of 8 grids are evenly distributed with the feed center as the origin and a grid spacing of 10mm. One structural feature point; the distribution density of the 81 feature points (9×9 grid) is set according to the physical size of the feed horn end face (approximately 80mm×80mm) and the spatial Nyquist sampling criterion of the half wavelength of the highest frequency elastic deformation, ensuring that the deformation spatial coordinate set is non-aliased; the initial spatial coordinates of these feature points in the reference state are multiplied one by one by the three-dimensional rotation matrix to perform coordinate mapping transformation, and the new spatial coordinates of each feature point after optical axis deflection are calculated, thereby obtaining a set of terminal local deformation spatial coordinates containing 81 spatial coordinate points at the sampling time; the terminal local deformation spatial coordinates of all sampling times are aggregated in chronological order to finally form the terminal local deformation spatial coordinate set.

[0037] Step 3.2: Perform a two-dimensional image plane projection operation based on the terminal local deformation spatial coordinate set and the preset camera intrinsic parameters to obtain the deformation projection transformation matrix. Specifically, this includes: after obtaining the terminal local deformation spatial coordinate set, performing homogeneous coordinate transformation on the three-dimensional spatial coordinates of each feature point at each sampling time, representing it as a four-dimensional homogeneous coordinate vector in the camera coordinate system. The first three components of this vector are the first axis coordinate value, the second axis coordinate value, and the third axis coordinate value of the feature point in the camera coordinate system, respectively, and the fourth component is always 1. Introduce the preset camera intrinsic parameter matrix, which is a 3x3 upper triangular matrix. Its main diagonal elements are, in order, the ratio of the equivalent focal length of the camera lens divided by the physical size of the pixel unit in the row direction, the ratio of the equivalent focal length divided by the physical size of the pixel unit in the column direction, and the value 1. The first row and second column elements of the matrix are the non-perpendicularity coefficients between the row direction and the column direction, and the first row and third column elements and the second row and third column elements are the row coordinates and column coordinates of the principal point of the image plane in the pixel coordinate system, respectively.

[0038] Further, based on perspective projection geometry, the four-dimensional homogeneous coordinate vector of each feature point is multiplied by the preset camera intrinsic parameter matrix to obtain the two-dimensional homogeneous coordinates of the feature point in the pixel coordinate system. After normalization to make the third component equal to 1, the first two components are extracted to obtain the actual imaging position of the feature point on the image plane. This projection operation is performed on all 81 feature points within the same sampling time to obtain a set of image plane projection points corresponding to that sampling time. Each projection point in the image plane projection point set is compared one by one with the reference imaging position of the feature point under no deformation conditions, and the coordinate difference between the two along the row direction in the pixel coordinate system is calculated. The coordinate difference along the column direction yields a two-dimensional vector describing the pixel position offset caused by deformation at that sampling moment. The two-dimensional offset vectors of the 81 feature points are arranged according to their corresponding 9x9 grid spatial positions and assembled into a 9x9 deformation projection transformation matrix. The element in the r-th row and c-th column of the matrix is ​​a two-dimensional vector containing both row and column offsets, quantifying the pixel-level displacement mapping relationship generated by the local deformation at the r-th row and c-th column grid node on the imaging terminal surface on the image plane. The above operation is performed sequentially for all sampling moments to obtain a series of deformation projection transformation matrices that change continuously along the time axis.

[0039] Step 3.3: Extract the preset camera single-frame exposure time window, calculate the pixel cumulative offset path within the exposure time window based on the deformation projection transformation matrix, and obtain the undersampled aliasing trajectory within the exposure period. Specifically, this includes: after obtaining the complete time series of the deformation projection transformation matrix, extract the preset camera single-frame exposure time window; this exposure time window is determined by the camera's operating parameters and includes an exposure start time accurate to microseconds and an exposure end time accurate to microseconds. The difference between the two is the exposure duration. For typical shipborne optical imaging equipment, this exposure duration is between 1ms and 20ms; for the exposure time window of each image frame, extract all matrices whose timestamps completely fall within the time window from the time series of the deformation projection transformation matrix. These matrices constitute the instantaneous sampling set of the deformation mapping relationship during the exposure of that frame.

[0040] Since the sampling period of high-frequency flexible deformation component sequences is typically 5ms to 10ms, and the exposure time of a single frame is on the same order of magnitude, under typical exposure time (≥10ms) configurations, deformation projection transformation matrices can usually be extracted for 2 to 5 instantaneous sampling moments. If there are fewer than 2 sampling points within the exposure time window, a first-order linear extrapolation method is used to complete the instantaneous sampling matrix, thereby constructing the undersampled aliasing trajectory within the exposure period. This means that the spatial position corresponding to any pixel on the image plane is not fixed between the start and end of the exposure. Instead, it moves continuously in the pixel coordinate system along a broken line segment trajectory determined by the deformation projection transformation matrix of these 2 to 5 discrete sampling times. According to this physical process, for each pixel on the image plane, the offset vector of the nearest neighbor grid node in the deformation projection transformation matrix corresponding to each instantaneous sampling time within the exposure time window is queried in chronological order. The precise pixel position of the pixel at each instantaneous sampling time is obtained by bilinear interpolation. These instantaneous pixel positions are connected in chronological order to form a continuous pixel cumulative offset path.

[0041] This path completely records the entire positional change process of the pixel due to elastic deformation during a single frame exposure. The starting point of the path corresponds to the pixel position at the beginning of the exposure, and the ending point corresponds to the pixel position at the end of the exposure. Performing this path extraction operation on all pixels yields the set of pixel cumulative offset paths for all pixels within the exposure period of that frame. This set is the undersampled aliasing trajectory within the exposure period. The trajectory is called undersampled aliasing because the sampling frequency of the high-frequency flexible deformation component sequence is limited by the output rate of the inertial measurement unit, with a maximum of only 200Hz. However, the local elastic deformation of the hull may contain vibration components as high as 50Hz to 100Hz. According to the Nyquist sampling theorem, the sampling frequency is insufficient to completely capture these high-frequency vibrations, resulting in spectral aliasing when reconstructing the continuous deformation process using discrete sampling points. This causes high-frequency components that cannot be correctly restored due to undersampling to be mixed into the pixel cumulative offset path.

[0042] Step 3.4: Calculate intra-frame pixel gradient diffusion based on the undersampled aliasing trajectory within the exposure period to obtain a pixel-level displacement distortion field characterizing aperiodic phase distortion. The aperiodic phase distortion refers to a pixel-level displacement distortion mode caused by local elastic deformation of the hull and lacking a fixed periodicity. Specifically, it includes: Based on the obtained undersampled aliasing trajectory within the exposure period, intra-frame pixel gradient diffusion calculation is performed on each frame image. Specifically, for each pixel cumulative offset path in the undersampled aliasing trajectory, the starting pixel coordinates and ending pixel coordinates are extracted, and the Euclidean distance between the two points is calculated. This distance quantifies the total net displacement of the pixel due to elastic deformation during a single frame exposure, in pixels. Simultaneously, the displacement vectors formed by all adjacent instantaneous sampling points on the path are extracted, and the root mean square value of the length of these displacement vectors is calculated. This root mean square value reflects the degree of non-stationary jitter of the pixel's displacement process during exposure and is used to characterize the high-frequency fluctuation amplitude of the displacement. The total net displacement is used as the magnitude of the comprehensive displacement distortion vector, the direction of the line connecting the starting point and the ending point is used as the direction of the vector, and the high-frequency fluctuation amplitude is added as the confidence scalar parameter of the vector, together forming a comprehensive displacement distortion vector for the pixel.

[0043] By traversing all pixels within an image frame, a corresponding comprehensive displacement distortion vector is assigned to each pixel, forming a two-dimensional vector field with the same number of rows and columns as the original image frame. In this vector field, the two components of each vector element are along the column and row directions of the pixel coordinate system, respectively. This vector field accurately describes the non-uniform offset direction and magnitude of the corresponding pixel position caused by the local elastic deformation of the hull within the exposure period. The vectors at different pixel positions exhibit spatial gradations in both direction and magnitude, and there is no regular periodic arrangement relationship between the vectors. Therefore, this two-dimensional vector field is a pixel-level displacement distortion field characterizing aperiodic phase distortion. This distortion field provides a pixel-by-pixel accurate distortion compensation basis for subsequent steps to correct the initial rigid body motion sequence of the hull and generate a high-precision steady-state motion vector field.

[0044] In a preferred embodiment of the present invention, step 4 above may include: Step 4.1: Analyze the spatial gradient distribution characteristics of the pixel-level displacement distortion field to obtain the reverse phase compensation instruction set; perform trajectory superposition correction on the time-synchronized initial rigid body motion sequence of the hull according to the reverse phase compensation instruction set to obtain the initial compensated motion trajectory, including: Step 4.11: Calculate the first-order spatial difference along the row and column directions for each vector element of the pixel-level displacement distortion field to obtain the distortion gradient vector field; mark the distortion-sensitive region according to the magnitude distribution of the distortion gradient vector field and the preset gradient threshold; take the inverse vector of the distortion vector for each pixel position within the distortion-sensitive region to generate an inverse phase compensation instruction set, specifically including: The pixel-level displacement distortion field is used as input. This distortion field is a two-dimensional vector field with the same number of rows and columns as the original image frame. Each element is a two-dimensional vector containing row and column offsets. For each vector element of this distortion field, the difference between the offset vectors of the right and left adjacent pixels in the row direction is taken and divided by twice the pixel spacing. The difference between the offset vectors of the lower and upper adjacent pixels in the column direction is taken and divided by twice the pixel spacing to obtain the row and column gradients at each pixel position. The two vectors are summed and merged into a distortion gradient vector. Its magnitude represents the degree of spatial change of the distortion field at that position, and its direction points in the direction of increasing distortion.

[0045] Based on the magnitude distribution of the distortion gradient vectors, pixel regions with magnitudes exceeding a preset gradient threshold are marked as distortion-sensitive regions. The preset gradient threshold is set to the mean of the magnitude distribution of all distortion gradient vectors in the current frame plus 1.5 to 2.5 times the standard deviation, with a typical value range of 0.05 pixels to 0.15 pixels. For each pixel position within the distortion-sensitive region, the inverse vector of its original distortion vector is taken as an inverse phase compensation instruction. The components of this inverse vector have the same value but opposite sign as the components of the original distortion vector, indicating the amount and direction of pixel displacement compensation required to offset the aperiodic phase distortion at that pixel. The inverse phase compensation instructions of all distortion-sensitive regions in the same frame are summarized to obtain the inverse phase compensation instruction set for that frame.

[0046] Step 4.12: The pixel-level compensation instructions in the reverse phase compensation instruction set are back-projected to the three-dimensional object space through the inverse mapping relationship of the preset camera intrinsic parameter matrix, and converted into angular displacement compensation in the camera coordinate system; the angular displacement compensation is then superimposed and corrected with the angular velocity data at the corresponding moment in the initial rigid body motion sequence of the ship synchronized with the time sequence after time differential transformation, to obtain the initial-order compensated motion trajectory, specifically including: Each pixel-level compensation instruction in the inverse phase compensation instruction set is back-projected to the three-dimensional object space through the inverse mapping relationship of the preset camera intrinsic parameter matrix, and converted into a set of angular displacement compensation quantities in the camera coordinate system; the conversion process is as follows: Subtracting the principal point coordinates from the compensated pixel coordinates yields the offset relative to the principal point. Dividing this offset by the ratio of the equivalent focal length to the physical size of the pixel unit, we obtain the normalized directional angular components along the second and first axes in the camera coordinate system. Using the camera's front node as the rotation center, the spatial direction determined by the normalized directional angular components is taken as the corrected object-side line-of-sight direction. An angle difference decomposition is performed between this direction and the nominal object-side line-of-sight direction of the pixel under the reference distortion-free state to obtain the pitch displacement compensation around the second axis and the azimuth displacement compensation around the third axis.

[0047] The angular displacement compensation amounts corresponding to all pixels within the same frame are merged and averaged to obtain a set of comprehensive angular displacement compensation amounts for that frame. The comprehensive angular displacement compensation amounts are then converted by time differentiation to obtain the angular velocity compensation increment, which is superimposed with the three-axis angular velocity sampling values ​​in the corresponding time segment of the initial rigid body motion sequence of the ship synchronized with the time sequence to obtain the corrected three-axis angular velocity sequence. The above operation is performed on all image frames to obtain the initial-order compensated motion trajectory after pixel-level distortion field correction.

[0048] Step 4.2: Extract the angular acceleration abrupt change nodes of the initial-order compensated motion trajectory to obtain the trajectory abrupt change feature set; perform weight mapping calculation based on the trajectory abrupt change feature set and the preset rigid body rotational inertia smoothing constraint to obtain the time-frequency domain joint filtering weight set, specifically including: for each axial angular velocity sequence in the obtained initial-order compensated motion trajectory, calculate the angular acceleration value point by point along the time axis, the angular acceleration value is obtained by dividing the difference in angular velocity between the current sampling point and the immediately preceding sampling point by the sampling period; perform abrupt change node detection on the angular acceleration sequences of the three axes respectively, specifically: set an angular acceleration jump detection threshold, in a preferred embodiment, the value of the angular acceleration jump detection threshold is based on: calculating the angular acceleration sequences of the three axes of the initial-order compensated motion trajectory within the entire time period. The threshold is set to the mean plus 3 to 5 times the standard deviation. For typical shipborne mobile communication terminals in medium to high sea states, the standard deviation of angular acceleration of the roll and pitch axes is usually between 2 radians per second squared and 8 radians per second squared, while the standard deviation of the bow axis is slightly lower. Correspondingly, the typical range of the angular acceleration jump detection threshold is 10 radians per second squared to 30 radians per second squared. This value can effectively identify abnormal jumps caused by high-frequency flutter or sensor noise, and can also avoid misjudging normal wave-induced motion fluctuations as abrupt change nodes. When the absolute value of angular acceleration at a certain sampling point exceeds the threshold and the excess is greater than a preset multiple of the absolute value of angular acceleration at the previous sampling point, the sampling point is marked as an angular acceleration abrupt change node, and the occurrence time, axis position, and jump amplitude of the node are recorded.

[0049] All angular acceleration abrupt change nodes detected along the three axes are merged and arranged in chronological order of occurrence to form the trajectory abrupt change feature set of the initial-order compensated motion trajectory. Based on this, a preset rigid body rotational inertia smoothing constraint is introduced. This constraint is set based on the actual physical characteristics of the shipborne mobile communication terminal. That is, as a rigid body with a certain mass and size distribution, the rotational inertia of the terminal about each axis determines that its angular acceleration cannot undergo infinitely large instantaneous jumps. Any angular acceleration abrupt change exceeding the physically achievable range is considered a high-frequency flutter or noise artifact that needs to be suppressed. Based on the jump amplitude of each abrupt change node in the trajectory abrupt change feature set and the preset rotational inertia value corresponding to its axis position, the suppression weight to be assigned to that abrupt change node in the time-frequency domain joint filtering is calculated. In a preferred embodiment, the formula for calculating the suppression weight is: ; In the formula, For the first The suppression weights for each node with abrupt angular acceleration change have values ​​ranging from 0 to 1. This is the index of the mutation node in the trajectory mutation feature set; It is the first The preset normalized moment of inertia value corresponding to the axis where each mutation node is located is obtained by normalizing the moment of inertia of the three axes of the terminal. The axis with the largest moment of inertia has a value of 1, and the other axes are proportionally converted. The threshold value for the physically achievable maximum angular acceleration increment, determined based on the maximum output torque of the terminal drive motor and the frictional torque of the transmission mechanism, is expressed in radians per second squared. It is the first The absolute value of the angular acceleration jump at each mutation node; it should be noted that when Less than the preset minimum positive number ,like At that time, directly ordered To avoid the result approaching infinity when the denominator approaches zero, division is not performed. This indicates taking the smaller of the two values; the suppression weights of all mutation nodes in the trajectory mutation feature set are reallocated according to the time axis and axial axis to obtain a time-frequency domain joint filtering weight set that is isomorphic to the initial compensated motion trajectory in the time and axial dimensions. The default value of the weights at non-mutation node positions in this weight set is 1.

[0050] Step 4.3: Perform time-frequency domain nonlinear fusion operation on the initial-order compensated motion trajectory according to the time-frequency domain joint filter weight set to obtain the candidate steady-state motion vector sequence after high-frequency flutter suppression. Specifically, this includes: firstly, performing time-frequency transformation on the initial-order compensated motion trajectory along the three axes to obtain the time-frequency amplitude spectrum and time-frequency phase spectrum of each axis; simultaneously, mapping the time-frequency domain joint filter weight set to the time-frequency domain according to the same time-frequency transformation parameters to form the time-frequency weight matrix of each axis; for the time-frequency amplitude spectrum of each axis, a weight-guided nonlinear fusion rule is used for processing: that is, at each time-frequency unit position, the original amplitude value and the reference amplitude value smoothed by low-pass filtering are mixed according to the weight value of that position; in a preferred embodiment, the calculation formula for this nonlinear fusion operation is: ; in, This indicates the fused time-frequency amplitude spectrum at the frequency index. and time frame index The amplitude value at that point; The amplitude value of the original time-frequency amplitude spectrum of the initial compensated motion trajectory at this time-frequency unit; It is the amplitude value of the reference amplitude spectrum at this time-frequency unit after low-pass filtering; This represents the weight value corresponding to the time-frequency unit in the joint filtering weight set of the time-frequency domain. For time-frequency units with a weight value close to 1, the fusion result tends to retain the original amplitude to preserve the high-frequency details of the motion trajectory. For time-frequency units with a weight value close to zero, the fusion result tends to use the smoothed reference amplitude to suppress the high-frequency flutter components at abrupt change nodes. After fusion, the corrected time-frequency amplitude spectrum is recombined with the original time-frequency phase spectrum, and the fused time-frequency spectrum of each axis is mapped back to the time domain through inverse time-frequency transformation to obtain the high-frequency flutter-suppressed angular velocity sequence of each axis. The angular velocity sequences of the three axes are aligned and merged in time to form the candidate steady-state motion vector sequence after high-frequency flutter suppression.

[0051] Step 4.4: Based on the candidate steady-state motion vector sequence after high-frequency flutter suppression, perform spectral continuity verification and phase alignment interpolation to obtain a continuous, non-aliased, high-precision steady-state motion vector field. Specifically, this includes: performing spectral continuity verification and phase alignment interpolation on the candidate steady-state motion vector sequence. The specific operation of spectral continuity verification is: transforming the candidate steady-state motion vector sequence back to the time-frequency domain, and checking along the frequency axis whether there are any sudden changes in the local peak frequency positions of the amplitude spectrum between adjacent time frames exceeding a preset frequency jump threshold. In a preferred embodiment, the value of the preset frequency jump threshold is based on: the output of step 4.3... The residual spectral aliasing in the candidate steady-state motion vector sequence mainly originates from the energy leakage of adjacent frequency bands caused by the imperfect weight transition region during the nonlinear fusion of the frequency domain components of the high-frequency structural response in the time-frequency domain. The peak frequency jump caused by this leakage usually does not exceed twice the analysis frequency resolution. Therefore, the preset frequency jump threshold is set to twice the frequency resolution of the time-frequency transform used in step 4.3. This frequency resolution is equal to the reciprocal of the length of the time-frequency transform analysis window. For a typical analysis window length of 64 sampling points and a sampling rate of 200Hz, the frequency resolution is 3.125Hz, and the corresponding preset frequency jump threshold is 6.25Hz.

[0052] If the detected offset of the peak frequency of the amplitude spectrum between adjacent time frames exceeds the threshold, it is determined that there is residual spectral aliasing in the motion vector of the local time region, and it is repaired by performing local median filtering on the time-frequency units of the region. After the verification is passed, the candidate steady-state motion vector sequence is compared with the time-synchronized initial rigid body motion sequence of the hull generated in step 1 by sampling point to detect the time deviation between the two at the zero-crossing point of the time-domain waveform. For the detected time deviation exceeding the sampling period, Within the sampling point range, a synchronous interpolation algorithm is used to insert or adjust the timing and amplitude of sampling points in the candidate sequence, so that the corrected sequence is precisely aligned with the zero-crossing point of the waveform of the original rigid body motion sequence in phase, eliminating the cumulative phase delay introduced by the aforementioned multi-stage processing. After the above-mentioned spectral continuity verification and phase alignment interpolation processing, the resulting motion vector sequence is phase-continuous without jumps in the time domain and smooth spectral amplitude without aliasing in the frequency domain, which is a continuous and aliased high-precision steady-state motion vector field. This vector field provides a motion compensation benchmark accurate to the sub-pixel level for subsequent steps of inverse geometric transformation and motion blur suppression processing.

[0053] In a preferred embodiment of the present invention, step 5 above may include: Step 5.1: Analyze the inter-frame spatial transformation parameters of the high-precision steady-state motion vector field to obtain the inverse geometric mapping reference set; perform frame-by-frame coordinate remapping operation on the original image dataset of the time-series matching based on the inverse geometric mapping reference set to obtain the geometrically corrected image sequence, including: Step 5.11: Integrate the three-axis angular velocity vectors at two adjacent sampling times in the high-precision steady-state motion vector field over time to obtain inter-frame spatial transformation parameters, which include three-axis rotation vectors; invert the three-axis rotation vectors to obtain inverse geometric mapping reference parameters, and collect them in chronological order to form an inverse geometric mapping reference set; construct an inverse mapping lookup table from the corrected reference frame pixel coordinates to the original image pixel coordinates of the current frame based on each set of parameters in the inverse geometric mapping reference set; based on the inverse mapping lookup table, perform sub-pixel interpolation coordinate remapping operations frame by frame on the temporally matched original image dataset to obtain a geometrically corrected image sequence, specifically including: The high-precision steady-state motion vector field is used as input. This vector field is a continuous and non-aliased sequence of motion vectors on the time axis. Each sampling moment contains a three-axis angular velocity vector, which describes the ideal steady-state pointing change of the camera optical axis after rigid-flexible separation compensation and rotational inertia smoothing constraint. The three-axis angular velocity vectors of two adjacent sampling moments in this vector field are integrated over time to obtain the inter-frame spatial transformation parameters. The inter-frame spatial transformation parameters include a three-axis rotation vector and a three-axis translation vector. The three-axis rotation vector describes the rigid body rotation angle of the camera coordinate system about the three coordinate axes between two frames, and the three-axis translation vector describes the displacement of the origin of the camera coordinate system in three-dimensional space between two frames. For scenarios where the shipborne mobile communication terminal camera and antenna are coaxially mounted and imaging at a long distance, the three-axis translation vector is set to zero.

[0054] Based on the inversion of the three-axis rotation vector in the inter-frame spatial transformation parameters, inverse geometric mapping reference parameters are obtained. The inverse geometric mapping reference parameters of all adjacent frame pairs are then collected in chronological order to form an inverse geometric mapping reference set. Based on each set of parameters in the inverse geometric mapping reference set, an inverse mapping lookup table is constructed, mapping the pixel coordinates of the corrected reference frame to the pixel coordinates of the current frame's original image. This lookup table specifies the corresponding sampling coordinates in the current frame's original image for each pixel position of the corrected reference frame with sub-pixel precision. Based on the inverse mapping lookup table, coordinate remapping operations are performed frame-by-frame on each original image in the time-matched original image dataset. A bicubic interpolation method is used to obtain resampled grayscale values ​​at the sub-pixel coordinate positions and fill them into the corresponding pixel positions of the corrected reference frames, generating geometrically corrected images that correspond one-to-one with the original image frames. All geometrically corrected images are arranged in their original chronological order to obtain a geometrically corrected image sequence.

[0055] The construction process of the inverse mapping lookup table is as follows: Using each integer pixel coordinate of the corrected reference frame as the starting position, the pixel coordinate is subtracted from the principal point coordinate to obtain the pixel offset relative to the principal point. This offset is then divided by the ratio of the equivalent focal length of the camera lens to the physical size of the pixel unit in the row direction and the ratio of the equivalent focal length to the physical size of the pixel unit in the column direction, respectively, to obtain the two-dimensional coordinates of the pixel on the normalized image plane. In the camera coordinate system, with the camera front node as the origin, the normalized image plane coordinates are extended into a three-dimensional direction vector. The first axis component of this direction vector takes the normalized focal length value, and the second and third axis components take the two coordinate values ​​of the normalized image plane coordinates, respectively, thereby determining the spatial orientation of the object-side line of sight corresponding to the pixel in the reference distortion-free state.

[0056] Based on the three-axis rotation vector in the inverse geometric mapping reference parameters corresponding to the current frame pair, a three-dimensional rotation matrix is ​​constructed. This three-dimensional rotation matrix is ​​applied to the three-dimensional direction vector to obtain a new direction vector after rotation compensation. This new direction vector describes the actual spatial orientation of the same object-side line of sight in the camera coordinate system at the exposure time of the original image of the current frame. The second and third axis components of the new direction vector are divided by the first axis component to obtain the normalized image plane coordinates. These coordinates are then multiplied by the ratio of the equivalent focal length to the physical size of the pixel unit and the principal point coordinates are added to obtain the sub-pixel sampling coordinates corresponding to the corrected reference frame pixel in the original image of the current frame.

[0057] Traverse all pixel positions of the corrected reference frame and perform the above mapping operation one by one. Record the mapping relationship between each pixel of the corrected reference frame and its corresponding sub-pixel sampling coordinates. In other words, construct the inverse mapping lookup table from the pixel coordinates of the corrected reference frame to the pixel coordinates of the original image of the current frame.

[0058] Step 5.2: Analyze the inter-frame instantaneous gradient changes of the high-precision steady-state motion vector field to obtain the residual velocity component; extract the local pixel grayscale gradient features of the geometrically corrected image sequence to obtain the intra-frame spatial frequency distribution set; construct a spatial degradation point diffusion model based on the intra-frame spatial frequency distribution set and the residual velocity component, and perform adaptive frequency domain deconvolution operation to obtain a set of corrected frame images, including: Step 5.21: Calculate the first-order forward difference along the time axis for the triaxial angular velocity vectors at each sampling moment in the high-precision steady-state motion vector field to obtain the angular acceleration vector; map the angular acceleration vector to the image plane through a preset camera intrinsic parameter matrix to obtain the residual velocity component; extract the local pixel grayscale gradient features of each frame in the geometrically corrected image sequence to obtain the intra-frame spatial frequency distribution set; determine the trail length and motion direction angle based on the magnitude and direction of the residual velocity component, construct a linear point spread function and its corresponding optical transfer function, and form a spatial degradation point spread model; based on the spatial degradation point spread model, perform frequency domain Wiener filtering deconvolution operation on the geometrically corrected image sequence frame by frame to obtain a set of corrected frame images; specifically including: For each sampling moment in the high-precision steady-state motion vector field, the first-order forward difference is calculated along the time axis to obtain the rate of change of angular velocity between adjacent sampling moments, i.e., the angular acceleration vector; the formula for calculating the first-order forward difference is: ; in, For the first The angular acceleration vector calculated at each sampling time point is a three-dimensional vector containing three scalar components; It is the first in a high-precision steady-state motion vector field The three-axis angular velocity vector at each sampling time; For the first The three-axis angular velocity vector at each sampling time; The time step between two adjacent sampling moments in a high-precision steady-state motion vector field; the angular acceleration vector The image is projected onto the image plane using the mapping relationship of the preset camera intrinsic parameter matrix, and converted into the residual velocity component at that moment in the pixel coordinate system. The residual velocity component is a two-dimensional vector, whose row direction component and column direction component represent the pixel motion velocity remaining on the image plane due to incomplete compensation of angular acceleration.

[0059] Local pixel grayscale gradient features are extracted from each frame of the geometrically corrected image sequence to obtain an intra-frame spatial frequency distribution set. Based on this, the residual velocity component and the intra-frame spatial frequency distribution set are jointly used to construct a spatial degradation point spread model. The core description of this model is: during a single-frame exposure, the pixel grayscale value undergoes a linear tailing along the motion direction driven by the residual velocity component. The tailing length is equal to the magnitude of the residual velocity component multiplied by the exposure duration, and the tailing direction is consistent with the direction of the residual velocity component. This tailing effect is spatially equivalent to a convolution operation between the original sharp image and a linear point spread function. The specific expression of the linear point spread function is: ; In the formula Represents coordinates in the image plane space The point spread function value at which, The coordinates are in the row direction spatial coordinates. The coordinates are in column direction space. This indicates the trail length, which is equal to the magnitude of the residual velocity component multiplied by the exposure time of a single frame, and is expressed in pixels. The motion direction angle representing the residual velocity component is defined as the angle between the residual velocity vector and the positive direction of the axis of motion, in radians; Let represent the Dirac impulse function, which takes the value of infinity when the independent variable is 0 and the value of 0 when the independent variable is non-zero, and its integral result is 1; This represents a rectangular window function, where the absolute value of its independent variable is less than... When the value is 1, the absolute value equals The time value is Absolute value greater than The value is 0 at that time; and These are the cosine and sine functions, respectively; the corresponding representation of the spatial degradation point spread model in the frequency domain is the optical transfer function, which is obtained by a two-dimensional Fourier transform of the linear point spread function; in a preferred embodiment, the specific expression of the optical transfer function is: ; In the formula Represents the coordinates in the frequency domain The optical transfer function value at the location, where For the spatial frequency in the row direction, The column-direction spatial frequency is represented by periods per pixel. The singer function is defined as follows: when the independent variable z≠0, When the independent variable z=0, The limit is defined as 1, i.e., sinc(0) = 1; and The meaning is the same as the definition in the linear point spread function expression; based on this optical transfer function, an adaptive frequency domain deconvolution operation is performed on each frame in the geometrically corrected image sequence. Specifically, the image is divided into blocks and then a two-dimensional Fourier transform is performed on each block. In the frequency domain, the image is divided by the optical transfer function corresponding to the block and multiplied by a Wiener filter correction factor to suppress noise amplification. Then, the image is returned to the spatial domain through an inverse Fourier transform to obtain the corrected frame image after removing residual motion blur.

[0060] Based on the optical transfer function, an adaptive frequency domain deconvolution operation is performed on each frame of the geometrically corrected image sequence: the image is divided into blocks, and a two-dimensional Fourier transform is performed on each block. In the frequency domain, the block is divided by the corresponding optical transfer function and multiplied by a Wiener filter correction factor to suppress noise amplification. Then, an inverse Fourier transform is performed to return to the spatial domain, resulting in a corrected frame image after removing residual motion blur. The Wiener filter correction factor is a frequency domain weighting term adaptively adjusted according to the local signal-to-noise ratio (SNR) of the image. Its value at each frequency coordinate is the ratio of the image signal power spectrum to the sum of the signal power spectrum and the noise power spectrum. When the SNR at that frequency is high, the Wiener filter correction factor approaches 1, and the deconvolution is approximately a direct inverse filter. When the SNR is low, the Wiener filter correction factor approaches 0, and the deconvolution operation is suppressed to avoid excessive noise amplification. All corrected frame images are arranged in chronological order to obtain a set of corrected frame images.

[0061] The specific construction process of the spatial degradation point diffusion model is as follows: The geometric parameters of the spatial degradation point diffusion model are determined by the residual velocity components. The trailing length, in pixels, is obtained by multiplying the magnitude of the residual velocity component by the single-frame exposure duration; the motion direction angle, in radians, is obtained by taking the angle between the residual velocity vector and the positive direction of the line axis. The trailing length and motion direction angle together define the linear trailing trajectory of the pixel grayscale value along the motion direction during a single-frame exposure.

[0062] The frequency domain adaptive partitioning parameters of the spatial degradation point diffusion model are determined by the intra-frame spatial frequency distribution set. The extraction process of the intra-frame spatial frequency distribution set is as follows: For each frame in the geometrically corrected image sequence, an analysis window of a preset size slides pixel by pixel along the row and column directions. At each window position, the two-dimensional gradient magnitude of the pixel grayscale value within the window is calculated. The calculation formula is as follows: with the center pixel of the window as the reference, the center difference gradient in the row and column directions is calculated respectively. The square root of the sum of the squares of the gradients in the two directions is used to obtain the spatial frequency response value of the window position. The spatial frequency response values ​​of all window positions are sorted according to the coordinates of the center pixel of the window to form the intra-frame spatial frequency distribution set of the frame. The value of each element in this distribution set represents the richness of image details in the local neighborhood of the corresponding pixel. The larger the value, the higher the spatial frequency of the region and the richer the edge and texture information.

[0063] Based on this, the intra-frame spatial frequency distribution set and the residual velocity components are jointly integrated into the construction of the spatial degradation point spread model: a linear point spread function is constructed based on the tail length and motion direction angle. The linear point spread function is then subjected to a two-dimensional Fourier transform to obtain the optical transfer function; the optical transfer function is then spatially adaptively weighted and corrected using the intra-frame spatial frequency distribution set. Specifically: The geometrically corrected image is divided into high-frequency texture regions and low-frequency flat regions based on the local statistical characteristics of the intra-frame spatial frequency distribution set. Within the intra-frame spatial frequency distribution set, pixel regions with spatial frequency response values ​​greater than the average of all spatial frequency response values ​​in the frame are classified as high-frequency texture regions, and the rest are classified as low-frequency flat regions. For high-frequency texture regions, the optical transfer function is directly adopted to fully recover the high-frequency details lost due to motion blur. For low-frequency flat regions, a regularized offset term is applied to the optical transfer function to appropriately increase the amplitude of the optical transfer function near the zero frequency point in this region, thereby suppressing the noise amplification effect caused by the denominator approaching zero in subsequent deconvolution operations.

[0064] The final spatial degradation point diffusion model consists of the following three parts: The model employs a geometric degradation parameter comprised of the tail length and motion direction angle determined by the residual velocity components, a spatial partitioning mapping of high-frequency texture regions and low-frequency flat regions determined by the intra-frame spatial frequency distribution set, and an optical transfer function adaptively corrected for each partition. In subsequent adaptive frequency-domain deconvolution operations, the model applies Wiener filtering deconvolution to the high-frequency texture regions and low-frequency flat regions respectively, using the corresponding partitioned optical transfer functions. This effectively suppresses noise amplification in flat regions while restoring image details weakened by motion blur, thus balancing deblurring effect with image signal-to-noise ratio.

[0065] Step 5.3: Extract the inter-frame key target spatial residual from the set of calibration frame images to obtain the target pose deviation sequence; calculate the offset between the actual pointing of the terminal optical axis and the preset beam tracking reference based on the target pose deviation sequence to obtain the pointing error fluctuation sequence. Specifically, this includes: selecting a preset key target region in each frame image according to the set of calibration frame images. This key target region is usually the circular edge contour region of the satellite antenna parabolic reflector in the image or the obvious corner feature region at the end of the feed support rod; performing edge detection and contour fitting on the image within the key target region to extract the precise pixel position coordinates of the key target in the current frame; comparing the pixel position of the key target in the current frame with the pixel position of the same key target in the preset reference frame; calculating the difference in pixel coordinates between the two in the row and column directions. This difference is the inter-frame key target spatial residual.

[0066] The inter-frame key target spatial residuals of all adjacent frame pairs in the entire set of calibrated frame images are arranged in chronological order to obtain the target pose deviation sequence. This sequence quantifies the residual pixel-level jitter amplitude of key targets in the imaging field of view on the time axis after geometric correction and motion blur suppression. Based on this, each deviation value in the target pose deviation sequence is back-projected onto the three-dimensional object space through the inverse mapping relationship of the preset camera intrinsic parameter matrix. Combined with the coaxial installation geometry of the camera and the satellite antenna, the angular offset of the actual pointing of the terminal optical axis relative to the preset beam tracking reference at that moment is calculated. The preset beam tracking reference is the ideal pointing of the antenna beam calculated in real time based on the satellite ephemeris and ship position attitude. The angular offset is decomposed into two components along the azimuth axis and the elevation axis, and arranged in chronological order to obtain the pointing error fluctuation sequence.

[0067] Step 5.4: Perform multi-level threshold closed-loop verification based on the pointing error fluctuation sequence to obtain a frame sequence temporal rearrangement instruction; perform dynamic frame sequence reorganization on the corrected frame image set according to the instruction to obtain the stable video sequence, specifically including: the multi-level threshold system includes three levels: the first level is the warning threshold, set to the half-power beamwidth of the satellite beam. When the azimuth or elevation error component of a frame in the pointing error fluctuation sequence exceeds the warning threshold, it is determined that the beam pointing accuracy of that frame has slightly deteriorated, and the frame is marked but not removed; the second level is the error correction threshold, which is set to the half-power beamwidth of the satellite beam. When the error component of a frame exceeds the error correction threshold, it is determined that the image of that frame has deviated significantly from the beam pointing. At this time, the frame sequence temporal reordering mechanism is triggered. While maintaining temporal continuity, the frame is locally ordered with the preceding and following frames to smooth the inter-frame abrupt changes in pointing error.

[0068] The third level is the rejection threshold, set to the half-power beamwidth of the satellite beam. When the error component of a frame exceeds the rejection threshold, the frame is determined to no longer meet the minimum requirements for satellite communication beam alignment. The frame is then directly removed from the correction frame image set, and interpolation is performed using the preceding and following valid frames. It should be noted that the above three-level threshold ratios are based on ITU-R. The graded requirements for pointing accuracy of shipborne mobile communication antennas in Recommendation S.465-6 are determined by ROC curve analysis and optimization based on the measured data of this system. The aforementioned multi-level threshold closed-loop verification generates a frame sequence time-series rearrangement instruction, which records in detail the frame indices that need to be marked, rearranged, or removed, and their corresponding processing operation types. According to this instruction, the set of corrected frame images is dynamically reordered. Marked frames are retained but with added quality marking information. The time position of rearranged frames in the output sequence is adjusted. Removed frames are replaced with the interpolation results of the preceding and following valid frames, thereby obtaining a stable video sequence that meets the dual constraints of the preset motion blur radius and satellite beam pointing error. The preset motion blur radius is pre-set according to the satellite communication link budget and image quality requirements, with a typical value of 1 pixel. This stable video sequence has clear and unblurred edge details and stable and jitter-free inter-frame transitions in terms of image quality, and ensures that the antenna beam is always accurately aligned with the target satellite throughout the entire video period in terms of communication performance.

[0069] like Figure 2 As shown, embodiments of the present invention also provide an image stabilization processing system for shipborne mobile satellite communication, comprising: The timing alignment module is used to perform timing alignment processing based on sampling clock deviation on the synchronously acquired three-axis angular velocity signals and the original image frame sequence to obtain the timing-synchronized initial rigid body motion sequence of the hull, and extract the corresponding image data to obtain the timing-matched original image dataset. The rigid-flexible frequency domain separation module is used to perform multi-scale fluid-structure rigid-flexible frequency domain separation operations on the initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum. Based on the mixed frequency domain spectrum, the low-frequency rigid body swaying component is filtered out to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. The distortion field analysis module is used to construct the optomechanical coupling spatial projection mapping relationship based on the high-frequency flexible deformation component sequence, obtain the deformation projection transformation matrix, analyze the undersampled aliasing trajectory within the exposure period, and obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. The motion trajectory fusion module is used to correct the initial rigid body motion sequence of the hull according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory. The first-order compensated motion trajectory is nonlinearly fused in the time and frequency domain using a preset rigid body rotation inertia smoothing constraint to obtain a continuous and non-aliased high-precision steady-state motion vector field. The image correction and closed-loop verification module is used to perform inverse geometric transformation and motion blur suppression processing on the original image dataset of the time-matched model frame by frame according to the high-precision steady-state motion vector field to obtain a set of corrected frame images, and then perform closed-loop verification of pointing error to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

[0070] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0071] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. An image stabilization processing method for shipborne mobile satellite communication, characterized in that, The method includes: The synchronously acquired three-axis angular velocity signals and the original image frame sequence are subjected to time alignment processing based on sampling clock deviation to obtain a time-synchronized initial rigid body motion sequence of the hull; the corresponding image data is extracted according to the initial rigid body motion sequence of the hull to obtain a time-matched original image dataset. Based on the initial rigid body motion sequence of the hull, a multi-scale fluid-structure rigid-flexible frequency domain separation operation is performed to obtain a mixed frequency domain spectrum; based on the mixed frequency domain spectrum, the low-frequency rigid body swaying component is filtered out to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. Based on the high-frequency flexible deformation component sequence, an optomechanical coupled spatial projection mapping relationship is constructed to obtain the deformation projection transformation matrix; based on the deformation projection transformation matrix, the undersampled aliasing trajectory within the exposure period is analyzed to obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. The initial rigid body motion sequence of the hull is corrected according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory; the first-order compensated motion trajectory is nonlinearly fused in the time and frequency domain using a preset rigid body rotation inertia smoothing constraint to obtain a continuous and non-aliased high-precision steady-state motion vector field. The original image dataset of the time-matched sequence is subjected to inverse geometric transformation and motion blur suppression processing frame by frame based on the high-precision steady-state motion vector field to obtain a set of corrected frame images; the pointing error closed-loop verification is performed based on the set of corrected frame images to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

2. The image stabilization processing method for shipborne mobile satellite communication according to claim 1, characterized in that, The synchronously acquired triaxial angular velocity signals and the original image frame sequence are subjected to timing alignment processing based on sampling clock deviation to obtain the timing-synchronized initial rigid body motion sequence of the hull. Based on the initial rigid body motion sequence of the hull, corresponding image data is extracted to obtain a time-matched original image dataset, including: The local time-frequency energy distribution characteristics of the triaxial angular velocity signal are analyzed, and a dynamic convolution constraint matrix is ​​constructed based on the local time-frequency energy distribution characteristics; Path optimization calculation is performed based on the dynamic curl constraint matrix and the sampling clock deviation to obtain the adaptive time remapping trajectory. The triaxial angular velocity signal is subjected to variable step-size interpolation based on the adaptive time remapping trajectory to obtain a resampled angular velocity sequence. Extract the frame time boundary of the resampled angular velocity sequence and the exposure window of the original image frame sequence; perform phase alignment matching on the frame time boundary and the exposure window to obtain the time-synchronized initial rigid body motion sequence of the hull; The original image data is extracted from the initial rigid body motion sequence of the ship in time synchronization, and then encapsulated by frame-level data serialization to obtain the original image dataset with time matching.

3. The image stabilization processing method for shipborne mobile satellite communication according to claim 2, characterized in that, Based on the initial rigid body motion sequence of the hull, multi-scale fluid-structure rigid-flexible frequency domain separation operations are performed to obtain a mixed frequency domain spectrum; based on the mixed frequency domain spectrum, low-frequency rigid body rocking components are filtered out to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation, including: The initial rigid body motion sequence of the ship's hull, which is synchronized in time, is subjected to multi-scale fluid-structure rigid-flexible frequency domain expansion operation to obtain a mixed frequency domain spectrum. The multi-scale fluid-structure rigid-flexible frequency domain separation operation refers to using multiple sets of bandpass filter banks with different time resolutions and frequency resolutions to expand the initial rigid body motion sequence of the ship's hull at multiple time and frequency scales in order to separate the low-frequency rigid body swaying component excited by wave fluid loads from the high-frequency flexible deformation component of the elastic response of the ship structure. The spectral energy distribution gradient of the mixed frequency domain spectrum is analyzed, and the motion component within the wave fluid load excitation frequency band is extracted based on the spectral energy distribution gradient to obtain the low-frequency interference characteristic spectrum; A frequency domain dynamic isolation mask is constructed based on the low-frequency interference characteristic spectrum. The low-frequency energy removal operation is performed on the mixed frequency domain spectrum through the frequency domain dynamic isolation mask to obtain the frequency domain component of the high-frequency structural response. Inverse time-domain phase reconstruction is performed based on the frequency domain components of the high-frequency structural response to obtain a sequence of high-frequency flexible deformation components characterizing local elastic deformation.

4. The image stabilization processing method for shipborne mobile satellite communication according to claim 3, characterized in that, Based on the high-frequency flexible deformation component sequence, an optomechanical coupled spatial projection mapping relationship is constructed to obtain the deformation projection transformation matrix; Based on the deformation projection transformation matrix, the undersampled aliasing trajectory within the exposure period is analyzed to obtain the pixel-level displacement distortion field characterizing the aperiodic phase distortion, including: The spatial deformation vector of the high-frequency flexible deformation component sequence is analyzed. Based on the spatial deformation vector, a camera optical axis deflection mapping model is established using the optomechanical coupling spatial projection mapping relationship. The optomechanical coupling spatial projection mapping relationship refers to the rigid connection and transmission relationship between the optical imaging system and the antenna mechanical support structure in the shipborne terminal, which couples and maps the angular displacement of the mounting base caused by the local elastic deformation of the hull into the deflection attitude of the camera optical axis. Then, through camera projection geometry, the optical axis deflection is mapped into a spatial coordinate transformation relationship of pixel displacement on the image plane. Spatial coordinate mapping calculations are performed using the camera optical axis deflection mapping model to obtain the spatial coordinate set of the terminal's local deformation. The deformation projection transformation matrix is ​​obtained by performing a two-dimensional image plane projection operation based on the local deformation space coordinate set of the terminal and the preset camera intrinsic parameters; Extract the preset single-frame exposure time window of the camera, calculate the pixel cumulative offset path within the exposure time window based on the deformation projection transformation matrix, and obtain the undersampling aliasing trajectory within the exposure period; Intra-frame pixel gradient diffusion is calculated based on the undersampled aliasing trajectory within the exposure period to obtain a pixel-level displacement distortion field characterizing aperiodic phase distortion; wherein, the aperiodic phase distortion refers to a pixel-level displacement distortion mode caused by local elastic deformation of the hull and without a fixed periodicity.

5. The image stabilization processing method for shipborne mobile satellite communication according to claim 4, characterized in that, The initial rigid body motion sequence of the hull is corrected based on the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory; the first-order compensated motion trajectory is nonlinearly fused in the time and frequency domain using a preset rigid body rotational inertia smoothing constraint to obtain a continuous, non-aliased, high-precision steady-state motion vector field, including: The spatial gradient distribution characteristics of the pixel-level displacement distortion field are analyzed to obtain the reverse phase compensation instruction set; the trajectory superposition correction of the time-synchronized initial rigid body motion sequence of the hull is performed according to the reverse phase compensation instruction set to obtain the first-order compensated motion trajectory. Extract the angular acceleration mutation nodes of the initial compensated motion trajectory to obtain the trajectory mutation feature set; perform weight mapping calculation based on the trajectory mutation feature set and the preset rigid body rotational inertia smoothing constraint to obtain the time-frequency domain joint filtering weight set; Based on the time-frequency domain joint filter weight set, a time-frequency domain nonlinear fusion operation is performed on the initial compensated motion trajectory to obtain a candidate steady-state motion vector sequence after high-frequency flutter suppression. Based on the candidate steady-state motion vector sequence after high-frequency flutter suppression, spectral continuity verification and phase alignment interpolation are performed to obtain a continuous and non-aliased high-precision steady-state motion vector field.

6. The image stabilization processing method for shipborne mobile satellite communication according to claim 5, characterized in that, Based on the high-precision steady-state motion vector field, the original image dataset of the time-matched sequence is subjected to inverse geometric transformation and motion blur suppression processing frame by frame to obtain a set of corrected frame images; based on the set of corrected frame images, a pointing error closed-loop verification is performed to obtain a stable video sequence that satisfies the preset motion blur radius and satellite beam pointing error constraints, including: The inter-frame spatial transformation parameters of the high-precision steady-state motion vector field are analyzed to obtain the inverse geometric mapping reference set; the frame-by-frame coordinate remapping operation is performed on the original temporally matched image dataset based on the inverse geometric mapping reference set to obtain the geometrically corrected image sequence; The inter-frame instantaneous gradient of the high-precision steady-state motion vector field is analyzed to obtain the residual velocity component; the local pixel gray-level gradient features of the geometrically corrected image sequence are extracted to obtain the intra-frame spatial frequency distribution set; a spatial degradation point diffusion model is constructed based on the intra-frame spatial frequency distribution set and the residual velocity component, and adaptive frequency domain deconvolution operation is performed to obtain the corrected frame image set. The spatial residuals of key targets between frames are extracted from the set of corrected frame images to obtain the target pose deviation sequence; the offset between the actual pointing of the terminal optical axis and the preset beam tracking reference is calculated based on the target pose deviation sequence to obtain the pointing error fluctuation sequence. Multi-level threshold closed-loop verification is performed based on the pointing error fluctuation sequence to obtain the frame sequence temporal reordering instruction; according to the instruction, the set of corrected frame images is dynamically reordered to obtain the stable video sequence.

7. The image stabilization processing method for shipborne mobile satellite communication according to claim 6, characterized in that, The spatial gradient distribution characteristics of the pixel-level displacement distortion field are analyzed to obtain a reverse phase compensation instruction set. Based on the reverse phase compensation instruction set, the initial rigid body motion sequence of the time-synchronized hull is subjected to trajectory superposition correction to obtain the initial-order compensated motion trajectory, including: The first-order spatial difference is calculated for each vector element of the pixel-level displacement distortion field along the row and column directions to obtain the distortion gradient vector field. The distortion sensitive region is marked according to the magnitude distribution of the distortion gradient vector field and the preset gradient threshold. The inverse vector of the distortion vector of each pixel position in the distortion sensitive region is taken to generate the inverse phase compensation instruction set. The pixel-level compensation commands in the reverse phase compensation command set are back-projected to the three-dimensional object space through the inverse mapping relationship of the preset camera intrinsic parameter matrix, and converted into angular displacement compensation in the camera coordinate system. The angular displacement compensation is then superimposed and corrected with the angular velocity data at the corresponding moment in the initial rigid body motion sequence of the ship synchronized with the time sequence after time differential transformation, so as to obtain the first-order compensated motion trajectory.

8. The image stabilization processing method for shipborne mobile satellite communication according to claim 7, characterized in that, The inter-frame spatial transformation parameters of the high-precision steady-state motion vector field are analyzed to obtain an inverse geometric mapping reference set; based on the inverse geometric mapping reference set, a frame-by-frame coordinate remapping operation is performed on the original temporally matched image dataset to obtain a geometrically corrected image sequence, including: Time integration is performed on the three-axis angular velocity vectors at two adjacent sampling times in the high-precision steady-state motion vector field to obtain the inter-frame spatial transformation parameters, which include the three-axis rotation vectors. The three-axis rotation vectors are then inverted to obtain the inverse geometric mapping reference parameters, which are then collected in chronological order to form the inverse geometric mapping reference set. Based on the parameters of each set in the inverse geometric mapping reference set, an inverse mapping lookup table is constructed from the pixel coordinates of the corrected reference frame to the pixel coordinates of the original image in the current frame. Based on the inverse mapping lookup table, sub-pixel interpolation coordinate remapping operation is performed frame by frame on the temporally matched original image dataset to obtain the geometrically corrected image sequence.

9. The image stabilization processing method for shipborne mobile satellite communication according to claim 8, characterized in that, The inter-frame instantaneous gradient of the high-precision steady-state motion vector field is analyzed to obtain the residual velocity component; local pixel grayscale gradient features of the geometrically corrected image sequence are extracted to obtain the intra-frame spatial frequency distribution set; a spatial degradation point diffusion model is constructed based on the intra-frame spatial frequency distribution set and the residual velocity component, and adaptive frequency domain deconvolution is performed to obtain a set of corrected frame images, including: The first-order forward difference of the three-axis angular velocity vectors at each sampling time in the high-precision steady-state motion vector field is calculated along the time axis to obtain the angular acceleration vector; the angular acceleration vector is then mapped to the image plane through a preset camera intrinsic parameter matrix to obtain the residual velocity components. Local pixel grayscale gradient features are extracted from each frame in the geometrically corrected image sequence to obtain the intra-frame spatial frequency distribution set; the trail length and motion direction angle are determined based on the magnitude and direction of the residual velocity component, and a linear point spread function and its corresponding optical transfer function are constructed to form a spatial degradation point spread model. Based on the spatial degradation point diffusion model, frequency domain Wiener filtering deconvolution operation is performed on the geometrically corrected image sequence frame by frame to obtain a set of corrected frame images.

10. An image stabilization processing system for shipborne satellite communication in motion, the system implementing the method as described in any one of claims 1 to 9, characterized in that, include: The timing alignment module is used to perform timing alignment processing based on sampling clock deviation on the synchronously acquired three-axis angular velocity signals and the original image frame sequence to obtain the timing-synchronized initial rigid body motion sequence of the hull, and extract the corresponding image data to obtain the timing-matched original image dataset. The rigid-flexible frequency domain separation module is used to perform multi-scale fluid-structure rigid-flexible frequency domain separation operations on the initial rigid body motion sequence of the hull to obtain a mixed frequency domain spectrum. Based on the mixed frequency domain spectrum, the low-frequency rigid body swaying component is filtered out to obtain a high-frequency flexible deformation component sequence characterizing local elastic deformation. The distortion field analysis module is used to construct the optomechanical coupling spatial projection mapping relationship based on the high-frequency flexible deformation component sequence, obtain the deformation projection transformation matrix, analyze the undersampled aliasing trajectory within the exposure period, and obtain the pixel-level displacement distortion field characterizing the non-periodic phase distortion. The motion trajectory fusion module is used to correct the initial rigid body motion sequence of the hull according to the pixel-level displacement distortion field to obtain the first-order compensated motion trajectory. The first-order compensated motion trajectory is nonlinearly fused in the time and frequency domain using a preset rigid body rotation inertia smoothing constraint to obtain a continuous and non-aliased high-precision steady-state motion vector field. The image correction and closed-loop verification module is used to perform inverse geometric transformation and motion blur suppression processing on the original image dataset of the time-matched model frame by frame according to the high-precision steady-state motion vector field to obtain a set of corrected frame images, and then perform closed-loop verification of pointing error to obtain a stable video sequence that meets the preset motion blur radius and satellite beam pointing error constraints.

Images