All-weather atmospheric temperature, humidity, wind and rain profile combined detection method and system

Abstract

The invention discloses an all-weather atmospheric temperature, humidity, wind and rain profile combined detection method and system, particularly relates to the technical field of meteorological radar signal processing, and is used for solving the problem of motion parameter inversion failure caused by phase nonlinear distortion of an existing coherent Doppler radar at a dielectric mutation interface. A radar original echo signal is obtained; extracting a phase change sequence and calculating a phase difference value of adjacent range gates; identifying a dielectric abrupt change interface location when the differential value exceeds a linear threshold range; analyzing the continuous period interface position recurrence plot, and extracting diagonal characteristics as a deterministic coefficient and vertical section distribution as a laminar flow parameter; determining a dynamic stability level based on the two, and executing path integral compensation when the dynamic stability level is lower than a threshold value; an original signal phase is corrected through a compensation strategy, high-precision motion parameters are output through coherent accumulation, dynamic compensation of phase distortion in a non-uniform medium is achieved, and robustness of atmospheric parameter inversion is guaranteed.

Description

Technical Field

[0001] The present invention relates to the technical field of weather radar signal processing, and more particularly to an all-weather atmospheric temperature, humidity, wind and rain profile joint detection method and system. Background Art

[0002] Millimeter-wave radar remote sensing of atmospheric media typically uses coherent Doppler processing to obtain target motion information. Existing technologies transmit pulse signals, receive backscattered waves, and invert radial velocity using signal phase changes. When radar waves propagate in a homogeneous medium, the phase changes exhibit a linear accumulation characteristic, enabling high-precision velocity measurement based on standard coherent integration algorithms. This method has been widely used to detect aerosol and particle motion in meteorology and environmental monitoring.

[0003] However, when the radar beam passes through the interface of a sudden medium with a significant gradient in dielectric constant (such as the alternating area between a meteorological target and a non-meteorological medium), the wavefront is distorted due to non-uniform phase shift, causing the phase accumulation characteristics of the echo signal to deviate from the ideal linear model. This phenomenon destroys the prerequisite for coherent processing, causing the spectral energy output of the traditional accumulation algorithm to disperse and the signal-to-noise ratio to deteriorate, ultimately resulting in systematic jump errors in the motion parameter inversion results. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method and system for joint detection of all-weather atmospheric temperature, humidity, wind and rain profiles to solve the problems raised in the above-mentioned background technology.

[0005] To achieve the above object, the present invention provides the following technical solutions: The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method includes the following steps: S1: Obtain the original echo signal received by the radar on the detection path; S2: Extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; S3: When the phase difference value exceeds the preset linear threshold range, the interface position of the dielectric constant mutation is identified; S4: Quantitative analysis of the dielectric constant mutation interface position is performed on the recursive graph for a preset number of consecutive cycles, the diagonal structural characteristics of the recursive graph are extracted as the deterministic coefficient, and the vertical segment distribution is calculated as the laminar flow parameter; S5: determining a dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and executing a path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; S6: Perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing to output high-precision motion parameter inversion results.

[0006] Furthermore, obtaining the original echo signal received by the radar on the detection path includes: The orthogonal dual-channel echo signals on the detection path are collected by the radar receiver; Performing time alignment processing on the orthogonal dual-channel echo signals to generate synchronized orthogonal dual-channel echo signals; Separate the real part and the imaginary part from the synchronous orthogonal dual-channel echo signal to form the in-phase component and the orthogonal component of the original echo signal; The in-phase component and the quadrature component are combined into a complex original echo signal.

[0007] Furthermore, the phase change sequence is extracted from the original echo signal, and the phase difference value of the adjacent range gate is calculated, including: Perform arc tangent operation on the in-phase component and quadrature component of the original echo signal to obtain the phase value of each range gate; Arrange the phase values ​​of all range gates in the order of detection time to form a phase change sequence; Calculate the difference between the phase values ​​of two adjacent range gates in the phase change sequence to obtain a phase difference score; Store the phase difference values ​​in the phase difference sequence array.

[0008] Furthermore, when the phase difference value exceeds a preset linear threshold range, identifying the dielectric constant mutation interface position includes: Traverse each phase difference value in the phase difference sequence array; Performing a double-boundary comparison between the current phase difference value and the upper and lower limits of a preset linear threshold range; When the phase difference value is greater than the upper limit of the preset linear threshold range or less than the lower limit of the preset linear threshold range, the spatial coordinates of the corresponding range gate are recorded; The projection position of the spatial coordinate on the radar detection path is marked as the position of the dielectric constant mutation interface.

[0009] Furthermore, a recursive graph quantitative analysis is performed on the interface positions of dielectric constant mutations for a preset number of consecutive cycles, the diagonal structural characteristics of the recursive graph are extracted as the deterministic coefficient, and the vertical segment distribution is calculated as the laminar flow parameter, including: Arranging the dielectric constant mutation interface positions of a preset number of consecutive cycles in chronological order as an interface position time series; Set the recursive threshold and calculate the Euclidean distance between the interface positions corresponding to every two time points in the interface position time series; When the Euclidean distance is less than the recursion threshold, the recursion point is marked on the corresponding coordinates of the two-dimensional grid to generate a recursion graph; Identify the line segments formed by continuous recursive points in the diagonal direction of the recursive graph, and calculate the characteristic value of the line segment length distribution as the deterministic coefficient; The line segments formed by the vertical continuous recursive points in the recursive graph are detected, and the characteristic values ​​of the line segment length distribution are statistically analyzed as laminar flow parameters.

[0010] Furthermore, a dynamic stability level is determined based on the deterministic coefficient and the laminar flow parameter. When the dynamic stability level is lower than a preset stability threshold, a path integral compensation strategy is executed, including: Input the deterministic coefficient and laminar flow parameters into the preset stability mapping model and output the dynamic stability level value; Read the preset stability threshold based on radar wavelength and medium uniformity calibration; Compare the dynamic stability level value with the preset stability threshold; When the dynamic stability level value is less than the preset stability threshold, the path integral compensation strategy is activated to generate instructions; Generate instructions based on the path integral compensation strategy to initialize the path integral compensation parameter matrix.

[0011] Furthermore, the preset stability mapping model is implemented in the following way: Establish a two-dimensional feature space of deterministic coefficients and laminar flow parameters; Divide the dynamic stability level value into contour areas in the two-dimensional feature space; The dynamic stability level value corresponding to the input deterministic coefficient and laminar flow parameter is calculated by linear interpolation.

[0012] Furthermore, the original echo signal is phase corrected based on the path integral compensation strategy to generate a corrected signal and perform coherent accumulation processing to output high-precision motion parameter inversion results, including: Extracting the range gate compensation vector and the time step compensation vector from the path integral compensation parameter matrix; Performing a Hadamard product operation on the range gate compensation vector and the phase component of the original echo signal to generate a primary compensation signal; Performing a Hadamard product operation on the time step compensation amount vector and the phase component of the primary compensation signal to generate a corrected signal; Performing discrete Fourier transform coherent accumulation processing on the corrected signal to obtain Doppler spectrum energy distribution; The radial velocity value is calculated according to the peak position of the Doppler spectrum energy distribution as the high-precision motion parameter inversion result.

[0013] Furthermore, the high-precision motion parameter inversion results are used for the joint inversion of atmospheric temperature, humidity, wind and rain profiles.

[0014] In another aspect, the present invention provides an all-weather atmospheric temperature, humidity, wind and rain profile joint detection system, comprising the following modules: The original echo acquisition module is used to obtain the original echo signal received by the radar on the detection path; Phase difference calculation module, used to extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; An interface position identification module is used to identify the interface position where the dielectric constant mutation occurs when the phase difference value exceeds a preset linear threshold range; The recursion graph analysis module is used to perform recursion graph quantitative analysis on the interface positions of dielectric constant mutations for a preset number of consecutive cycles, extract the diagonal structural characteristics of the recursion graph as the deterministic coefficient, and calculate the vertical segment distribution as the laminar flow parameter; A compensation strategy decision module is used to determine the dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and to execute the path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; The signal correction processing module is used to perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing, and output high-precision motion parameter inversion results.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. By dynamically identifying interfaces with dielectric constant mutations and quantifying atmospheric dynamic stability, the inversion accuracy of motion parameters in complex atmospheric environments is significantly improved. Unlike traditional linear accumulation models, this approach precisely locates dielectric interfaces based on phase differential mutation detection. Combining the diagonal structure and vertical segment distribution of the recursive graph to extract deterministic coefficients and laminar flow parameters, this approach constructs a quantitative assessment system for atmospheric dynamic stability. This effectively addresses the issue of nonlinear phase distortion caused by dielectric mutations and provides a physical basis for subsequent compensation strategies.

[0016] Through adaptive correction of the path integral compensation strategy, the linear characteristics of the phase accumulation are reconstructed while retaining the coherence of the original signal. The dielectric mutation detection, stability classification and phase compensation closed-loop are linked, so that the coherent accumulation algorithm can still maintain the spectral energy concentration in inhomogeneous media, and ultimately output high-precision motion parameters. Compared with existing technologies, the error of the atmospheric temperature, humidity, wind and rain profile inversion in strong gradient interface scenarios is reduced without increasing hardware costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 This is a flow chart of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection method of the present invention; Figure 2 It is a structural schematic diagram of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection system of the present invention. DETAILED DESCRIPTION

[0018] The following will provide a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0019] Example 1: Figure 1 The present invention provides an all-weather atmospheric temperature, humidity, wind and rain profile joint detection method, which includes the following steps: S1: Obtain the original echo signal received by the radar on the detection path; S2: Extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; S3: When the phase difference value exceeds the preset linear threshold range, the interface position of the dielectric constant mutation is identified; S4: Quantitative analysis of the dielectric constant mutation interface position is performed on the recursive graph for a preset number of consecutive cycles, the diagonal structural characteristics of the recursive graph are extracted as the deterministic coefficient, and the vertical segment distribution is calculated as the laminar flow parameter; S5: determining a dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and executing a path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; S6: Perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing to output high-precision motion parameter inversion results.

[0020] S1: Obtain the original echo signal received by the radar on the detection path. The specific implementation is as follows: The radar system transmits a linear frequency-modulated pulse signal through an antenna array fixed to a meteorological observation station. This signal has a carrier frequency of, for example, 35 GHz and a pulse repetition frequency of, for example, 2000 Hz. This system is the core component of the TK001THWR atmospheric temperature, humidity, wind, and rain profiler radar. Designed for simultaneous acquisition of atmospheric temperature, humidity, wind speed, wind direction, and vertical profiles of precipitation liquid water, this system is capable of all-weather operation (stable operation in cloud, rain, and fog conditions), with a target vertical resolution of 10 meters and a temporal resolution of ≤1 minute. Electromagnetic waves propagate through the atmosphere, generating backscatter when they encounter aerosol particles or cloud droplets. The radar receiver synchronously samples electromagnetic wave energy along the detection path via two orthogonal mixing channels (labeled Channel I and Channel Q). The local oscillator signal is generated by a phase-locked loop circuit and phase-locked to the transmitted signal. The receiver's intermediate frequency output is connected to a dual-channel analog-to-digital converter, which simultaneously samples two analog voltage signals at a sampling rate of, for example, 5 MHz per second, forming the initial orthogonal dual-channel echo signal data stream.

[0021] Due to differences in the physical paths of the hardware circuits, the orthogonal dual-channel echo signals exhibit time offsets. The processing unit performs a time alignment operation: first, the cross-correlation function between the Channel I and Channel Q signals is calculated to locate the offset point of the cross-correlation function's maximum value. A polynomial interpolation algorithm (for example, a cubic interpolation function constructed using four adjacent sampling points) is then applied to the delayed channel to resample the signals, aligning the rising edges of the two signals to the same sampling point index. This process generates synchronized orthogonal dual-channel echo signals, with a time synchronization error within, for example, 1% of the sampling interval (less than 2 nanoseconds for a 0.2 microsecond sampling interval).

[0022] The synchronized quadrature dual-channel echo signals contain independent real and imaginary components. The processing unit performs signal separation: the voltage value output by synchronized channel I is defined as the in-phase component of the original echo signal; the voltage value output by channel Q is defined as the quadrature component of the original echo signal. The in-phase component voltage range is, for example, -2.5 volts to +2.5 volts, while the quadrature component has the same voltage range. Both are stored in a two's complement format (e.g., a 16-bit signed integer) and stored in a random access memory buffer.

[0023] The process of combining the in-phase and quadrature components into a complex raw echo signal is achieved through numerical reconstruction: for each range gate index n, the in-phase component value In and the quadrature component value Qn are extracted. The real value In and the imaginary value jQn together form the complex number Sn (j represents the imaginary unit). This complex sequence is arranged in ascending order by range gate index in a storage array with dimensions of, for example, 2048 rows x 1 column (corresponding to 2048 range gates). Each element contains two 32-bit floating-point numbers, respectively, storing the real and imaginary values. This complex raw echo signal serves as the input data source for the subsequent phase extraction process.

[0024] S2: Extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates. The specific implementation is as follows: The complex raw echo signal is stored in the memory of the digital signal processor. This signal contains data for multiple range gates generated in the previous step. Each range gate data consists of an in-phase component and a quadrature component. For each range gate, the corresponding in-phase and quadrature component values ​​are extracted. The phase value for that range gate is calculated using a four-quadrant inverse tangent function. First, the sign combination of the in-phase and quadrature component values ​​is determined. When both the in-phase and quadrature components are positive, the phase value is equal to the inverse tangent of the quadrature component divided by the in-phase component value. When the in-phase component value is negative, the phase value is equal to the inverse tangent of the quadrature component divided by the in-phase component value plus the radian value corresponding to pi. When both the in-phase component value is positive and the quadrature component value is negative, the phase value is equal to the inverse tangent of the quadrature component divided by the in-phase component value plus the radian value corresponding to twice pi. This method ensures that the calculated phase value for each range gate ranges from zero radians to twice pi. For example, if the in-phase component value of a range gate is 0.8 and the quadrature component value is 0.6, then its phase value is approximately 0.6435 radians.

[0025] Arrange all range gates in chronological order based on their corresponding detection times, with the earliest detection gate at the front and the latest detection gate at the back. The detection time order is directly determined by the range gate number: smaller range gates correspond to earlier detection times, and larger range gates correspond to later detection times. Following this order, the phase values ​​for each range gate are arranged in sequence, forming an ordered sequence called a phase change sequence. This sequence fully records the phase variation along the detection path as it changes with detection time (or distance). For example, if there are 2048 range gates, the phase change sequence contains 2048 phase values, with the first value corresponding to range gate 1 and the last value corresponding to range gate 2048.

[0026] In the phase change sequence, the phase difference corresponding to each pair of adjacent range gates is calculated. Specifically, the phase value of the previous range gate is subtracted from the phase value of the next range gate, resulting in a phase difference. This value reflects the phase change between the two adjacent range gates. During the calculation process, if the obtained phase difference value is less than a negative pi radian, the phase difference value is increased by twice the pi radian; if the obtained phase difference value is greater than a positive pi radian, the phase difference value is reduced by twice the pi radian. This processing step ensures that each phase difference value ultimately falls between negative and positive pi radians, which is consistent with physical reality. For example, the phase value of the previous range gate is 5.5 radians, and the phase value of the next range gate is 0.1 radians. Direct subtraction results in negative 5.4 radians. Because negative 5.4 is less than negative pi (approximately negative 3.14), it is corrected to negative 5.4 plus twice pi (approximately 6.28), which is approximately 0.88 radians.

[0027] Each phase difference value, after the calculation and correction described above, is stored sequentially in a dedicated array, in the order of its corresponding adjacent range gate pairs. This array is called the phase difference sequence array. The length of this array is one less than the total number of range gates, because each pair of adjacent range gates produces a phase difference value. Each element in the array stores the phase difference between two adjacent range gates at the corresponding position. For example, the first element of the array stores the phase difference between range gates 1 and 2, and the last element stores the phase difference between range gates 2047 and 2048. This array is stored in memory in floating-point format for use in subsequent processing steps.

[0028] S3: When the phase difference value exceeds the preset linear threshold range, the dielectric constant mutation interface position is identified, which is specifically implemented as follows: The phase difference sequence array is stored in a designated memory area of ​​the digital signal processor. This array contains the phase difference values ​​of adjacent range gates calculated in the previous step. The array length is one less than the total number of range gates. For example, if the total number of range gates is 2048, the array contains 2047 phase difference values. The processing unit starts at the beginning of the array and reads each phase difference value sequentially in index order. The read operation is performed using direct memory access, reading one 32-bit floating-point number at a time until all elements in the array are traversed. During the traversal process, the adjacent range gate number corresponding to the currently processed phase difference value is recorded. This number is calculated by adding 1 to the array index to obtain the previous range gate number and adding 2 to the array index to obtain the next range gate number.

[0029] The preset linear threshold range is defined by two boundary values: a lower limit of -0.3 radians and an upper limit of +0.3 radians. This threshold range is set based on the atmospheric dielectric constant gradient model, and the phase difference value is generally within this range when the atmospheric parameters change continuously. The threshold parameters are stored in non-volatile memory and loaded into the register bank when the system starts. The processing unit performs a double-boundary comparison between the currently traversed phase difference value and the lower and upper limits of the preset linear threshold range: first, it determines whether the phase difference value is less than -0.3 radians, and second, it determines whether it is greater than +0.3 radians. The comparison operation is implemented using a floating-point comparator hardware circuit, which outputs a binary status flag.

[0030] When the double-boundary comparison results show a phase difference less than -0.3 radians or greater than +0.3 radians, a sudden change in the dielectric constant is determined. The location of the sudden change is recorded: the two range gate numbers corresponding to the current phase difference are obtained, and the range gate with the smaller number is used as the marker location. The range gate number is used to query a pre-stored spatial coordinate mapping table. This table, generated during radar initialization, stores the three-dimensional spatial coordinates of each range gate. The spatial coordinates are based on the northeast celestial coordinate system, with the origin at the radar antenna phase center. The east coordinate is calculated by multiplying the sine of the radar azimuth angle by the slant range, the north coordinate by multiplying the cosine of the radar azimuth angle by the slant range, and the celestial coordinate by multiplying the sine of the radar elevation angle by the slant range. For example, for range gate number 100, with a radar azimuth of 30 degrees, an elevation of 5 degrees, and a slant range of 3000 meters, its spatial coordinates are 1500 meters east, 2598 meters north, and 261 meters celestial.

[0031] The acquired spatial coordinates are projected onto the radar detection path. The radar detection path is defined as a unit vector pointing from the radar position to the current scanning direction. The projection operation is implemented using a vector dot product calculation: a position vector is constructed with the radar position as the origin and the target's spatial coordinates as the endpoint. The dot product of this position vector and the radar beam pointing unit vector is calculated. The resulting scalar value is the projected slant distance of the target point along the detection path. This projected slant distance is accurate to 1 meter resolution. The calculated slant distance value is finally marked as the location of the dielectric constant abrupt change interface and stored in the abrupt change location record array, along with the corresponding millisecond detection timestamp. For example, the spatial coordinates (1500, 2598, 261) projected slant distance to a beam direction with an azimuth angle of 30 degrees and an elevation angle of 5 degrees is 3000 meters.

[0032] S4: Perform quantitative analysis of the recursive graph of the dielectric constant mutation interface position for a preset number of consecutive cycles, extract the diagonal structural characteristics of the recursive graph as the deterministic coefficient, and calculate the vertical segment distribution as the laminar flow parameter. The specific implementation is as follows: The radar system's memory unit retrieves records of the locations of the dielectric constant abrupt change interface for a predetermined number of consecutive cycles. This number is determined based on the typical timescale of atmospheric motion, for example, 10 consecutive complete radar scan cycles. Each cycle corresponds to a complete sweep of the radar beam along the detection path. The cycle duration is determined by the radar pulse repetition frequency (PRF). For example, at a 2000 Hz PRF, the cycle is 0.0005 seconds multiplied by the number of range gates. All dielectric constant abrupt change interface locations detected within each cycle are sorted chronologically, with the locations detected in the first cycle first and those detected in the last cycle last, forming a time series of interface locations. Each location in the sequence contains three-dimensional spatial coordinate information (easting, northing, and celestial coordinates) and a time stamp accurate to milliseconds. For example, the total length of the sequence might be 150 locations, with location 1 corresponding to the earliest detection time and location 150 corresponding to the latest detection time.

[0033] The recursive threshold is used to determine the degree of position similarity. It is set by first calculating the straight-line distances between all pairs of position points in the interface position time series and finding the maximum distance. Then, 5% of this maximum distance is used as the recursive threshold. For example, if the maximum distance is 500 meters, the recursive threshold is set to 25 meters. This threshold parameter is stored in a floating-point register. Straight-line distance calculation uses the three-dimensional distance formula: take the square of the difference in the easting coordinates of the two points, add the square of the difference in the northing coordinates, and add the square of the difference in the celestial coordinates. The square root of the sum is taken to obtain the distance. The calculation is performed to one decimal place.

[0034] A square two-dimensional grid is constructed to generate the recurrence graph. The horizontal and vertical axes of the grid represent the index numbers of the interface position time series, ranging from the first to the last position in the sequence. For any two position points at different time points in the sequence (for example, position points 50 and 100), the straight-line distance between them is calculated. When this distance value is less than a recurrence threshold (for example, 25 meters), the intersection of the corresponding grid row and column coordinates is marked as a recurrence point (usually represented by a value of 1); otherwise, it is marked as a non-recurrence point (usually represented by a value of 0). After traversing all possible combinations of position point pairs, a complete recurrence graph matrix is ​​formed. For example, 150 position points will generate a grid matrix with 150 rows and 150 columns.

[0035] Identify line segments formed by consecutive recursive points in a recursive graph along a diagonal direction. The diagonal direction is defined as the direction parallel to the main diagonal of the grid, i.e., the direction with a constant difference between the row and column numbers. Scan along each line parallel to the main diagonal: Whenever at least two recursive points appear consecutively, this is counted as a line segment. Calculate the length distribution of all such line segments: set the length statistics interval to short segments (2-3 consecutive points), medium segments (4-6 consecutive points), long segments (7-10 consecutive points), and very long segments (11 or more consecutive points). Calculate the proportion of line segments in each interval to the total number of segments. Generate a eigenvector containing the four proportions as a deterministic coefficient. For example, the statistical result might be that short segments account for 15%, medium segments account for 35%, long segments account for 45%, and very long segments account for 5%.

[0036] Detect the line segments formed by consecutive recursive points in the vertical direction of the recursive graph. The vertical direction is defined as the change in the direction of increasing column numbers while the row number remains unchanged. For each fixed row number (for example, row 80), scan along the column number direction: when at least two recursive points appear consecutively in this row, it is recorded as a vertical line segment. Count the length distribution of all vertical line segments: use the same four length intervals as the diagonal lines; calculate the proportion of line segments in each interval; and form a feature vector containing four proportion values ​​as the laminar flow parameter. For example, the result may be a 25% proportion of short vertical segments, a 40% proportion of medium vertical segments, a 30% proportion of long vertical segments, and a 5% proportion of very long vertical segments. The final output of the deterministic coefficient and laminar flow parameter vector is stored in the register group for subsequent stability level calculation.

[0037] S5: Determine the dynamic stability level based on the deterministic coefficient and the laminar flow parameter. When the dynamic stability level is lower than the preset stability threshold, execute the path integral compensation strategy. The specific implementation is as follows: The determinism coefficient and laminarity parameter are stored in the processing unit registers. Both are eigenvectors containing four proportional values, generated by the recursive graph analysis in the previous step. These two eigenvectors are input into the preset stability mapping model. The model is constructed as follows: First, the four proportional values ​​of the determinism coefficient eigenvector are weighted and summed according to the weight coefficients 0.4, 0.3, 0.2, and 0.1 to obtain a composite value (for example, input [0.15, 0.35, 0.45, 0.05]); The result is 0.4 × 0.15 + 0.3 × 0.35 + 0.2 × 0.45 + 0.1 × 0.05 = 0.32. The same method is used to calculate the composite value of the laminar flow parameters (for example, [0.25, 0.40, 0.30, 0.05] yields 0.35). A two-dimensional plane, with the composite value of the deterministic coefficient on the horizontal axis and the composite value of the laminar flow parameters on the vertical axis, is divided into five dynamic stability levels (Level 1 is the most stable, Level 5 is the least stable). The boundaries of these zones are determined based on meteorological experimental data. When a eigenvector is input, the composite value is calculated and the coordinate point is located (for example, 0.32, 0.35). The corresponding level value of the zone in which this point is located is then found (for example, if it is in the Level 3 zone, the value 3 is output).

[0038] Read the preset stability threshold, which is determined by the radar operating wavelength (e.g., 5 cm) and dielectric uniformity. Dielectric uniformity is calculated by calculating the ratio of the standard deviation of the dielectric constant to the mean over the last 10 minutes (e.g., a ratio of 0.25 corresponds to uniformity level 3). In the pre-stored data table, the threshold for uniformity level 3 at a 5 cm wavelength is set to 2.5. The current threshold is updated every minute and stored in a dedicated register.

[0039] Compare the dynamic stability level value to the preset stability threshold: Use a digital comparator to determine the numerical relationship. If the dynamic stability level value is less than the preset stability threshold (for example, 3.0 < 2.5 is not true), a logically false signal is output; otherwise, a logically true signal is output. The comparison result is recorded in the status flag.

[0040] When the status flag indicates a logical false signal (i.e., the dynamic stability level is less than the threshold), a 32-bit control instruction is generated. The first 16 bits of this instruction are fixed to the strategy code 64521, and the last 16 bits are the current scan cycle number (e.g., 258). This instruction is transmitted to the compensation unit via the bus, simultaneously generating an interrupt request with a priority level of 2.

[0041] Initialize the compensation parameter matrix according to the instructions. The number of rows in the matrix is ​​equal to the total number of range gates (for example, 2048 rows), and it contains two columns of data. The first column, the range gate compensation vector, is initialized to 2048 1.0 values, and the second column, the time step compensation vector, is also initialized to 2048 1.0 values. The following operations are performed sequentially in memory: first, the current cycle number is written; then, 2048 initial range gate compensation values ​​are written continuously; and finally, 2048 initial time step compensation values ​​are written continuously. This results in a 2048-row x 2-column matrix for subsequent use.

[0042] S6: Perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing, and output high-precision motion parameter inversion results. The specific implementation is as follows: The path integral compensation parameter matrix is ​​read from a dual-port random access memory (DPRM). This matrix contains 2048 rows and two columns. The first column stores the range gate compensation vector, which consists of 2048 floating-point values; the second column stores the time step compensation vector, which also consists of 2048 floating-point values. This fetch is performed using direct memory access: first, the matrix base address is located, and the 2048 range gate compensation values ​​are read consecutively starting at an offset of 0 bytes from the base address. Then, the 2048 time step compensation values ​​are read consecutively starting at an offset of 16384 bytes from the base address. For example, the 100th element of the range gate compensation vector corresponds to range gate number 100 and has an initial value of 1.0.

[0043] Perform an element-by-element multiplication of the range gate compensation vector and the phase component of the original echo signal. The original echo signal, stored in the digital signal processor's buffer, consists of 2048 complex numbers, each with its imaginary part representing the phase component. The calculation proceeds as follows: For each range gate index (1 to 2048), multiply the value at the corresponding position in the range gate compensation vector by the phase component value for that range gate in the original echo signal. For example, for range gate 50, a compensation value of 0.98 multiplied by a phase component of 0.75 radians yields 0.735 radians. This result forms the primary compensation signal, which is stored in a newly allocated memory area and has the same data structure as the original echo signal.

[0044] The time step compensation vector is element-wise multiplied by the phase component of the primary compensation signal. The values ​​in the time step compensation vector are arranged in chronological order. For each range gate index (1 to 2048), the value at the corresponding position in the time step compensation vector is multiplied by the phase component value of the primary compensation signal for that range gate. For example, for range gate 150, a time step compensation of 1.02 is multiplied by the primary phase component of 0.82 radians, resulting in a value of 0.8364 radians. This result is a corrected signal, which retains the real component of the original echo signal unchanged and only modifies the phase component. This signal is stored in the digital signal processor's output buffer.

[0045] The corrected signal is subjected to a discrete Fourier transform (DFT) coherent accumulation process. This process involves three steps: first, the corrected signal for 32 consecutive pulse periods is arranged in time order to form a 32 × 2048 complex matrix. Then, for each range gate index (from 1 to 2048), a discrete Fourier transform (DFT) is performed on the 32 complex numbers along the time dimension. Finally, the squared modulus of the transform result is calculated to obtain the Doppler spectrum energy distribution for each range gate. For example, the spectrum output by range gate 1000 contains 32 energy values, and the frequency resolution is determined by the pulse repetition frequency (a 2000 Hz repetition frequency corresponds to a 62.5 Hz resolution).

[0046] The radial velocity value is calculated based on the peak position of the Doppler spectrum energy distribution. First, the index of the maximum energy value in each range gate spectrum is located. For example, the peak value appears at index 15 for range gate 800. The velocity calculation formula is: radial velocity equals the peak position index multiplied by the wavelength multiplied by the pulse repetition frequency divided by the number of discrete Fourier transform points. For example, with a wavelength of 0.1 meters, a pulse repetition frequency of 2000 Hz, and 32 transform points, the velocity corresponding to index 15 is 15 × 0.1 × 2000 / 32 = 93.75 meters per second. The calculation result is stored in the result register group as a high-precision motion parameter inversion result. This high-precision motion parameter inversion result is one of the key inputs for achieving single-device synchronous inversion of atmospheric temperature, humidity, wind direction, and rainfall intensity vertical profiles. Its accuracy directly affects the quality of the final profile product.

[0047] The high-precision motion parameter inversion results are transmitted to the atmospheric parameter solution unit for joint inversion of the vertical distribution of atmospheric temperature, humidity, wind direction, and rainfall intensity. Specifically, the vertical distribution of atmospheric temperature profiles, relative humidity profiles, liquid water profiles, wind speed profiles, wind direction profiles, raindrop spectrum, and rainfall rate are analyzed. The solution process is based on a set of meteorological fluid dynamics equations, using the radial velocity field as a key input parameter. The atmospheric state estimate is optimized using a variational assimilation algorithm. For example, in boundary layer wind field inversion, radial velocity data is used to constrain the horizontal momentum term in the three-dimensional wind field reconstruction equation.

[0048] Example 2: Figure 2 The structural diagram of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection system of the present invention is given. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection system includes the following modules: The original echo acquisition module is used to obtain the original echo signal received by the radar on the detection path; Phase difference calculation module, used to extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; An interface position identification module is used to identify the interface position where the dielectric constant mutation occurs when the phase difference value exceeds a preset linear threshold range; The recursion graph analysis module is used to perform recursion graph quantitative analysis on the interface positions of dielectric constant mutations for a preset number of consecutive cycles, extract the diagonal structural characteristics of the recursion graph as the deterministic coefficient, and calculate the vertical segment distribution as the laminar flow parameter; A compensation strategy decision module is used to determine the dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and to execute the path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; The signal correction processing module is used to perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing, and output high-precision motion parameter inversion results.

[0049] The calculations involved in the embodiments are all dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to actual conditions.

[0050] The above embodiments may be implemented in whole or in part through software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments may be implemented in whole or in part in the form of a computer program product.

[0051] Those skilled in the art will appreciate that the modules and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application of the technical solution and the invention constraints. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0052] In addition, each functional module in each embodiment of the present application may be integrated into one processing module, or each module may exist physically separately, or two or more modules may be integrated into one module.

[0053] In the several embodiments provided in this application, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are merely schematic. For example, the division of the modules is only a logical function division. In actual implementation, there may be other division methods, such as multiple modules or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or modules, which can be electrical, mechanical or other forms.

[0054] The above description is merely a specific embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of this application. Therefore, the scope of protection of this application should be based on the scope of protection of the claims.

[0055] Finally: The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. All-weather atmospheric temperature, humidity, wind and rain profile joint detection method, characterized by: The steps include: S1: Obtain the original echo signal received by the radar on the detection path; S2: Extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; S3: When the phase difference value exceeds the preset linear threshold range, the interface position of the dielectric constant mutation is identified; S4: Quantitative analysis of the dielectric constant mutation interface position is performed on the recursive graph for a preset number of consecutive cycles, the diagonal structural characteristics of the recursive graph are extracted as the deterministic coefficient, and the vertical segment distribution is calculated as the laminar flow parameter; S5: determining a dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and executing a path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; S6: Perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing to output high-precision motion parameter inversion results.

2. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 1 is characterized in that: Obtain the original echo signal received by the radar on the detection path, including: The orthogonal dual-channel echo signals on the detection path are collected by the radar receiver; Performing time alignment processing on the orthogonal dual-channel echo signals to generate synchronized orthogonal dual-channel echo signals; Separate the real part and the imaginary part from the synchronous orthogonal dual-channel echo signal to form the in-phase component and the orthogonal component of the original echo signal; The in-phase component and the quadrature component are combined into a complex original echo signal.

3. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 2 is characterized in that: Extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates, including: Perform arc tangent operation on the in-phase component and quadrature component of the original echo signal to obtain the phase value of each range gate; Arrange the phase values ​​of all range gates in the order of detection time to form a phase change sequence; Calculate the difference between the phase values ​​of two adjacent range gates in the phase change sequence to obtain a phase difference score; Store the phase difference values ​​in the phase difference sequence array.

4. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 3 is characterized in that: When the phase difference value exceeds the preset linear threshold range, the dielectric constant mutation interface position is identified, including: Traverse each phase difference value in the phase difference sequence array; Performing a double-boundary comparison between the current phase difference value and the upper and lower limits of a preset linear threshold range; When the phase difference value is greater than the upper limit of the preset linear threshold range or less than the lower limit of the preset linear threshold range, the spatial coordinates of the corresponding range gate are recorded; The projection position of the spatial coordinate on the radar detection path is marked as the position of the dielectric constant mutation interface.

5. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 4 is characterized in that: The recursive graph is used to quantitatively analyze the interface position of the dielectric constant mutation for a preset number of consecutive cycles. The diagonal structural characteristics of the recursive graph are extracted as the deterministic coefficient, and the vertical segment distribution is calculated as the laminar flow parameter, including: Arranging the dielectric constant mutation interface positions of a preset number of consecutive cycles in chronological order as an interface position time series; Set the recursive threshold and calculate the Euclidean distance between the interface positions corresponding to every two time points in the interface position time series; When the Euclidean distance is less than the recursion threshold, the recursion point is marked on the corresponding coordinates of the two-dimensional grid to generate a recursion graph; Identify the line segments formed by continuous recursive points in the diagonal direction of the recursive graph, and calculate the characteristic value of the line segment length distribution as the deterministic coefficient; The line segments formed by the vertical continuous recursive points in the recursive graph are detected, and the characteristic values ​​of the line segment length distribution are statistically analyzed as laminar flow parameters.

6. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 5 is characterized in that: The dynamic stability level is determined based on the deterministic coefficient and the laminar flow parameter. When the dynamic stability level is lower than the preset stability threshold, a path integral compensation strategy is executed, including: Input the deterministic coefficient and laminar flow parameters into the preset stability mapping model and output the dynamic stability level value; Read the preset stability threshold based on radar wavelength and medium uniformity calibration; Compare the dynamic stability level value with the preset stability threshold; When the dynamic stability level value is less than the preset stability threshold, the path integral compensation strategy is activated to generate instructions; Generate instructions based on the path integral compensation strategy to initialize the path integral compensation parameter matrix.

7. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 6, characterized in that: The preset stability mapping model is implemented in the following way: Establish a two-dimensional feature space of deterministic coefficients and laminar flow parameters; Divide the dynamic stability level value into contour areas in the two-dimensional feature space; The dynamic stability level value corresponding to the input deterministic coefficient and laminar flow parameter is calculated by linear interpolation.

8. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 6, characterized in that: Based on the path integral compensation strategy, the original echo signal is phase corrected to generate the corrected signal and perform coherent accumulation processing to output high-precision motion parameter inversion results, including: Extracting the range gate compensation vector and the time step compensation vector from the path integral compensation parameter matrix; Performing a Hadamard product operation on the range gate compensation vector and the phase component of the original echo signal to generate a primary compensation signal; Performing a Hadamard product operation on the time step compensation amount vector and the phase component of the primary compensation signal to generate a corrected signal; Performing discrete Fourier transform coherent accumulation processing on the corrected signal to obtain Doppler spectrum energy distribution; The radial velocity value is calculated according to the peak position of the Doppler spectrum energy distribution as the high-precision motion parameter inversion result.

9. The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to claim 8, characterized in that: in, The high-precision motion parameter inversion results are used for the joint inversion of atmospheric temperature, humidity, wind and rain profiles.

10. An all-weather atmospheric temperature, humidity, wind and rain profile joint detection system, used to implement the all-weather atmospheric temperature, humidity, wind and rain profile joint detection method according to any one of claims 1 to 9, characterized in that: Includes the following modules: The original echo acquisition module is used to obtain the original echo signal received by the radar on the detection path; Phase difference calculation module, used to extract the phase change sequence from the original echo signal and calculate the phase difference value of adjacent range gates; An interface position identification module is used to identify the interface position where the dielectric constant mutation occurs when the phase difference value exceeds a preset linear threshold range; The recursion graph analysis module is used to perform recursion graph quantitative analysis on the interface positions of dielectric constant mutations for a preset number of consecutive cycles, extract the diagonal structural characteristics of the recursion graph as the deterministic coefficient, and calculate the vertical segment distribution as the laminar flow parameter; A compensation strategy decision module is used to determine the dynamic stability level according to the deterministic coefficient and the laminar flow parameter, and to execute the path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; The signal correction processing module is used to perform phase correction on the original echo signal based on the path integral compensation strategy, generate the corrected signal and perform coherent accumulation processing, and output high-precision motion parameter inversion results.

Images