A method for suppressing noise of a detection signal of a laser wind-radar

Abstract

The present application relates to the technical field of signal processing, in particular to a detection signal noise suppression method of a laser wind measurement radar. Based on the spatial physical redundancy of multi-range gate observation of the laser radar, the present application utilizes the physical characteristics that the platform rigid body motion is a full field common mode and the atmospheric turbulence is a spatial differential mode, to solve the technical problem that the real wind field and the platform noise are difficult to separate under the condition of frequency domain aliasing. By calculating the spatial frequency standard deviation sequence, the present application uses the redundancy of spatial information to eliminate the ambiguity of pure time dimension analysis. This technical means enables the motion inversion algorithm to adaptively adjust the trajectory smoothing weight coefficient according to the environmental turbulence intensity, automatically reduces the inertia constraint to retain the high frequency wind speed characteristics at strong turbulence time, and enhances the constraint to smooth the observation noise at weak turbulence time, so as to realize the accurate stripping of the platform common mode motion trajectory under complex weather conditions, effectively suppress the detection signal noise, and improve the wind speed inversion precision.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of signal noise processing, and particularly relates to a detection signal noise suppression method of a laser wind radar. BACKGROUND

[0002] The laser wind radar inverts the wind field velocity by emitting a laser beam and receiving the backscattering echo of atmospheric aerosols. When the radar is carried on a dynamic platform such as a ship, a buoy or an airborne platform, the high-frequency rigid-body motion of the platform will introduce additional Doppler frequency shift in the echo signal. In order to realize long-distance detection, the radar usually uses non-coherent accumulation technology to improve the signal-to-noise ratio of weak signals. However, during the accumulation period, the rapid motion of the platform will cause the Doppler frequency of the echo signal at different times to drift. This non-stationarity causes the energy of the accumulated signal spectrum to be dispersed and widened, resulting in the signal-to-noise ratio being unable to be effectively improved with the accumulation time, which seriously restricts the long-range detection capability of the radar in a dynamic environment.

[0003] The inventor found in practice that the prior art has the following defects:

[0004] The prior art usually relies on an inertial measurement unit (IMU) for motion compensation, but in low-cost applications or severe shaking conditions, the measurement residual of the IMU will still remain in the signal. Extracting motion parameters from the radar echo signal itself is an effective supplementary means, but in a real atmospheric environment, the wind field itself has turbulent changes. The rapid wind speed fluctuations caused by turbulence and the velocity changes caused by platform motion are highly mixed in the frequency domain characteristics, and the traditional low-pass filtering method cannot distinguish between the two, which is easy to misjudge the real wind shear as platform noise for filtering, or fail to completely remove the platform interference, resulting in a decrease in the accuracy of the final wind speed inversion. SUMMARY

[0005] In order to solve the technical problem that the prior art cannot distinguish between the rapid wind speed fluctuations caused by turbulence and the velocity changes caused by platform motion in the frequency domain characteristics, the purpose of the present application is to provide a detection signal noise suppression method of a laser wind radar, and the technical solution adopted is as follows:

[0006] The present application provides a detection signal noise suppression method of a laser wind radar, which comprises the following steps:

[0007] Obtaining the original echo signal of the radar in the full observation period, short-time slicing the original echo signal and performing spectrum transformation to obtain a power spectral density matrix;

[0008] reference distance gates are screened out from the carrier-to-noise ratio of each initial distance gate in all short-time slices; a frequency fluctuation matrix is obtained based on the frequency fluctuation of each reference distance gate in each short-time slice; a hypothesis platform frequency shift grid is constructed based on the frequency fluctuation matrix; a spatial frequency standard deviation sequence is determined according to the additive motion characteristics of the frequency distribution platform of all reference distance gates of each short-time slice and the spatial non-uniform motion characteristics of turbulence;

[0009] a cost function of a dynamic programming algorithm is defined based on the hypothesis platform frequency shift grid, the cost function including a frequency shift mutation penalty; a trajectory smoothing weight coefficient in the dynamic programming algorithm is determined based on the spatial frequency standard deviation sequence; and a globally optimal platform common mode frequency shift estimation sequence is searched for by performing dynamic programming recursion and backtracking according to the frequency shift mutation penalty in the hypothesis platform frequency shift grid by using the trajectory smoothing weight coefficient;

[0010] measured spectral data is extracted from the power spectral density matrix and corrected by using the platform common mode frequency shift estimation sequence, and high-precision wind speed data free from platform motion interference is output.

[0011] Further, the power spectral density matrix is obtained in the following manner:

[0012] For data in each short-time slice, a windowing method is used to suppress spectral sidelobe leakage.

[0013] Fast Fourier transform is performed on the windowed data, the square of the modulus is calculated, and a Doppler power spectral density is obtained; all Doppler power spectral densities form the power spectral density matrix.

[0014] Further, the reference distance gates are obtained in the following manner:

[0015] The average carrier-to-noise ratio of each initial distance gate in all short-time slices is calculated.

[0016] All initial distance gates are sorted in descending order of the average carrier-to-noise ratio value, and a plurality of initial distance gates at the top of the sorting are selected to form a reference distance gate set.

[0017] Further, the frequency fluctuation matrix is obtained in the following manner:

[0018] For each reference distance gate in the reference distance gate set, a Doppler center frequency is extracted from the spectral data of each short-time slice to obtain an instantaneous frequency value.

[0019] The frequency mean value of the reference distance gate in the entire observation period is calculated, the instantaneous frequency value is subtracted from the frequency mean value to obtain a frequency fluctuation value after removing the direct current component, and the frequency fluctuation values of all short-time slices and reference distance gates are combined to construct a frequency fluctuation matrix.

[0020] Further, the manner of obtaining the hypothesis platform frequency shift grid comprises:

[0021] Windowing the data in each short-time slice, performing fast Fourier transform on the windowed data to determine the frequency resolution;

[0022] Determining the maximum possible frequency shift boundary according to the platform physical motion limit calibrated by the inertial navigation unit;

[0023] Setting the frequency resolution as the step size, dividing a plurality of discrete states in the maximum possible frequency shift boundary interval to construct the hypothesis platform frequency shift grid, and the hypothesis platform frequency shift grid comprises a plurality of discrete state hypothesis platform common mode frequency shifts.

[0024] Further, the manner of obtaining the spatial frequency standard deviation sequence comprises:

[0025] Calculating the normalized signal-to-noise ratio weight of each reference range gate in the reference range gate set;

[0026] Based on the normalized weight, calculating the frequency weighted mean of all reference range gates of each short-time slice, and then calculating the weighted spatial standard deviation of each short-time slice, and the weighted spatial standard deviations of all short-time slices form the spatial frequency standard deviation sequence.

[0027] Further, the manner of obtaining the trajectory smoothing weight coefficient comprises:

[0028] Calculating the global standard deviation of all elements in the frequency fluctuation matrix, and setting it as the reference turbulence reference value;

[0029] Taking the ratio of the elements of the spatial frequency standard deviation sequence and the reference turbulence reference value and the preset minimum positive number as the denominator, and taking the preset basic inertia constant as the numerator, calculating the fraction to obtain the trajectory smoothing weight coefficient.

[0030] Further, the manner of obtaining the platform common mode frequency shift estimation sequence comprises:

[0031] Based on the frequency shift mutation penalty and the hypothesis platform frequency shift grid, recursively and backtracking using the Viterbi algorithm to obtain the globally optimal platform common mode frequency shift estimation sequence.

[0032] Further, the manner of obtaining the frequency shift mutation penalty comprises:

[0033] Calculating the square of the trajectory smoothing weight coefficient and the hypothesis platform common mode frequency shift difference value of adjacent short-time slices to obtain the frequency shift mutation penalty.

[0034] Further, the manner of obtaining the high-precision wind speed data after removing the platform motion interference comprises:

[0035] Extracting target spectrum data from the power spectrum density matrix; correcting the target spectrum data of the platform common mode frequency shift estimation sequence frame by frame; obtaining the motion compensated cumulative spectrum based on the corrected spectrum data, and performing spectrum peak search based on the cumulative spectrum to determine the frequency index corresponding to the maximum capability peak, and calculating the final wind speed data by using the Doppler calculation formula.

[0036] The present application has the following beneficial effects:

[0037] Based on the spatial physical redundancy of the multi-range gate observation of the laser radar, the present application solves the technical problem that the real wind field and the platform noise are difficult to separate under the condition of frequency domain aliasing by using the physical characteristics that the platform rigid body motion is a common mode of the whole field and the atmospheric turbulence is a spatial differential mode. By calculating the standard deviation sequence of the spatial frequency, the present application uses the redundancy of the spatial information to eliminate the ambiguity of the pure time dimension analysis. This technical means enables the motion inversion algorithm to adaptively adjust the trajectory smoothing weight coefficient according to the environmental turbulence intensity, automatically reduces the inertia constraint to retain the high-frequency wind speed characteristics at the strong turbulence moment, and enhances the constraint to smooth the observation noise at the weak turbulence moment, thereby realizing the accurate stripping of the platform common mode motion trajectory under complex weather conditions, effectively suppressing the detection signal noise, and improving the wind speed inversion precision. BRIEF DESCRIPTION OF DRAWINGS

[0038] In order to more clearly illustrate the technical solutions and advantages of the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiment or prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and those skilled in the art can also obtain other drawings according to these drawings without any creative effort.

[0039] Figure 1 A flow chart of a detection signal noise suppression method of a laser wind measuring radar provided by an embodiment of the present application. DETAILED DESCRIPTION

[0040] In order to further illustrate the technical means and effects adopted by the present application to achieve the predetermined invention purpose, the following will combine the drawings and the preferred embodiments to specifically describe the specific implementation, structure, features and effects of the detection signal noise suppression method of the laser wind measuring radar according to the present application, and the detailed description is as follows. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. In addition, the specific features, structures or characteristics in one or more embodiments can be combined in any suitable form.

[0041] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the present application belongs.

[0042] The application provides a detection signal noise suppression method of a laser wind-radar.

[0043] Embodiment 1

[0044] The application provides a detection signal noise suppression method of a laser wind-radar. Figure 1 Fig. 1 shows a flowchart of a detection signal noise suppression method of a laser wind-radar according to an embodiment of the application, and the method comprises the following steps.

[0045] In step S101, original echo signals of a radar in a full observation period are acquired, and the original echo signals are subjected to short-time slicing and spectrum transformation to acquire a power spectrum density matrix.

[0046] The application scenario of the embodiment of the application can be a scenario of laser wind measurement on a dynamic platform such as a ship, a buoy or an airborne platform, original signals can be collected by a laser wind-radar, and the original signals are subjected to signal preprocessing, and after signal processing, original echo signals are obtained and input into a memory for calling and acquisition.

[0047] The original echo signals are continuous long-time observation signals, in order to capture the high-frequency transient motion characteristics of the platform, the continuous long-time observation signals need to be decomposed into a series of short-time observation units.

[0048] In step S101, original echo signals of a radar in a full observation period are acquired, and the original echo signals are subjected to short-time slicing and spectrum transformation to acquire a power spectrum density matrix.

[0049] In an implementation manner of the embodiment of the application, the time length T of the short-time slicing is usually 1 / 10 to 1 / 20 of the platform motion period, for example, 50 ms, and the sliding step length t is usually 50% of T to maintain data continuity.

[0050] In step S101, original echo signals of a radar in a full observation period are acquired, and the original echo signals are subjected to short-time slicing and spectrum transformation to acquire a power spectrum density matrix.

[0051] In step S102, a reference distance gate is screened out based on the carrier-to-noise ratio of each initial distance gate in all short-time slices; a frequency fluctuation matrix is acquired based on the frequency fluctuation of each reference distance gate in each short-time slice; a hypothetical platform frequency shift grid is constructed based on the frequency fluctuation matrix; and a spatial frequency standard deviation sequence is determined according to the frequency distribution of all reference distance gates of each short-time slice.

[0052] The initial distance gate refers to all spatial sampling points divided by the radar on the full detection path, which includes both the near-field strong signal region and the far-field weak signal region. The initial distance gate is filtered to obtain a reference distance gate.

[0053] The reference distance gate is a high carrier-to-noise ratio subset selected from all distance gates and is usually located in the near-field region. These distance gates have strong echo signals and are less affected by noise, and can accurately reflect the frequency fluctuations caused by platform motion and atmospheric turbulence. They are used as reference observation sources for inverting the common-mode motion trajectory of the platform to decouple and compensate for motion errors in the far-field weak signal.

[0054] Due to the attenuation characteristics of laser in atmospheric transmission, the echo signal-to-noise ratios of different distance gates differ significantly. Low signal-to-noise ratio signals are greatly affected by shot noise and cannot accurately reflect the fine motion of the platform, while strong echo signals and weak signals are modulated by platform rigid body motion simultaneously. Therefore, high signal-to-noise ratio signals need to be selected as references for inverting platform motion. The average carrier-to-noise ratio of each distance gate in all short-time slices is calculated by traversing the Doppler power spectral density matrix. All distance gates are sorted in descending order of average carrier-to-noise ratio values. A plurality of distance gates at the top of the sorting are selected to form a reference distance gate set.

[0055] The frequency fluctuation matrix is a two-dimensional space-time data matrix representing the instantaneous frequency AC change characteristics of the radar near-field high signal-to-noise ratio region after removing the average wind speed DC component. Each row of the matrix represents a short-time slice, and each column represents a selected reference distance gate.

[0056] The absolute Doppler frequency measured by the radar includes a DC component generated by the average wind speed and an AC component generated by the turbulence and platform motion. The platform motion is only reflected as an AC fluctuation that changes over time. In order to separate the motion characteristics, the DC frequency bias of each distance gate needs to be removed.

[0057] For each reference distance gate in the reference distance gate set, the Doppler center frequency is extracted from the spectrum data of each short-time slice to obtain the instantaneous frequency value. The frequency average of the reference distance gate over the entire observation period is calculated. The instantaneous frequency value is subtracted from the frequency average to obtain the frequency fluctuation value after removing the DC component. The frequency fluctuation values of all short-time slices and reference distance gates are combined to construct a frequency fluctuation matrix. This matrix only contains AC change information reflecting platform motion and atmospheric turbulence, and is used as input data for subsequent motion inversion algorithms.

[0058] The platform frequency shift grid is a preset state space constructed for discrete optimal path search in this scheme.

[0059] The Doppler shift caused by the platform motion is a continuously changing analog quantity in physics, and directly performing nonlinear trajectory optimization on the observed data containing random noise and turbulence disturbance in the continuous domain is a mathematical problem with extremely high computational complexity or even unsolvable; therefore, the continuous frequency shift space is first discretized into a platform frequency shift grid based on the maximum dynamic range of the platform. Through the above operation, the complex trajectory inversion is converted into a calculable discrete state search problem, and the necessary normalization scale is provided for the adaptive inertia constraint mechanism, thereby ensuring the feasibility and robustness of the motion inversion algorithm.

[0060] wherein the spatial frequency standard deviation sequence is a time sequence constructed for quantifying the spatial non-uniformity of the atmospheric wind field at each time.

[0061] Since the flight time of the laser pulse is much smaller than the change period of the platform motion and the turbulence, it can be considered that within the same short time slice, the platform motion exerts an equal Doppler shift (i.e. common-mode interference) on all distance gates. According to the shift-invariance principle of variance, the dispersion degree of the frequency distribution of the full-field reference distance gate is not affected by the common-mode shift, but only depends on the spatial non-uniformity of the atmospheric turbulence. Therefore, calculating the spatial standard deviation can effectively separate the atmospheric turbulence intensity from the platform mixed motion.

[0062] In an implementation manner of the embodiment of the application, the reference distance gate screening range is usually 5 to 10 distance gates ranked in the front according to the average carrier-to-noise ratio, so as to ensure the robustness of the statistics.

[0063] In an implementation manner of the embodiment of the application, the extraction method of the Doppler center frequency is usually a spectrum peak search method or a barycenter method.

[0064] In an implementation manner of the embodiment of the application, the maximum dynamic range of the platform can be determined by obtaining the maximum possible frequency shift boundary through the platform physical motion limit calibrated by the inertial navigation unit.

[0065] Step S103: based on the assumed platform frequency shift grid, a cost function of a dynamic programming algorithm is defined, the cost function containing a frequency shift mutation penalty; based on the spatial frequency standard deviation sequence, a trajectory smoothing weight coefficient in the dynamic programming algorithm is determined; and by using the trajectory smoothing weight coefficient, a global optimal platform common-mode frequency shift estimation sequence is searched by performing dynamic programming recursion and backtracking according to the frequency shift mutation penalty in the assumed platform frequency shift grid.

[0066] wherein the frequency shift mutation penalty is a cost function item for constraining the smoothing degree of the calculated trajectory in the dynamic programming algorithm of the scheme, and is intended to reflect the physical inertia of the object motion.

[0067] In order to recover the platform motion trajectory conforming to the physical law from the observation data containing noise, a cost evaluation system containing motion continuity constraints needs to be established.

[0068] Among them, the motion continuity constraint cannot be fixed: it should be enhanced to suppress noise under weak turbulence, and it should be relaxed to allow the algorithm to track the real wind speed mutation under strong turbulence. Further, in order to solve the problem of discontinuity of single-point estimation, a dynamic programming algorithm is used to find a globally optimal path with the minimum cumulative cost on the entire time axis. The path is the platform motion trajectory.

[0069] Among them, the trajectory smoothing weight coefficient is the core control parameter in this scheme for dynamically adjusting the inertia constraint strength of the motion inversion algorithm.

[0070] The purpose of calculating the trajectory smoothing weight coefficient is to dynamically adjust the constraint strength of the motion inversion algorithm on the frequency mutation according to the real-time perceived environmental turbulence intensity, in order to balance noise smoothing and high-frequency feature tracking. The algorithm is given environmental adaptive ability, solving the problem that real wind shear and platform noise are difficult to separate under frequency domain aliasing, ensuring that the platform motion trajectory can be accurately inverted under complex working conditions such as strong turbulence or severe shaking, thereby significantly improving the compensation accuracy and signal-to-noise ratio of far-field weak signals.

[0071] Among them, the platform common mode frequency shift estimation sequence is the final motion inversion result output after spatio-temporal joint solution.

[0072] The purpose of calculating the platform common mode frequency shift estimation sequence is to accurately separate the additive Doppler shift component caused only by the platform rigid body motion from the mixed echo disturbed by atmospheric turbulence. Its beneficial effect is to provide a high-precision motion compensation reference, so that the system can eliminate the spectrum broadening and drift caused by the dynamic platform through reverse frequency rearrangement without external sensors, thereby realizing coherent focusing of signal energy and significantly improving the far-field detection distance and wind speed inversion accuracy of the radar in a dynamic environment.

[0073] Using the trajectory smoothing weight coefficient, dynamic programming recursion and backtracking are performed according to the frequency shift mutation penalty in the assumed platform frequency shift grid, and a globally optimal platform common mode frequency shift estimation sequence is searched.

[0074] Step S104: Extract the to-be-tested frequency spectrum data from the power spectrum density matrix, correct it using the platform common mode frequency shift estimation sequence, and output high-precision wind speed data free of platform motion interference.

[0075] From the power spectrum density matrix after time-frequency transformation, extract the to-be-tested frequency spectrum data to construct the to-be-corrected target spectrum matrix. At this time, due to the modulation of the dynamic platform non-stationary motion, the spectral peak energy is dispersed and broadened in the frequency domain, and the signal-to-noise ratio is very low.

[0076] With the platform common mode frequency shift estimation sequence as the compensation reference, the Doppler frequency shift amount caused by platform motion at each time is calculated, and the target spectrum matrix is accordingly subjected to frame-by-frame inverse frequency shift. This operation aims to build a servo coordinate system in the frequency domain, and correct and align the originally time-drifting signal energy to an equivalent stationary base frequency position.

[0077] The calibrated spectrum data is subjected to coherent accumulation, and the linear focusing of energy and the improvement of signal-to-noise ratio are realized by using the signal phase alignment feature. The spectrum peak center frequency is searched based on the accumulated high signal-to-noise ratio spectrum, the line-of-sight wind speed data excluding the platform motion interference component is obtained, and high-precision far-field detection in a dynamic environment is realized.

[0078] After the processing of steps S101 to S104, the current high-precision wind speed data excluding platform motion interference is obtained. In the process of calculating the wind speed data, the collected data is first subjected to short-time slicing and time-frequency domain conversion. Considering the significant difference in echo signal-to-noise ratio of different distance gates, the carrier-to-noise ratio is used to select a reference distance gate from the initial distance gates. Since the absolute Doppler frequency measured by the radar contains a direct current component generated by the average wind speed and an alternating current component generated by the turbulence and platform motion, the frequency fluctuation matrix is obtained according to the frequency fluctuation characteristics of the reference distance gate in each short-time slice. Before the dynamic programming algorithm is executed, since the dynamic programming algorithm cannot directly perform path optimization in a continuous infinite numerical space, and the subsequent adaptive weight calculation lacks a unified quantization standard, the continuous frequency shift space is discretized into a hypothetical platform frequency shift grid based on the maximum dynamic range of the platform, that is, a reasonable path optimization space is limited to ensure the rationality of the scheme. The discrete degree of the frequency distribution of the full-field reference distance gate is not affected by the common mode frequency shift, but only depends on the spatial non-uniformity of atmospheric turbulence. Therefore, calculating the spatial standard deviation can effectively separate the atmospheric turbulence intensity from the platform mixed motion. In order to recover the platform motion trajectory conforming to the physical law from the observation data containing noise, a cost evaluation system containing motion continuity constraints is established, and a dynamic programming algorithm is used to find a globally optimal path with the minimum cumulative cost on the entire time axis. Finally, the collected spectrum data is corrected by using the platform common mode frequency shift estimation sequence, and high-precision wind speed data excluding platform motion interference is obtained.

[0079] Preferably, in some possible implementation manners of the embodiments of the present application, the power spectrum density matrix is obtained in the following manner:

[0080] Due to the finite length time domain truncation of non-integer period signal, the discontinuity occurs at the boundary of FFT period extension, which leads to the spectrum side lobe leakage; this phenomenon can significantly raise the side lobe level and broaden the main lobe, so that the side lobe of strong signal covers the adjacent weak signal, and the system dynamic range is seriously reduced. Therefore, the time domain windowing method is needed to smooth the signal edge to eliminate the boundary mutation, so as to effectively suppress the side lobe leakage and improve the weak signal detection ability.

[0081] For the data in each short-time slice, the Hanning window or Hamming window is applied to suppress the spectrum side lobe leakage. The fast Fourier transform is performed on the windowed data, and the square of the modulus is calculated to obtain the Doppler power spectrum density. All Doppler power spectrum densities constitute the power spectrum density matrix.

[0082] Preferably, in some possible implementation manners of the embodiment of the application, the reference distance gate acquisition method comprises:

[0083] Due to the non-uniformity of atmospheric aerosol distribution and the geometric attenuation effect of laser transmission, the echo signals of the laser radar at different detection distances have significantly different carrier-to-noise ratios. The echo signals in the far field or low aerosol area are weak and are seriously polluted by shot noise and speckle noise, and the extracted instantaneous frequency has a great random error, which cannot accurately reflect the fine motion characteristics of the platform. The echo signals of all distance gates are physically modulated by the rigid body motion of the platform, which is common mode synchronous.

[0084] Step S201: calculating the average carrier-to-noise ratio of each initial distance gate in all short-time slices.

[0085] The average carrier-to-noise ratio is selected to filter the reference distance gate, so as to select the near-field data with the smallest noise pollution and the highest signal quality from the full-field echo as the observation reference. Only the high carrier-to-noise ratio signal can accurately reflect the fine frequency fluctuation of the platform motion, and the measurement noise can be reduced to the maximum extent to interfere with the turbulence frequency characteristics, so as to ensure the physical authenticity and calculation accuracy of the subsequent motion inversion algorithm.

[0086] Step S202: sorting all initial distance gates according to the average carrier-to-noise ratio value from large to small, and selecting a plurality of initial distance gates at the front of the sorting to form a reference distance gate set.

[0087] According to the sorting from large to small according to the average carrier-to-noise ratio, the reference source as the input of the algorithm is ensured to be the data with the best quality and the most reliable in the whole observation period.

[0088] The reference distance gates are obtained by screening the initial distance gates and are usually located in the near-field region. These distance gates have strong echo signals and are less affected by noise, and can reflect the frequency fluctuations caused by platform motion and atmospheric turbulence with high fidelity. They are used as a reference observation source for inverting the common-mode motion trajectory of the platform to decouple and compensate for motion errors in the weak signals in the far field.

[0089] The embodiment ranks the average carrier-to-noise ratio in the entire observation period to select several reference distance gates with high confidence. By using the high signal-to-noise ratio advantage of the near-field strong echo signal to replace the full-field data for motion feature extraction, the influence of observation noise on frequency estimation accuracy is significantly reduced, ensuring that the frequency fluctuation matrix constructed subsequently can retain the true details of platform motion with high fidelity, thereby laying a low-noise data foundation for high-precision motion trajectory inversion.

[0090] It should be noted that in one specific embodiment of the present embodiment, the reference distance gate is obtained by traversing the power spectral density matrix, calculating the average carrier-to-noise ratio (CNR) of each distance gate in all short-time slices, and sorting all distance gates in descending order of average CNR value. The top N distance gates are selected to form a reference distance gate set, where N is an integer greater than or equal to 1. In all short-time slices. All distance gates are sorted in descending order of average CNR value. The top N distance gates are selected to form a reference distance gate set, where N is an integer greater than or equal to 1. The value of N is usually in the range of 5 to 10 to ensure the robustness of the statistics. The elements in the set are the selected reference distance gate indexes, where i is the sequence number of the reference distance gate in the set, and N is the number of reference distance gates in the set. The elements in the set are the selected reference distance gate indexes, where i is the sequence number of the reference distance gate in the set, and N is the number of reference distance gates in the set. The elements in the set are the selected reference distance gate indexes, where i is the sequence number of the reference distance gate in the set, and N is the number of reference distance gates in the set. The elements in the set are the selected reference distance gate indexes, where i is the sequence number of the reference distance gate in the set, and N is the number of reference distance gates in the set. The elements in the set are the selected reference distance gate indexes, where i is the sequence number of the reference distance gate in the set, and N is the number of reference distance gates in the set.

[0091] It should be noted that in one specific embodiment of the present embodiment, the reference distance gate is obtained by traversing the power spectral density matrix, calculating the average carrier-to-noise ratio (CNR) of each distance gate in all short-time slices, and sorting all distance gates in descending order of average CNR value. The top N distance gates are selected to form a reference distance gate set, where N is an integer greater than or equal to 1.

[0092] Preferably, in some possible implementations of the present embodiment, the frequency fluctuation matrix is obtained in the following manner:

[0093] The frequency fluctuation matrix is a two-dimensional space-time data matrix that represents the instantaneous frequency AC variation characteristics of the radar near-field high signal-to-noise ratio region after removing the average wind speed direct current component. Each row of the matrix represents a short-time slice, and each column represents a preferred reference distance gate.

[0094] The absolute Doppler frequencies measured by radar include a DC component generated by the average wind speed and an AC component generated by turbulence and platform motion, with the platform motion manifesting only as time-varying AC fluctuations. To separate the motion characteristics, the DC frequency bias of each range gate needs to be removed. This decouples the DC background wind field from the AC disturbance signal in the frequency domain, eliminates the bias introduced by wind shear, and allows the residual fluctuation quantities in the matrix to accurately characterize the time-varying features of the platform's rigid body motion and atmospheric turbulence. This provides precise zero-mean observation input for subsequent trajectory inversion based on spatial common-mode characteristics.

[0095] Specifically, for each reference range gate in the reference range gate set, the Doppler center frequency is extracted from the spectral data of each short-time slice to obtain the instantaneous frequency value. The frequency mean of this reference range gate over the entire observation period is calculated: the instantaneous frequency value is subtracted from the frequency mean to obtain the frequency fluctuation value after removing the DC component. The frequency fluctuation values ​​of all short-time slices and reference range gates are combined to construct a frequency fluctuation matrix. This matrix contains only AC variation information reflecting platform motion and atmospheric turbulence, and serves as the input data for subsequent motion inversion algorithms.

[0096] It should be noted that, in one specific embodiment of this example, the frequency fluctuation matrix is ​​obtained as follows:

[0097] For each reference range gate in the reference range gate set In each short-time slice Spectrum data The Doppler center frequency is extracted. Common extraction methods include the spectral peak search method or the centroid method, to obtain the instantaneous frequency value. .

[0098] Calculate the reference distance gate over the entire observation period. frequency mean within .

[0099] Instantaneous frequency value Subtract the frequency mean This yields the frequency fluctuation value after removing the DC component. All A short slice, A frequency fluctuation matrix is ​​constructed by combining the frequency fluctuation values ​​of each reference distance gate. .

[0100] The frequency fluctuation matrix The dimension is Elements in the matrix Characterized the first The reference distance gate at the 1st Doppler frequency fluctuations within a short-time slice.

[0101] Preferably, in some possible implementations of the embodiments of the present invention, it is assumed that the method for obtaining the platform frequency shift grid includes:

[0102] In this scheme, it is assumed that the platform frequency shift grid is a preset state space constructed for discretized optimal path search.

[0103] Since the Doppler frequency shift caused by platform motion is physically a continuously changing analog quantity, directly performing nonlinear trajectory optimization on observation data containing random noise and turbulence disturbances in the continuous domain is a mathematical problem with extremely high computational complexity, or even unsolvable. Furthermore, the dynamic programming approach used in this scheme is essentially a globally optimal path search algorithm based on discrete state space. Therefore, this embodiment, based on the maximum dynamic range and frequency resolution of the radar system, discretizes the continuous Doppler frequency shift interval into a finite number of state nodes, constructing a hypothetical platform frequency shift grid. This transforms the continuous function fitting problem into a shortest path search problem in graph theory. While ensuring the frequency quantization accuracy meets application requirements, it significantly reduces the algorithm's search dimensionality and computational overhead, making it possible to obtain the globally optimal platform motion trajectory with limited computational resources.

[0104] It should be noted that, in one specific implementation of this embodiment, the platform frequency shift grid is obtained as follows: the maximum possible frequency shift boundary is determined based on the platform's physical motion limits calibrated by the inertial navigation unit. Frequency resolution determined by FFT The step size is within the interval Inner division A discrete state is used to construct a hypothetical platform frequency shift grid. .

[0105] The physical motion limit of the platform refers to the maximum instantaneous line-of-sight velocity boundary that the platform can achieve under the constraints of its own dynamic characteristics and external environmental loads (such as sea state). This parameter is usually determined based on the platform's structural design specifications, seakeeping indicators, or the statistical distribution of historical measured data under typical operating conditions. For example, the heave velocity limit of an ocean buoy under severe sea conditions is usually set within ±5 m / s.

[0106] The maximum possible frequency shift boundary is based on the Doppler effect principle, mapping the aforementioned physical velocity limit to a search space constraint value in the frequency domain. This boundary defines the effective state space of the motion inversion algorithm, ensuring that the solution results conform to physical reality. For conventional 1.5μm band shipborne or buoy radar systems, the value of this frequency shift boundary is usually set between ±2MHz and ±10MHz, with the specific value depending on the platform type and the maximum expected sea state.

[0107] Preferably, in some possible implementation manners of the embodiment of the present application, the acquisition manner of the spatial frequency standard deviation sequence comprises:

[0108] The spatial frequency standard deviation sequence is a time sequence constructed for quantifying the spatial non-uniformity of the atmospheric wind field at each moment.

[0109] Since the Doppler shift in the radar echo contains both the platform rigid body motion and the atmospheric turbulence, and both are highly aliasing (both are high-frequency random fluctuations) in the spectral characteristics of the single time dimension, the traditional time-domain filtering or frequency-domain low-pass method cannot strip the platform motion interference while preserving the true high-frequency wind field characteristics. Based on the physical fact, the platform rigid body motion as an additive common mode only causes the overall shift of the frequency distribution, without changing its dispersion, while the atmospheric turbulence as a spatial non-uniform quantity directly leads to the increase of the dispersion of the frequency distribution.

[0110] Therefore, the embodiment utilizes the observation redundancy in the spatial dimension to calculate the spatial frequency standard deviation sequence of the full-field reference distance gate frequency distribution at each moment. This operation constructs a turbulence intensity quantification index independent of the time dimension, successfully realizes the independent perception and decoupling of the atmospheric turbulence under unknown platform motion state, and provides a key environmental prior criterion for subsequent adaptive adjustment of the inertial constraint weight in the dynamic programming algorithm.

[0111] The normalized signal-to-noise ratio weight of each reference distance gate in the reference distance gate set is calculated Based on the normalized weight , the frequency weighted mean of all reference distance gates of each short-time slice is calculated, and then the weighted spatial standard deviation of each short-time slice is calculated. The weighted spatial standard deviations of all short-time slices constitute the elements in the spatial frequency standard deviation sequence, and all the elements are combined to form the spatial frequency standard deviation sequence.

[0112] It should be noted that in the embodiment of the present application, the signal ratio weight of all reference distance gates can be maximum minimum normalized by maximum minimum normalization processing, so that the value range is between [0, 1].

[0113] It should be noted that in one specific embodiment of the present application, the acquisition manner of the spatial frequency standard deviation sequence is:

[0114] For the data in the th row of the frequency fluctuation matrix (where ), based on the normalized weight , the weighted mean of the reference distance gate frequency at this moment is calculated .

[0115] Then, the weighted spatial standard deviation at this moment is calculated, that is, the element in the spatial frequency standard deviation sequence .

[0116] The above calculation is repeated for all short-time slices to generate a complete spatial frequency standard deviation sequence . The greater the mean value of the sequence, the worse the spatial consistency of the wind field at the current moment (the stronger the turbulence); the smaller the mean value, the more stable the wind field.

[0117] Preferably, in some possible implementation manners of the embodiment of the application, the trajectory smoothing weight coefficient is obtained in the following manner:

[0118] The trajectory smoothing weight coefficient is a core control parameter in the scheme for dynamically adjusting the inertial constraint strength of the motion inversion algorithm.

[0119] Since the dynamic radar platform has non-stationary motion characteristics under different sea states or weather conditions: under a stationary sea state, the platform motion is slow, and observation noise dominates, and the smoothing constraint needs to be enhanced; and under a severe sea state (strong turbulence), the platform often accompanies high-frequency violent shaking, and forced smoothing will cause the algorithm to produce lag errors and lose high-frequency motion details. The traditional fixed-parameter filter cannot meet the needs of these two opposite scenarios.

[0120] Therefore, by establishing an adaptive mapping mechanism based on the spatial frequency standard deviation, the trajectory smoothing weight coefficient in the dynamic programming algorithm is calculated in real time. The inertial constraint term in the state transition cost function is dynamically adjusted: the weight is automatically reduced to relax the constraint in strong turbulence to allow the algorithm to track the rapidly changing platform trajectory; the weight is increased to suppress observation noise in weak turbulence. Thus, the optimal trade-off between denoising ability and dynamic tracking ability is achieved in a non-stationary environment.

[0121] It should be noted that, in one specific implementation manner of the embodiment, the trajectory smoothing weight coefficient is obtained in the following manner: first, a reference turbulence reference value is calculated:

[0122] The reference turbulence reference value is a statistical normalization factor for normalizing the turbulence intensity measure.

[0123] The global standard deviation of all elements in the frequency fluctuation matrix is calculated, which is set as the reference turbulence reference value .

[0124] Further calculation is performed using the spatial frequency standard deviation sequence and the reference turbulence reference value .

[0125] wherein, a base inertia constant (typically in the range 10 to 100, depending on the sampling rate of the radar system and the mechanical inertia characteristics of the platform), a regulation sensitivity index (typically taken in the range 0.1 to 0.9), , a very small positive number to prevent division by zero errors (typically taken in the range 10-6 to 10-8), The formula makes the weight automatically decrease when the local turbulence is significantly higher than the reference value , thus reducing the penalty for frequency jumps.

[0126] Preferably, in some possible implementations of embodiments of the application, the way in which the sequence of platform common-mode frequency shifts is obtained comprises:

[0127] wherein the sequence of platform common-mode frequency shifts is the final motion inversion result output after joint space-time resolution.

[0128] Since the non-stationary motion of a dynamic platform introduces time-varying Doppler shift noise in the radar echo signal, causing energy broadening and peak ambiguity of the signal spectrum after long-time accumulation, the signal-to-noise ratio and detection range of a far-field weak signal are severely reduced; and the platform motion modulation received by different range gates is common-mode, so an accurate time sequence trajectory describing the common-mode motion must be obtained to uniformly correct the full-field data. Therefore, a dynamic programming algorithm is used to search for a globally optimal sequence of platform common-mode frequency shifts in a discrete state space. This sequence, as a high-precision parameterized representation of the instantaneous radial velocity of the platform, not only removes the environmental turbulence interference from the mixed signal, but also provides an accurate inverse compensation reference for subsequent processing, making it possible to construct a servo coordinate system in the frequency domain and eliminate the Doppler non-stationarity caused by platform motion.

[0129] In one specific implementation of embodiments of the application, the sequence of platform common-mode frequency shifts is obtained in the following way:

[0130] First, the frequency shift fitting residual is calculated:

[0131] wherein the frequency shift fitting residual is a measure of the deviation of the actual observation data from the assumption that the platform common-mode frequency shift is at the th time instant. Its calculation method is: .

[0132] The frequency shift fitting residual and the frequency shift jump penalty are input into the Viterbi algorithm for forward recursion calculation, and the specific steps are as follows:

[0133] 1. Establish the cumulative cost matrix (dimension ) and the path backtracking pointer matrix​ (dimension ).

[0134] 2. Initialize the cumulative cost matrix: for each state ( ), the cumulative cost depends only on the fitting residual: .

[0135] 3. Recursive computation: for each state at the current time, iterate over all possible states ( ) at the previous time, find the predecessor state that minimizes the previous time cumulative cost + transition cost.

[0136] Compute the minimum path cost: .

[0137] Record the optimal predecessor state index: .

[0138] After the forward recursion is completed, perform a global path backtracking, the specific process is as follows:

[0139] When recursion is performed to the last short observation window , perform a backtracking operation to extract the optimal path.

[0140] At time , find the terminal state index with the minimum cumulative cost: .

[0141] Using the path backtracking pointer matrix , deduce from to , to obtain the optimal state index at each time in turn: .

[0142] Finally, map the optimal state index sequence back to the corresponding frequency value in the hypothesis platform frequency shift grid , that is, at each time m, use the index value to obtain the numerical value in , thereby obtaining the corresponding frequency value data, and further forming the final platform common mode frequency shift estimation sequence .

[0143] Preferably, in some possible implementation manners of the embodiment of the application, the frequency shift mutation penalty obtaining manner comprises:

[0144] Wherein, the frequency shift mutation penalty is a cost function item used to constrain the smoothness of the calculated trajectory in the dynamic programming algorithm of the scheme, and is intended to reflect the physical inertia of the object motion.

[0145] In order to recover the platform motion trajectory conforming to the physical law from the observation data containing noise, a cost evaluation system containing motion continuity constraint needs to be established.

[0146] Wherein, the motion continuity constraint cannot be fixed: it should be enhanced to suppress noise under weak turbulence, and it should be relaxed to allow the algorithm to track the real wind speed mutation under strong turbulence. Further, in order to solve the problem of discontinuity of single-point estimation, a dynamic programming algorithm is used to find a globally optimal path with the minimum cumulative cost on the entire time axis. The path is the platform motion trajectory recovered.

[0147] In a specific implementation of an embodiment of the present application, the calculation method of the frequency shift mutation penalty may be:

[0148] Wherein, is the trajectory smoothness weight coefficient, is the assumed platform common mode frequency shift state at the current time and the previous time, respectively.

[0149] Preferably, in some possible implementation manners of the embodiment of the present application, the method for extracting the to-be-detected frequency spectrum data from the power spectrum density matrix, correcting by using the platform common mode frequency shift estimation sequence, and outputting high-precision wind speed data removed from platform motion interference is:

[0150] The target data to be detected is extracted from the power spectrum density matrix. According to the user-specified detection distance or any distance gate index other than the reference distance gate, the frequency spectrum data of the distance gate in all short-time slices is extracted, and a target frequency spectrum matrix to be corrected is constructed.

[0151] The dimension of the matrix is , and the row vector of the matrix represents the instantaneous power spectrum in the th short-time slice. At this time, due to the existence of platform motion, the spectral peak of the same wind speed target exists severe frequency drift between different time rows, and direct accumulation will cause the dispersion of spectral peak energy.

[0152] The target frequency spectrum is corrected frame by frame by using the platform common mode frequency shift estimation sequence . Since characterizes the additive Doppler shift noise introduced by platform motion, the correction operation is the inverse process of subtracting the frequency shift amount.

[0153] For each short-time slice ( Read the corresponding platform common-mode frequency shift value. Frequency resolution utilized Calculate the discrete frequency index offset corresponding to this moment. :

[0154] in This indicates the rounding operation.

[0155] Next, regarding the first Spectrum vector of the row Perform a cyclic shift operation, treating spectral energy exceeding the boundary as noise truncation or using zero-padding for shifting. Index the spectrum as... Data moved to new location The calculation formula is:

[0156] For out of index range The part uses cyclic boundary handling (i.e., modulo). (operation).

[0157] For all The above operation is repeated for each row of data to generate a corrected spectral matrix. Through this step, the signal energy that originally drifted over time is realigned to a constant frequency index corresponding to its true line-of-sight wind speed.

[0158] The frequency-corrected spectral matrix along the time dimension Perform an arithmetic mean to generate the motion-compensated cumulative spectrum. : .

[0159] During this accumulation process, the amplitude of the signal spectrum peak increases linearly with the number of accumulations, while the amplitude of background noise increases at a lower rate due to its random phase characteristics. Therefore, the accumulated spectrum after motion compensation... It exhibits a significant improvement in signal-to-noise ratio.

[0160] Finally, for Perform a peak search to determine the frequency index corresponding to the maximum energy peak. The final line-of-sight wind speed was calculated using the Doppler formula. :

[0161] in The operating wavelength of the lidar (typically 1000 nm) This wind speed is the high-precision wind measurement result after removing interference from platform motion.

[0162] In summary, the embodiment of the application is based on the spatial physical redundancy of the multi-distance gate observation of the laser radar, utilizes the physical characteristics that the platform rigid body motion is the full field common mode and the atmospheric turbulence is the spatial difference mode, constructs the spatial frequency standard deviation sequence as the independent measurement of the environmental turbulence intensity, further implements the adaptive dynamic programming inversion operation, adjusts the trajectory smoothing weight coefficient in real time according to the turbulence intensity, solves the platform common mode frequency shift estimation sequence considering the physical inertia and high frequency dynamic characteristics in the discrete frequency shift grid, and finally utilizes the sequence to perform the reverse frequency domain rearrangement and coherent accumulation on the far field echo. Through the above technical means, the scheme effectively solves the decoupling problem of the real wind field and the platform noise under the time-frequency aliasing condition, realizes the elimination of the platform motion caused Doppler spectrum broadening without the dependence on external attitude sensors, and significantly improves the signal-to-noise ratio and detection precision of the far field weak signal in the dynamic environment.

[0163] It should be noted that the above-mentioned embodiment sequence of the application is only for description, and does not represent the advantages and disadvantages of the embodiments. The processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multi-task processing and parallel processing are also possible or can be advantageous.

[0164] Each of the embodiments in the specification is described in a progressive manner, and the same and similar parts between the embodiments can be referred to each other. Each embodiment mainly describes the difference from other embodiments.

Claims

1. A method for suppressing noise in the detection signal of a laser wind-measuring radar, characterized in that, The method includes: The raw echo signal of the radar over the entire observation period is obtained, and the power spectral density matrix is ​​obtained by short-time slicing and spectral transformation of the raw echo signal. Reference range gates are selected based on the carrier-to-noise ratio of each initial range gate across all short-time slices; frequency fluctuation matrices are obtained based on the frequency fluctuations of each reference range gate in each short-time slice; a hypothetical platform frequency shift grid is constructed based on the frequency fluctuation matrix; and spatial frequency standard deviation sequences are determined based on the additive motion characteristics of the platform and the spatial non-uniform motion characteristics of turbulence, according to the frequency distribution of all reference range gates in each short-time slice. Based on the hypothetical platform frequency shift grid, a cost function for the dynamic programming algorithm is defined, which includes a frequency shift mutation penalty. Based on the spatial frequency standard deviation sequence, the trajectory smoothing weight coefficient in the dynamic programming algorithm is determined. Using the trajectory smoothing weight coefficient, dynamic programming recursion and backtracking are performed in the hypothetical platform frequency shift grid according to the frequency shift mutation penalty to search for the globally optimal platform common mode frequency shift estimation sequence. The measured spectrum data is extracted from the power spectral density matrix, and corrected using the platform common-mode frequency shift estimation sequence to output high-precision wind speed data free from platform motion interference.

2. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The power spectral density matrix is ​​obtained in the following ways: For the data within each short-time slice, windowing is used to suppress spectral sidelobe leakage; Perform a fast Fourier transform on the windowed data, calculate the square of its modulus, and obtain the Doppler power spectral density; all Doppler power spectral densities constitute the power spectral density matrix.

3. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining the reference distance gate include: Calculate the average carrier-to-noise ratio for each initial distance gate across all short-time slices; Sort all initial range gates in descending order of average carrier-to-noise ratio (CNR) values, and select the top-ranked initial range gates to form a reference range gate set.

4. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The frequency fluctuation matrix is ​​obtained in the following ways: For each reference range gate in the set of reference range gates, the Doppler center frequency is extracted from the spectral data of each short-time slice to obtain the instantaneous frequency value; The frequency mean of the reference range gate over the entire observation period is calculated. The instantaneous frequency value is subtracted from the frequency mean to obtain the frequency fluctuation value after removing the DC component. The frequency fluctuation values ​​of all short-time slices and the reference range gate are combined to construct the frequency fluctuation matrix.

5. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining the hypothetical platform frequency shift grid include: Windowing is applied to the data within each short-time slice, and a fast Fourier transform is performed on the windowed data to determine the frequency resolution. Based on the platform's physical motion limits calibrated by the inertial navigation unit, the maximum possible frequency shift boundary is determined; The frequency resolution is set to the step size, and multiple discrete states are divided within the maximum possible frequency shift boundary interval. A hypothetical platform frequency shift grid is constructed, which contains the hypothetical platform common-mode frequency shift of multiple discrete states.

6. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The spatial frequency standard deviation sequence is obtained through the following methods: Calculate the normalized signal-to-noise ratio weight for each reference range gate in the reference range gate set; Based on normalized weights, the frequency-weighted mean of all reference distance gates for each short-time slice is calculated. Then, the weighted spatial standard deviation of each short-time slice is calculated. The weighted spatial standard deviations of all short-time slices form a spatial frequency standard deviation sequence.

7. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining the trajectory smoothing weight coefficients include: Calculate the global standard deviation of all elements in the frequency fluctuation matrix and set it as the reference turbulence baseline value; The trajectory smoothing weight coefficient is obtained by calculating the fraction using the ratio of the elements of the spatial frequency standard deviation sequence to the reference turbulence baseline value and the preset minimum positive sum value as the denominator, and the preset basic inertial constant as the numerator.

8. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining the platform common-mode frequency shift estimation sequence include: Based on the frequency shift abruptness penalty and the assumed platform frequency shift grid, the Viterbi algorithm is used for recursion and backtracking to obtain the globally optimal platform common mode frequency shift estimation sequence.

9. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining the frequency shift mutation penalty include: The frequency shift abrupt penalty is obtained by multiplying the trajectory smoothing weight coefficient by the square of the hypothetical platform common-mode frequency shift difference between adjacent short-time slices.

10. The method for suppressing noise in the detection signal of a laser wind-measuring radar according to claim 1, characterized in that, The methods for obtaining high-precision wind speed data after removing platform motion interference include: Extract the target spectrum data from the power spectral density matrix; perform frame-by-frame correction using the target spectrum data from the platform common-mode frequency shift estimation sequence; obtain the motion-compensated cumulative spectrum based on the corrected spectrum data, and perform peak search based on the cumulative spectrum to determine the frequency index corresponding to the maximum capability peak; and calculate the final wind speed data using the Doppler calculation formula.

Images