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

Abstract

The application discloses an all-weather atmospheric temperature, humidity, wind and rain profile combined detection method and system, and particularly relates to the technical field of weather radar signal processing, and is used for solving the problem of invalidation of movement parameter inversion caused by nonlinear distortion of phase at the dielectric mutation interface of the existing coherent Doppler radar; the original echo signal of the radar is acquired; the phase change sequence is extracted and the phase difference value of adjacent distance gates is calculated; the dielectric mutation interface position is identified when the difference value exceeds the linear threshold range; the recursive graph analysis of the continuous period interface position is carried out, the diagonal line feature is extracted as the determinacy coefficient, and the vertical segment distribution is extracted as the laminar parameter; the dynamic stability level is determined based on the two, and the path integral compensation is executed when the level is lower than the threshold value; the original signal phase is corrected through the compensation strategy, high-precision movement parameters are output through coherent accumulation, the dynamic compensation of the phase distortion in the non-uniform medium is realized, and the robustness of the atmospheric parameter inversion is ensured.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of meteorological radar signal processing, more particularly, the present application relates to all-weather atmospheric temperature, humidity, wind and rain profile joint detection method and system. BACKGROUND

[0002] In the remote sensing detection of atmospheric medium by millimeter wave radar, the coherent Doppler processing mechanism is usually used to obtain target motion information. The prior art uses the phase change of the signal to invert the radial velocity by transmitting a pulse signal and receiving the backscattered wave. When the radar wave propagates in a homogeneous medium, the phase change shows a linear cumulative characteristic, and high-precision velocity measurement can be realized based on the standard coherent accumulation algorithm. Such a method has been widely used in the detection of aerosol and particle motion in the fields of meteorology, environmental monitoring, etc.

[0003] However, when the radar beam passes through a sudden medium interface with a significant gradient of dielectric constant (such as an alternating region of meteorological targets and non-meteorological media), the wave front is distorted due to non-uniform phase shift, resulting in deviation of the phase accumulation characteristics of the echo signal from the ideal linear model. This phenomenon destroys the precondition of coherent processing, causing the spectral energy output by the traditional accumulation algorithm to be dispersed and the signal-to-noise ratio to deteriorate, ultimately resulting in systematic jump errors in the motion parameter inversion results. SUMMARY

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

[0005] To achieve the above-mentioned purpose, the present application provides the following technical solutions:

[0006] The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method comprises the following steps:

[0007] S1: obtaining the original echo signal received by the radar on the detection path;

[0008] S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates;

[0009] S3: identifying the dielectric constant sudden interface position when the phase difference value exceeds the preset linear threshold range;

[0010] S4: performing recursive graph quantization analysis on the dielectric constant sudden interface positions for a continuous preset number of periods, extracting the recursive graph diagonal line structure feature as the determinacy coefficient, and calculating the vertical segment distribution as the laminar flow parameter;

[0011] S5: determining the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and executing the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold;

[0012] S6: phase correction is performed on the original echo signal based on the path integral compensation strategy, a corrected signal is generated, coherent accumulation processing is performed, and a high-precision motion parameter inversion result is output.

[0013] Further, the original echo signal received by the radar on the detection path is obtained, including:

[0014] The orthogonal double-channel echo signal on the detection path is collected by the radar receiver;

[0015] The time alignment processing is performed on the orthogonal double-channel echo signal to generate a synchronous orthogonal double-channel echo signal;

[0016] The real part and the imaginary part are separated from the synchronous orthogonal double-channel echo signal to form the in-phase component and the quadrature component of the original echo signal;

[0017] The in-phase component and the quadrature component are combined into a complex form original echo signal.

[0018] Further, the phase change sequence is extracted from the original echo signal, and the phase difference value of adjacent range gates is calculated, including:

[0019] The in-phase component and the quadrature component of the original echo signal are subjected to arctangent operation to obtain the phase value of each range gate;

[0020] The phase values of all range gates are arranged in the order of detection time to form a phase change sequence;

[0021] The phase values of adjacent two range gates in the phase change sequence are subtracted to calculate the phase difference value;

[0022] The phase difference value is stored in a phase difference sequence array.

[0023] Further, when the phase difference value exceeds a preset linear threshold range, the dielectric constant mutation interface position is identified, including:

[0024] Each phase difference value in the phase difference sequence array is traversed;

[0025] The current phase difference value is compared with the upper limit value and the lower limit value of the preset linear threshold range;

[0026] When the phase difference value is greater than the upper limit value of the preset linear threshold range or less than the lower limit value of the preset linear threshold range, the spatial coordinates of the corresponding range gate are recorded;

[0027] The projection position of the spatial coordinates on the radar detection path is marked as the dielectric constant mutation interface position.

[0028] Further, the dielectric constant mutation interface positions of a continuous preset number of periods are subjected to recurrence plot quantitative analysis, the diagonal line structure features of the recurrence plot are extracted as the determinacy coefficient, and the vertical segment distribution is calculated as the laminar flow parameter, including:

[0029] The dielectric constant mutation interface positions of a continuous preset number of periods are arranged in time sequence as an interface position time sequence;

[0030] A recurrence threshold is set, and the Euclidean distance between the interface positions corresponding to each two time points in the interface position time sequence is calculated;

[0031] When the Euclidean distance is less than the recurrence threshold, a recurrence point is marked in the corresponding coordinates of the two-dimensional grid, and a recurrence plot is generated;

[0032] The line segments formed by the diagonal line direction continuous recurrence points in the recurrence plot are identified, and the line segment length distribution characteristic value is calculated as the determinacy coefficient;

[0033] The line segments formed by the vertical direction continuous recurrence points in the recurrence plot are detected, and the line segment length distribution characteristic value is calculated as the laminar flow parameter.

[0034] Further, the dynamic stability level is determined according to the determinacy coefficient and the laminar flow parameter, and when the dynamic stability level is lower than a preset stability threshold, a path integral compensation strategy is executed, including:

[0035] The determinacy coefficient and the laminar flow parameter are input into a preset stability mapping model, and a dynamic stability level value is output;

[0036] A preset stability threshold calibrated based on the radar wavelength and the medium uniformity is read;

[0037] The size relationship between the dynamic stability level value and the preset stability threshold is compared;

[0038] When the dynamic stability level value is less than the preset stability threshold, a path integral compensation strategy generation instruction is activated;

[0039] The path integral compensation parameter matrix is initialized according to the path integral compensation strategy generation instruction.

[0040] Further, the preset stability mapping model is implemented by the following method:

[0041] A two-dimensional feature space of the determinacy coefficient and the laminar flow parameter is established;

[0042] The contour area of the dynamic stability level value is divided in the two-dimensional feature space;

[0043] The dynamic stability level value corresponding to the input determinacy coefficient and laminar flow parameter is calculated by linear interpolation.

[0044] Further, the original echo signal is phase-corrected based on a path integral compensation strategy, a corrected signal is generated and coherent accumulation processing is performed, and a high-precision motion parameter inversion result is output, including:

[0045] A distance gate compensation quantity vector and a time step compensation quantity vector are extracted from the path integral compensation parameter matrix;

[0046] The distance gate compensation quantity vector and the phase component of the original echo signal are subjected to Hadamard product operation to generate a primary compensation signal;

[0047] The time step compensation quantity vector and the phase component of the primary compensation signal are subjected to Hadamard product operation to generate a corrected signal;

[0048] Discrete Fourier transform coherent accumulation processing is performed on the corrected signal to obtain a Doppler spectrum energy distribution;

[0049] A radial velocity value is calculated according to a peak position of the Doppler spectrum energy distribution, as a high-precision motion parameter inversion result.

[0050] Further, the high-precision motion parameter inversion result is used for joint inversion of atmospheric temperature, humidity, wind and rain profiles.

[0051] On the other hand, the present application provides an all-weather atmospheric temperature, humidity, wind and rain profile joint detection system, comprising the following modules:

[0052] An original echo acquisition module is configured to acquire an original echo signal received by a radar on a detection path;

[0053] A phase difference calculation module is configured to extract a phase change sequence from the original echo signal and calculate a phase difference value of adjacent distance gates;

[0054] An interface position identification module is configured to identify a dielectric constant mutation interface position when the phase difference value exceeds a preset linear threshold range;

[0055] A recurrence plot analysis module is configured to perform recurrence plot quantization analysis on dielectric constant mutation interface positions in a continuous preset period of time, extract a recurrence plot diagonal line structure feature as a determinacy coefficient, and calculate a vertical segment distribution as a laminarity parameter;

[0056] A compensation strategy decision module is configured to determine a dynamic stability level according to the determinacy coefficient and the laminarity parameter, and execute a path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold;

[0057] A signal correction processing module is configured to perform phase correction on the original echo signal based on the path integral compensation strategy, generate a corrected signal and perform coherent accumulation processing, and output a high-precision motion parameter inversion result.

[0058] Compared with the prior art, the present application has the following beneficial effects:

[0059] 1. By dynamically identifying the dielectric constant mutation interface and quantifying the atmospheric dynamic stability, the inversion accuracy of the motion parameters in a complex atmospheric environment is significantly improved. Unlike the traditional linear cumulative model, the present scheme accurately locates the dielectric interface position based on phase difference mutation detection, extracts the deterministic coefficient and laminar flow parameter by combining the recursive graph diagonal structure and vertical segment distribution, and constructs a quantitative evaluation system of atmospheric dynamic stability, effectively solving the problem of phase nonlinear distortion caused by dielectric mutation, and providing physical basis for subsequent compensation strategy.

[0060] Through the adaptive correction of the path integral compensation strategy, the linear characteristics of phase accumulation are reconstructed while preserving the coherence of the original signal, the dielectric mutation detection, stability grading and phase compensation are closed-loop linked, the coherent accumulation algorithm can still maintain the spectral energy concentration in non-uniform medium, and finally the high-precision motion parameters are output. Compared with the prior art, the error of atmospheric temperature, humidity, wind and rain profile inversion in the strong gradient interface scene is reduced, and the hardware cost is not increased. BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 The flowchart of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection method of the present application;

[0062] Figure 2 The structural schematic diagram of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection system of the present application. DETAILED DESCRIPTION

[0063] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0064] Embodiment 1: Figure 1 The all-weather atmospheric temperature, humidity, wind and rain profile joint detection method of the present application is given, which includes the following steps:

[0065] S1: obtaining the original echo signal received by the radar on the detection path;

[0066] S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates;

[0067] S3: identifying the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range;

[0068] S4: Recurrence plot quantification analysis is performed on the dielectric constant mutation interface positions of a continuous preset number of periods, the diagonal line structure features of the recurrence plot are extracted as the determinacy coefficient, and the vertical segment distribution is calculated as the laminarity parameter;

[0069] S5: The dynamic stability level is determined according to the determinacy coefficient and the laminarity parameter, and the path integral compensation strategy is executed when the dynamic stability level is lower than a preset stability threshold;

[0070] S6: The original echo signal is phase-corrected based on the path integral compensation strategy, a corrected signal is generated and coherent accumulation processing is performed, and a high-precision motion parameter inversion result is output.

[0071] S1: Obtain the original echo signal received by the radar on the detection path, which is specifically implemented as:

[0072] The radar system transmits a linear frequency modulation pulse signal through an antenna array fixed on a meteorological observation station, the carrier frequency of the signal is, for example, 35 gigahertz, and the pulse repetition frequency is, for example, 2000 hertz. The system is the core component of a TK001THWR type atmospheric temperature, humidity, wind, rain profile detection radar, which is designed to realize synchronous acquisition of atmospheric temperature, humidity, wind speed, wind direction and precipitation liquid water vertical profile by a single device, and has all-weather (stable operation under cloud, rain, fog conditions) working ability, target vertical resolution up to 10 meters, and time resolution up to minutes (≤1 minute). When electromagnetic waves encounter aerosol particles or cloud droplets during propagation in the atmosphere, backscattering occurs. The radar receiver synchronously collects electromagnetic wave energy on the detection path through two orthogonal mixing channels (labeled as channel I and channel Q, respectively), wherein the local oscillator signal is generated by a phase-locked loop circuit and maintains phase locking with the transmitted signal. The intermediate frequency output end of the receiver is connected to a dual-channel analog-to-digital converter to synchronously collect two analog voltage signals at a sampling rate of, for example, 5 megahertz per second, forming an initial orthogonal dual-channel echo signal data stream.

[0073] Due to the physical path difference of the hardware circuit, the orthogonal dual-channel echo signal has a time offset. The processing unit performs a time alignment operation: first, the cross-correlation function of the channel I and channel Q signals is calculated to locate the offset point number of the maximum value of the cross-correlation function; then, a polynomial interpolation algorithm (for example, a cubic interpolation function is constructed using four adjacent sampling points) is applied to the delayed channel for signal resampling, so that the pulse rising edges of the two signals are aligned to the same sampling point index. The synchronous orthogonal dual-channel echo signal is generated after this processing, and the time synchronization error is controlled within, for example, 1% of the sampling interval (corresponding to an error of less than 2 nanoseconds when the sampling interval is 0.2 microseconds).

[0074] The in-phase component is defined as the voltage output of the I channel, and the quadrature component is defined as the voltage output of the Q channel. The in-phase component voltage range is, for example, -2.5 volts to +2.5 volts, and the quadrature component has the same voltage range. Both the in-phase component and the quadrature component are stored in binary-complement format (e.g., 16-bit signed integer) and are stored in a random access memory buffer.

[0075] The in-phase component and the quadrature component are combined to form a complex-valued raw echo signal. For each range gate index n, the in-phase component value In and the quadrature component value Qn are extracted, and together they form a complex number Sn (j represents the imaginary unit). The complex number sequence is arranged in ascending order of the range gate number in a storage array, and the array size is, for example, 2048 rows x 1 column (corresponding to 2048 range gates). Each element contains two 32-bit floating-point numbers to store the real part and the imaginary part values, respectively. The complex-valued raw echo signal serves as the input data source for the subsequent phase extraction process.

[0076] S2: Extracting a phase change sequence from the raw echo signal and calculating the phase difference value of adjacent range gates, which is implemented as follows:

[0077] The complex-valued raw echo signal is stored in the memory of a digital signal processor, and the signal contains data of multiple range gates generated by the previous step. For each range gate, its corresponding in-phase component value and quadrature component value are extracted. The phase value of the range gate is calculated using the four-quadrant arctangent function. First, the positive and negative sign combination relationship of the in-phase component value and the quadrature component value is determined. When the in-phase component value is positive and the quadrature component value is positive, the phase value is equal to the arctangent value of the quotient of the quadrature component value divided by the in-phase component value. When the in-phase component value is negative, the phase value is equal to the arctangent value of the quotient of the quadrature component value divided by the in-phase component value plus the radian value corresponding to the constant pi. When the in-phase component value is positive and the quadrature component value is negative, the phase value is equal to the arctangent value of the quotient of the quadrature component value divided by the in-phase component value plus twice the radian value corresponding to the constant pi. In this way, the range of the phase value of each range gate calculated is ensured to be between zero radians and twice the constant pi radians. For example, if the in-phase component value of a certain range gate is 0.8 and the quadrature component value is 0.6, then the phase value is approximately 0.6435 radians.

[0078] All distance gates are arranged in order of their corresponding detection time, with the earliest detected distance gate at the front and the latest detected distance gate at the back. The detection time order is directly determined by the distance gate number, with a smaller distance gate number corresponding to an earlier detection time and a larger distance gate number corresponding to a later detection time. In this order, the phase values of each distance gate are arranged in sequence to form an ordered sequence, referred to as a phase change sequence. This sequence completely records the changes in phase along the detection path over time (or distance). For example, if there are a total of 2048 distance gates, the phase change sequence contains 2048 phase values, with the first value corresponding to distance gate 1 and the last value corresponding to distance gate 2048.

[0079] In the phase change sequence, the phase values corresponding to each pair of adjacent distance gates are calculated by subtraction. Specifically, the phase value of the latter distance gate is subtracted from the phase value of the former distance gate to obtain a result referred to as a phase difference value. This value reflects the amount of phase change between the two adjacent distance gates. During the calculation process, if the obtained phase difference value is less than negative circular constant radians, the phase difference value is increased by twice the circular constant radians; if the obtained phase difference value is greater than positive circular constant radians, the phase difference value is decreased by twice the circular constant radians. This processing step ensures that each phase difference value ultimately falls between negative circular constant radians and positive circular constant radians, in line with physical reality. For example, if the phase value of the former distance gate is 5.5 radians and the phase value of the latter distance gate is 0.1 radians, direct subtraction gives negative 5.4 radians, which is less than negative circular constant (about negative 3.14), so it is corrected to negative 5.4 plus twice the circular constant (about 6.28), resulting in about 0.88 radians.

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

[0081] S3: When the phase difference value exceeds the preset linear threshold range, the dielectric constant abrupt interface position is identified, specifically implemented as:

[0082] The phase difference sequence array is stored in a designated area of the memory of the digital signal processor, and the array contains the phase difference values of the adjacent range gates calculated in the previous step, and the length of the array is one less than the total number of range gates. For example, when the total number of range gates is 2048, the array contains 2047 phase difference values. The processing unit starts from the beginning of the array and reads each phase difference value in index order. The reading operation is performed by direct memory access, and each time a 32-bit floating-point number is read until all elements in the array are traversed. During the traversal process, the number of the adjacent range gate corresponding to the currently processed phase difference value is recorded, which is obtained by adding 1 to the array index value to obtain the front range gate number and adding 2 to the array index value to obtain the rear range gate number.

[0083] The preset linear threshold range is defined by two boundary values: the lower limit value is set to -0.3 radian, and the upper limit value is set to +0.3 radian. The threshold range is set according to the atmospheric dielectric constant gradient model, and the phase difference value is usually in this interval when the atmospheric parameters change continuously. The threshold parameter is stored in the non-volatile memory and loaded into the register group when the system starts. The processing unit compares the currently traversed phase difference value with the lower limit value and the upper limit value of the preset linear threshold range: first, it is judged whether the phase difference value is less than -0.3 radian, and second, it is judged whether it is greater than +0.3 radian. The comparison operation is realized by using a floating-point comparator hardware circuit, and a binary state flag is output.

[0084] When the double boundary comparison result shows that the phase difference value is less than -0.3 radian or greater than +0.3 radian, it is determined that there is a dielectric constant mutation at this position. At this time, the position information of the mutation is recorded: the two range gate numbers corresponding to the current phase difference value are obtained, and the range gate with the smaller number is taken as the marker position. According to the range gate number, the pre-stored spatial coordinate mapping table is queried, which is generated when the radar is initialized and stores the three-dimensional spatial coordinates of each range gate. The spatial coordinates are based on the East-North-Sky coordinate system, with the origin at the phase center of the radar antenna. The east coordinate value is calculated by multiplying the radar azimuth angle sine value by the slant range, the north coordinate value is calculated by multiplying the radar azimuth angle cosine value by the slant range, and the sky coordinate value is calculated by multiplying the radar elevation angle sine value by the slant range. For example, the range gate with number 100, when the radar azimuth angle is 30 degrees, the elevation angle is 5 degrees, and the slant range is 3000 meters, its spatial coordinates are east 1500 meters, north 2598 meters, and sky 261 meters.

[0085] The obtained 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 realized by vector dot product calculation: a position vector is formed with the radar position as the origin and the target point spatial coordinates as the terminal point; the dot product of the position vector and the radar beam pointing unit vector is calculated, and the obtained scalar value is the projection slant range of the target point on the detection path. The projection slant range is accurate to 1 meter resolution. Finally, the calculated slant range value is marked as the dielectric constant discontinuity interface position, and stored in the discontinuity position record array, while the corresponding millisecond-level detection timestamp is recorded. For example, the spatial coordinates (1500, 2598, 261) in the beam direction with an azimuth angle of 30 degrees and an elevation angle of 5 degrees have a projection slant range of 3000 meters.

[0086] S4: Recursive graph quantization analysis is performed on the dielectric constant discontinuity interface positions of a continuous preset number of periods, the diagonal line structure features of the recursive graph are extracted as the determinacy coefficient, and the vertical segment distribution is calculated as the laminar flow parameter, which is specifically implemented as:

[0087] The dielectric constant discontinuity interface position records of a continuous preset number of periods are extracted from the storage unit of the radar system. The preset number of periods is determined according to the typical time scale of atmospheric motion, for example, 10 complete radar scanning periods are selected continuously. Each period corresponds to a complete scan of the radar beam along the detection path, and the period length is determined by the radar pulse repetition frequency, for example, when the pulse repetition frequency is 2000 Hz, the period is 0.0005 seconds multiplied by the total number of range gates. All the dielectric constant discontinuity interface positions detected in each period are sorted in chronological order: the positions detected in the first period are placed first, and the positions detected in the last period are placed last, forming an interface position time sequence. Each position point in the sequence contains three-dimensional spatial coordinate information (east coordinate value, north coordinate value, and sky coordinate value) and an accurate time marker to the nearest millisecond. For example, the total length of the sequence may be 150 position points, with the first position point corresponding to the earliest detection time and the 150th position point corresponding to the latest detection time.

[0088] The recursive threshold is used to determine the degree of position similarity, and its setting method is as follows: first, calculate the straight-line distance between all pairs of position points in the interface position time sequence, and find the maximum distance value; then set 5% of the maximum distance value as the recursive threshold. For example, when the maximum distance is 500 meters, the recursive threshold is set to 25 meters. This threshold parameter is stored in a floating-point register. The straight-line distance is calculated using the three-dimensional space distance formula: take the square of the difference between the east coordinates of two position points, add the square of the difference between the north coordinates, and add the square of the difference between the sky coordinates, then take the square root of the sum to get the distance value. The calculation process retains one decimal place of accuracy.

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

[0090] Diagonal line segments formed by consecutive recurrence points in the recurrence plot are identified. The diagonal direction is defined as the direction parallel to the main diagonal of the grid, i.e. the direction where the row number and column number difference is constant. Along each straight line parallel to the main diagonal, when at least two consecutive recurrence points appear, it is recorded as a line segment. The length distribution of all such line segments is counted: set the length statistical interval as short line segments (2-3 consecutive points), medium line segments (4-6 consecutive points), long line segments (7-10 consecutive points), and super-long line segments (more than 11 consecutive points); count the proportion of the number of line segments in each interval to the total number of line segments; form a feature vector containing four proportion values as the determinacy coefficient. For example, the statistical results may be that short line segments account for 15%, medium line segments account for 35%, long line segments account for 45%, and super-long line segments account for 5%.

[0091] Vertical line segments formed by consecutive recurrence points in the recurrence plot are detected. The vertical direction is defined as the direction where the row number remains constant and changes along the column number in the increasing direction. For each fixed row number (e.g. the 80th row), scan along the column number direction: when at least two consecutive recurrence points appear in the row, it is recorded as a vertical line segment. The length distribution of all vertical line segments is counted: the same four length intervals as the diagonal line are used; the proportion of the number of line segments in each interval is calculated; a feature vector containing four proportion values is formed as the laminarity parameter. For example, the results may be that short vertical line segments account for 25%, medium vertical line segments account for 40%, long vertical line segments account for 30%, and super-long vertical line segments account for 5%. The final output determinacy coefficient and laminarity parameter vector are stored in the register group for subsequent stability level calculation.

[0092] S5: Determine the dynamic stability level according to the determinacy coefficient and the laminarity parameter, and execute the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold, which is specifically implemented as:

[0093] The determinacy coefficient and the laminarity parameter are stored in the processing unit registers, both of which are eigenvectors containing four proportion values generated by the recursive graph analysis of the previous step. The two eigenvectors are input into a preset stability mapping model. The model is constructed as follows: first, the four proportion values of the determinacy coefficient eigenvector are weighted and summed according to the weight coefficients 0.4, 0.3, 0.2, and 0.1 to obtain a comprehensive value (for example, input [0.15, 0.35, 0.45, 0.05];

[0094] The calculation is 0.4 x 0.15 + 0.3 x 0.35 + 0.2 x 0.45 + 0.1 x 0.05 = 0.32); and the laminarity parameter comprehensive value is calculated in the same way (for example, [0.25, 0.40, 0.30, 0.05] = 0.35). In a two-dimensional plane with the determinacy coefficient comprehensive value as the horizontal coordinate and the laminarity parameter comprehensive value as the vertical coordinate, five dynamic stability level areas are divided (1st level is the most stable, and 5th level is the least stable), and the area boundaries are determined based on meteorological experimental data. When the input eigenvector is input, the comprehensive value is calculated and the coordinate point is located (for example, 0.32, 0.35), and the level value corresponding to the area where the point is located is searched (for example, if it is located in the 3rd level area, the output value is 3).

[0095] The preset stability threshold is read, which is determined according to the radar operating wavelength (for example, 5 cm wavelength) and the medium uniformity. The medium uniformity is obtained by calculating the ratio of the standard deviation to the mean of the dielectric constant in the last 10 minutes (for example, the ratio is 0.25, which belongs to the 3rd level of uniformity). In the pre-stored data table, the threshold value corresponding to the 5 cm wavelength and the 3rd level of uniformity is set to 2.5. The current threshold value is updated every minute and stored in a special register.

[0096] The dynamic stability level value is compared with the preset stability threshold: a digital comparator is used to determine the size relationship of the values. If the dynamic stability level value is less than the preset stability threshold (for example, 3.0 < 2.5 is not true), a logical false signal is output; otherwise, a logical true signal is output. The comparison result is recorded in the state flag bit.

[0097] When the state flag bit shows a logical false signal (i.e., the dynamic stability level value is less than the threshold value), a 32-bit control instruction is generated. The first 16 bits of the instruction are fixed as the strategy code 64521, and the last 16 bits are the current scan period number (for example, 258). The instruction is transmitted to the compensation unit through the bus, and an interrupt request with a priority of 2 is generated at the same time.

[0098] According to the instruction, the compensation parameter matrix is initialized: the number of matrix rows is equal to the total number of distance gates (for example, 2048 rows), and contains two columns of data. The first column of the distance gate compensation vector is initialized to 2048 1.0 values, and the second column of the time step compensation vector is also initialized to 2048 1.0 values. The following operations are sequentially performed in the memory: first, write the current cycle number; then, write 2048 distance gate compensation initial values in succession; then, write 2048 time step compensation initial values in succession. Finally, a 2048 row × 2 column matrix is formed for subsequent use.

[0099] S6: Based on the path integral compensation strategy, the original echo signal is phase corrected, a corrected signal is generated, and coherent accumulation processing is performed, and a high-precision motion parameter inversion result is output. The specific implementation is as follows:

[0100] The path integral compensation parameter matrix is read from the dual-port random access memory, which contains 2048 rows and 2 columns of data. The first column stores the distance gate compensation amount vector, which contains 2048 floating point values; the second column stores the time step compensation amount vector, which also contains 2048 floating point values. The extraction operation is performed by direct memory access: first, locate the matrix base address, and start reading 2048 distance gate compensation amounts from the base address offset 0 bytes; then, start reading 2048 time step compensation amounts from the base address offset 16384 bytes. For example, the 100th element of the distance gate compensation amount vector corresponds to the 100th distance gate, and the initial value is 1.0.

[0101] The distance gate compensation amount vector is multiplied element by element with the phase component of the original echo signal. The original echo signal is stored in the cache area of the digital signal processor and contains 2048 complex number data, and the imaginary part of each complex number is the phase component. The operation process is as follows: for each distance gate index (from 1 to 2048), multiply the value at the corresponding position in the distance gate compensation amount vector by the phase component value of the same distance gate in the original echo signal. For example, for the 50th distance gate, the compensation amount 0.98 is multiplied by the phase component 0.75 radians to obtain 0.735 radians. The operation result forms a primary compensation signal, which is stored in a newly opened memory area and has the same data structure as the original echo signal.

[0102] The time step compensation amount vector is multiplied element by element with the phase component of the primary compensation signal. The values in the time step compensation amount vector are arranged in the order of detection time. The operation process is as follows: for each distance gate index (from 1 to 2048), multiply the value at the corresponding position in the time step compensation amount vector by the phase component value of the same distance gate in the primary compensation signal. For example, for the 150th distance gate, the time step compensation amount 1.02 is multiplied by the primary phase component 0.82 radians to obtain 0.8364 radians. The operation result generates a corrected signal, which retains the real part of the original echo signal unchanged and only the phase component is corrected, and is stored in the output buffer area of the digital signal processor.

[0103] The coherent accumulation processing of the discrete Fourier transform is performed on the corrected signal. The processing process includes three steps: first, the corrected signals of 32 consecutive pulse periods are arranged in time sequence to form a 32*2048 complex matrix; then, for each distance gate index (from 1 to 2048), the discrete Fourier transform is performed on 32 complex numbers along the time dimension; finally, the modulus square value of the transform result is calculated to obtain the Doppler spectrum energy distribution of each distance gate. For example, the spectrum output by the 1000th distance gate contains 32 energy values, and the frequency resolution is determined by the pulse repetition frequency (2000 Hz repetition frequency corresponds to 62.5 Hz resolution).

[0104] The radial velocity value is calculated according to the peak position of the Doppler spectrum energy distribution. First, the position index of the maximum energy value in each distance gate spectrum is located, for example, the 800th distance gate appears a peak at index 15. The velocity calculation formula is: the radial velocity is equal to the peak position index multiplied by the wavelength and multiplied by the pulse repetition frequency divided by the number of discrete Fourier transform points. For example, when the wavelength is 0.1 meters, the pulse repetition frequency is 2000 Hz, and the number of transform points is 32, the velocity corresponding to index 15 is 15*0.1*2000 / 32=93.75 m / s. The calculation result is stored in the result register group as the high-precision motion parameter inversion result, which is one of the key inputs for realizing the synchronous inversion of the vertical profiles of atmospheric temperature, humidity, wind direction and rain intensity by a single device, and the accuracy directly affects the quality of the final profile product.

[0105] The high-precision motion parameter inversion result is transmitted to the atmospheric parameter solving unit for joint inversion of the vertical distribution of atmospheric temperature, humidity, wind direction and rain intensity, specifically: the vertical distribution of atmospheric temperature profile, relative humidity profile, liquid water profile, wind speed profile, wind direction profile, raindrop spectrum and rainfall rate, etc. The solving process is based on the meteorological fluid mechanics equation set, taking the radial velocity field as the key input parameter, and optimizing the atmospheric state estimation by the variational assimilation algorithm. For example, in the boundary layer wind field inversion, the radial velocity data is used to constrain the horizontal momentum term of the three-dimensional wind field reconstruction equation.

[0106] Embodiment 2: Figure 2 The structure diagram of the all-weather atmospheric temperature, humidity, wind and rain profile joint detection system is given, and the all-weather atmospheric temperature, humidity, wind and rain profile joint detection system comprises the following modules:

[0107] An original echo acquisition module is configured to acquire original echo signals received by the radar on a detection path;

[0108] A phase difference calculation module is configured to extract a phase change sequence from the original echo signals and calculate phase difference values of adjacent distance gates;

[0109] An interface position recognition module is configured to recognize the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range.

[0110] A recurrence plot analysis module is configured to perform recurrence plot quantitative analysis on the dielectric constant mutation interface positions in a continuous preset period of time, extract the diagonal line structure features of the recurrence plot as the determinacy coefficient, and calculate the vertical segment distribution as the laminar flow parameter.

[0111] A compensation strategy decision module is configured to determine the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and execute the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold.

[0112] A signal correction processing module is configured to perform phase correction on the original echo signal based on the path integral compensation strategy, generate a corrected signal, and perform coherent accumulation processing to output a high-precision motion parameter inversion result.

[0113] The calculations involved in the embodiments are all de-dimensioned numerical calculations, and the preset parameters and threshold values in the calculations are set by a person skilled in the art according to actual conditions.

[0114] The above embodiments can be realized wholly or partially by software, hardware, firmware or any other combination. When realized by software, the above embodiments can be realized wholly or partially in the form of a computer program product.

[0115] Those skilled in the art can realize that the modules and algorithm steps of the examples described in combination with the embodiments disclosed herein can be realized by electronic hardware or a combination of computer software and electronic hardware. Whether the functions are realized in hardware or software depends on the specific application and the constraints of the technical solution. A person skilled in the art can use different methods to realize the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0116] In addition, the functional modules in each of the embodiments of the present application can be integrated in one processing module, or each module can exist physically independently, or two or more modules can be integrated in one module.

[0117] In several embodiments provided in the present application, it should be understood that the disclosed system, device and method can be implemented in other manners. For example, the division of the above-described device embodiment is only a logical function division, and there can be another division manner for actual implementation, for example, multiple devices or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the displayed or discussed mutual couplings or direct couplings or communication connections between different parts can be indirect couplings or communication connections through some interfaces, devices or modules, and can be in electrical, mechanical or other forms.

[0118] The above describes only the specific embodiments of the present application, but the protection scope of the present application is not limited thereto, and any modification or replacement within the technical range disclosed by the present application can be easily thought by any person skilled in the art, and should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

[0119] Finally: the above described is only the preferred embodiment of the present application, and is not used to limit the present application, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the protection scope of the present application.

Claims

1. A method for all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding, characterized in that, The method comprises the following steps: S1: obtaining the original echo signal received by the radar on the detection path; S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates; S3: identifying the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range; S4: performing recursive graph quantization analysis on the dielectric constant mutation interface positions of a continuous preset number of periods, extracting the diagonal line structure features of the recursive graph as the determinacy coefficient, and calculating the vertical segment distribution as the laminar flow parameter; S5: determining the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and executing the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold; S6: performing phase correction on the original echo signal based on the path integral compensation strategy, generating the corrected signal and performing coherent accumulation processing, and outputting the high-precision motion parameter inversion result.

2. The all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding method according to claim 1, characterized in that, The method comprises the following steps: S1: obtaining the original echo signal received by the radar on the detection path; S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates; S3: identifying the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range; S4: performing recursive graph quantization analysis on the dielectric constant mutation interface positions of a continuous preset number of periods, extracting the diagonal line structure features of the recursive graph as the determinacy coefficient, and calculating the vertical segment distribution as the laminar flow parameter; 3. The all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding method according to claim 2, characterized in that, S5: determining the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and executing the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold; S6: performing phase correction on the original echo signal based on the path integral compensation strategy, generating the corrected signal and performing coherent accumulation processing, and outputting the high-precision motion parameter inversion result. The method comprises the following steps: S1: obtaining the original echo signal received by the radar on the detection path; S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates; 4. The all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding method according to claim 3, characterized in that, S3: identifying the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range; S4: performing recursive graph quantization analysis on the dielectric constant mutation interface positions of a continuous preset number of periods, extracting the diagonal line structure features of the recursive graph as the determinacy coefficient, and calculating the vertical segment distribution as the laminar flow parameter; S5: determining the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and executing the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold; S6: performing phase correction on the original echo signal based on the path integral compensation strategy, generating the corrected signal and performing coherent accumulation processing, and outputting the high-precision motion parameter inversion result. The method comprises the following steps:

5. The all-weather atmospheric temperature, humidity, wind, and precipitation profile combined sounding method according to claim 4, characterized in that, S1: obtaining the original echo signal received by the radar on the detection path; S2: extracting the phase change sequence from the original echo signal and calculating the phase difference value of adjacent distance gates; S3: identifying the dielectric constant mutation interface position when the phase difference value exceeds the preset linear threshold range; S4: performing recursive graph quantization analysis on the dielectric constant mutation interface positions of a continuous preset number of periods, extracting the diagonal line structure features of the recursive graph as the determinacy coefficient, and calculating the vertical segment distribution as the laminar flow parameter; S5: determining the dynamic stability level according to the determinacy coefficient and the laminar flow parameter, and executing the path integral compensation strategy when the dynamic stability level is lower than the preset stability threshold; S6: performing phase correction on the original echo signal based on the path integral compensation strategy, generating the corrected signal and performing coherent accumulation processing, and outputting the high-precision motion parameter inversion result.

6. The all-weather atmospheric temperature, humidity, wind, and precipitation profile combined sounding method according to claim 5, characterized in that, The dynamic stability level is determined according to the determinacy coefficient and the laminarity parameter, and the path integral compensation strategy is executed when the dynamic stability level is lower than a preset stability threshold, including: The determinacy coefficient and the laminarity parameter are input into a preset stability mapping model to output a dynamic stability level value; A preset stability threshold is read based on radar wavelength and medium uniformity calibration; The size relationship between the dynamic stability level value and the preset stability threshold is compared; When the dynamic stability level value is less than the preset stability threshold, a path integral compensation strategy generation instruction is activated; The path integral compensation parameter matrix is initialized according to the path integral compensation strategy generation instruction.

7. The all-weather atmospheric temperature, humidity, wind, and precipitation profile combined sounding method according to claim 6, characterized in that, The preset stability mapping model is implemented by the following methods: A two-dimensional feature space of the determinacy coefficient and the laminarity parameter is established; The contour area of the dynamic stability level value is divided in the two-dimensional feature space; The dynamic stability level value corresponding to the input determinacy coefficient and the laminarity parameter is calculated by linear interpolation.

8. The all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding method according to claim 6, characterized in that, The original echo signal is phase-corrected based on the path integral compensation strategy, a corrected signal is generated and coherent accumulation processing is performed, and a high-precision motion parameter inversion result is output, including: The distance gate compensation amount vector and the time step compensation amount vector are extracted from the path integral compensation parameter matrix; The distance gate compensation amount vector and the phase component of the original echo signal are subjected to Hadamard product operation to generate a primary compensation signal; The time step compensation amount vector and the phase component of the primary compensation signal are subjected to Hadamard product operation to generate a corrected signal; The corrected signal is subjected to discrete Fourier transform coherent accumulation processing to obtain a 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 precipitation profile combined sounding method according to claim 8, characterized in that, Wherein, The high-precision motion parameter inversion result is used for atmospheric temperature, humidity, wind and rain profile joint inversion.

10. An all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding system for implementing the all-weather atmospheric temperature, humidity, wind and precipitation profile combined sounding method according to any one of claims 1-9, characterized in that, It includes the following modules: An original echo acquisition module for acquiring an original echo signal received by a radar on a detection path; A phase difference calculation module for extracting a phase change sequence from the original echo signal and calculating a phase difference value of adjacent distance gates; An interface position identification module for identifying a dielectric constant mutation interface position when the phase difference value exceeds a preset linear threshold range; A recurrence plot analysis module for performing recurrence plot quantization analysis on dielectric constant mutation interface positions for a continuous preset number of periods, extracting recurrence plot diagonal line structure features as a determinacy coefficient, and calculating a vertical segment distribution as a laminarity parameter; A compensation strategy decision module for determining a dynamic stability level according to the determinacy coefficient and the laminarity parameter, and executing a path integral compensation strategy when the dynamic stability level is lower than a preset stability threshold; A signal correction processing module for phase-correcting the original echo signal based on the path integral compensation strategy, generating a corrected signal and performing coherent accumulation processing, and outputting a high-precision motion parameter inversion result.

Images