Multi-mode multi-frequency GNSS-IR soil humidity inversion method based on PID search and dynamic time warping

Abstract

The invention discloses a multimode multi-frequency GNSS-IR soil humidity inversion method based on PID search and dynamic time warping, and belongs to the technical field of GNSS interference reflection remote sensing and soil humidity monitoring. The method is based on GNSS original observation data of a geodetic measurement type receiver, and comprises the following steps: firstly, extracting signal-to-noise ratios, elevation angles and azimuth angles of satellites, and separating reflected signals through polynomial fitting; optimizing a nonlinear least square algorithm by adopting a PID (Proportion Integration Differentiation) search algorithm to obtain accurate amplitude and phase; and finally, selecting available satellites by adopting a dynamic time warping satellite selection method and an amplitude weighted fusion strategy, and carrying out effective fusion in combination with the amplitude. According to the method, phase jump can be effectively inhibited, low-quality satellites are eliminated, accurate inversion of the soil humidity can be realized without depending on massive prior information, and the method has higher applicability.

Description

Technical Field

[0001] The present invention belongs to the technical field of GNSS interferometric reflection remote sensing and soil moisture monitoring, and more specifically relates to a multi-mode and multi-frequency GNSS-IR soil moisture inversion method based on PID search and dynamic time warping. Background Art

[0002] Soil moisture is an important indicator of the water content in the soil. It is of great significance for agricultural production, water resources management, ecological and environmental protection, etc. Soil moisture is an important basis for preventing disasters such as floods, droughts and climate change. In the past, traditional soil moisture measurement methods had disadvantages such as locality, discontinuity, high equipment cost and high requirements for professionals. In recent years, the inversion of surface environmental parameters using the characteristic variables (amplitude and phase) of GNSS reflection signals has become a new type of microwave remote sensing method. This technology is called global navigation satellite system interferometric reflectometry (GNSS-IR). It has the advantages of non-destructive, real-time continuous inversion, high temporal and spatial resolution, automated inversion, strong anti-interference ability, high data accuracy and versatility. It has achieved effective inversion of parameters such as vegetation cover, soil moisture, snow thickness and sea level.

[0003] A large number of studies have shown that GNSS-IR technology has shown significant effectiveness in the field of soil moisture monitoring using satellite data. This technology conventionally relies on the standard cosine model, extracts the frequency information of the multipath effect from the signal-to-noise ratio reflection component through the Lomb-Scargle spectrum analysis method, and further uses the nonlinear least squares (NLS) method to analyze the characteristics of SNRm, namely phase and amplitude, and then uses the phase as a parameter to track the change of soil moisture. However, in the process of solving the characteristic parameters, negative amplitude and abnormal phase jump problems often occur, which becomes a major challenge for high-precision inversion. The current response strategy is mostly to remove abnormal phase observations or abandon the affected satellite trajectories, which not only greatly reduces the amount of available satellite data, limits the ability to capture the key change trajectory of soil moisture, but also breaks the natural continuity of the phase parameter time series, affecting the consistency and accuracy of data analysis. In addition, in the process of multi-satellite fusion, the problem of satellite selection and the missing of some phase data are also key factors affecting the inversion accuracy. The reception of GNSS signals is affected by many factors such as the environment and weather. Satellite signals may be lost or incomplete due to obstruction, noise or multipath effects. How to effectively select suitable satellites has become an important issue. Summary of the invention

[0004] To solve the above problems, the present invention proposes a method for providing NLS initial values ​​through a PID search algorithm (PID-based search algorithm, PSA), automatically solving the amplitude and phase of the reflected signal component in the GNSS-IR inversion of soil moisture, and then inverting through the dynamic time warping (Dynamic Time Warping, DTW) satellite selection and amplitude weighted fusion strategy. The method is based on the frequency band data of L1, L2 of GPS and G1, G2 of GLONASS. First, for a certain frequency band of a certain system, after obtaining the satellite signal-to-noise ratio data, the reflected signal component is obtained by low-order polynomial fitting and removing the direct signal, and the window smoothing is performed to obtain the smoothed reflected component. The upper and lower limits of the amplitude are determined according to the peak point of the smoothed reflected component, and it is applied to the PSA for optimization to obtain the amplitude and phase. Afterwards, in order to improve the stability and accuracy of soil moisture inversion, the DTW distance between satellite phases is measured. Since DTW can process data of unequal lengths, similarity measurement can be performed when there is a partial loss in the satellite phase sequence. Then, satellites with poor similarity or missing data are eliminated to ensure that high-quality satellite data is retained for fusion. Amplitude normalization weighting is adopted to fuse multiple satellites in different frequency bands, and then they are converted through inversion formulas to obtain stable inversion results.

[0005] The present invention achieves the above-mentioned purpose through the following technical solutions. A multi-mode and multi-frequency GNSS-IR soil moisture inversion method based on PID search (PSA) and dynamic time warping (DTW) comprises the following steps:

[0006] Step 1: Obtain GNSS data and soil moisture. This experiment selects the observation data of the PBO network P041 station in Colorado, USA as the experimental data. The soil moisture of the station is provided by the Marshall Colorado 006 station of COSMOS (http: / / cosmos.hwr.arizona.edu / ).

[0007] Step 2: GNSS data preprocessing. Calculate GNSS satellite data through Rtklib to obtain the signal-to-noise ratio (SNR) data, satellite elevation angle, azimuth angle and other information of each satellite, and then select all available satellites of a single system and a single frequency band to participate in the subsequent calculations.

[0008] Step 3: Separation of GNSS satellite reflected signals. Use a third-order polynomial to fit the SNR data, remove the trend term, and retain the reflected component.

[0009] Step 4: Use PSA and find the best parameters. The algorithm is repeatedly driven through a preset number of cycles to explore the parameter space. After each run, the optimal objective function value (i.e., the lowest fitness value) and corresponding parameters of the current iteration are obtained and recorded. Through multiple iterations of the PSA algorithm, the appropriate initial values ​​of amplitude and phase are obtained and passed as input to NLS for precise fitting, and finally the phase and amplitude are solved.

[0010] Step 5: The time of each satellite phase is unified. The initial phase of each satellite has a corresponding SNR data sequence. The middle moment of the sequence period is used as the satellite phase acquisition time and recorded. Finally, all amplitude and phase values ​​of a single satellite within the specified time range are obtained. According to the acquisition time of a single satellite phase, it is divided into corresponding annual days, and the phase value within a single day is processed by the equal-weighted average method as the phase value corresponding to the annual day. After that, the time, phase and other related information are saved to obtain the relevant data of a single system single frequency band multi-satellite.

[0011] Step 6: Select appropriate satellites. By calculating the DTW distance between the satellite phase sequences in each frequency band and normalizing them, appropriate satellites can be selected. This method is suitable for processing missing data.

[0012] Step 7: Single-band normalized amplitude weighted fusion: For the selected satellites in each frequency band, normalized amplitude weighted fusion is performed to improve the accuracy of fusion.

[0013] Step 8: Multi-frequency fusion. The moisture estimation value sequence of each frequency band obtained in step 7 is equally weighted to reduce possible systematic deviations, and the multi-frequency band soil moisture estimation value is obtained, which is compared with the measured soil moisture accuracy evaluation.

[0014] This invention proposes a comprehensive solution to the limitations of global navigation satellite system interferometry (GNSS-IR) technology in soil moisture monitoring. First, PSA optimizes the phase solution process and suppresses phase jumps by accurately optimizing the initial value. After that, the normalized dynamic time warping (DTW) distance selection strategy is used to measure the similarity between satellite phase sequences and eliminate satellites with poor quality to improve the reliability of fusion in the case of missing data or large signal interference. Finally, combined with the amplitude normalization weighted method, the selected multiple satellites are weightedly fused, and the stability and accuracy of the inversion results are further improved through the reasonable integration of multi-satellite and multi-frequency data. Therefore, this research scheme has a clear theoretical basis and scientific research ideas, which can effectively overcome the phase jump problem of the current GNSS-IR technology in soil moisture inversion, and can make full use of existing data to improve the inversion accuracy and stability without almost any prior conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a schematic diagram of the process of the present invention.

[0016] Figure 2 Schematic diagram of the multipath effect geometric model.

[0017] Figure 3 Schematic diagram of the relationship between SNR (dB-Hz) and satellite elevation angle (0°-90°).

[0018] Figure 4 This is a schematic diagram of the relationship between the signal-to-noise ratio (Volts) and the sine value of the altitude angle (sin0° - sin25°).

[0019] Figure 5 Schematic diagram for optimization simulation comparison.

[0020] Figure 6 Schematic diagram of the comparison of the absolute values ​​of the correlation between the phase sequence and humidity solved by the two methods.

[0021] Figure 7 This is a schematic diagram of the four-band fusion effect in two time periods. DETAILED DESCRIPTION

[0022] In order to make the purpose, technical scheme and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0023] Figure 1 As shown, it is the implementation process of the present invention, and the specific steps are as follows:

[0024] Step 1: Obtain GNSS observation data and soil moisture. This experiment selected the GPS L1, L2 and GLONASS G1, G2 frequency band data of the PBO network P041 station in Colorado, USA in 2016 DOY2-321 (time period 1) and 2017 DOY82-365 (time period 2) as experimental data. The reference soil moisture is provided by the Marshall Colorado 006 weather station of COSMOS (http: / / cosmos.hwr.arizona.edu / ), which provides the measured rainfall data and soil moisture observation values ​​of the station. The weather station is located about 105 meters west of the P041 station.

[0025] Step 2: GNSS observation data preprocessing: The original observation data obtained by the GNSS receiver is used to extract the SNR, altitude angle, azimuth and other information of all satellites using Rtklib software, and the generated file is saved and then read using Matlab.

[0026] Step 3: Separation of GNSS satellite reflected signals. The composite signal recorded in the GNSS receiver mainly appears in the form of SNR. SNR is easily affected by the external environment, including the gain characteristics of the receiver antenna, the multipath effect of the signal, and background noise, such as Figure 2 shown. Figure 3 The relationship between SNR and satellite elevation angle is shown. When the satellite elevation angle is between 0° and 20°, the multipath effect has a significant interference on the SNR. When the satellite elevation angle continues to increase, the data fluctuation caused by the multipath effect gradually decreases, and finally is eliminated under the antenna gain effect. This phenomenon reveals the significant relationship between SNR and multipath effect. The direct signal components of GPS and GLONASS satellites are fitted by third-order polynomials, and the reflected signals of each satellite are obtained by subtracting the direct signal from the SNR, and then PSA is used to automatically solve the amplitude and phase of the satellite. Figure 4 This is a schematic diagram of the relationship between the signal-to-noise ratio (Volts) and the sine value of the altitude angle (sin0° - sin25°). In the figure, we can see the reflected signal after the trend term is removed.

[0027] Step 4: Use PSA and find the best parameters. The core of PSA is to control the update direction and amplitude of the solution through proportional, integral and differential adjustment. The update of each solution is given by the following formula:

[0028] (1)

[0029] In the formula, Indicates that the solution is in the iteration The deviation is defined as the difference between the current solution and the optimal solution. , , are the proportional, integral and differential coefficients respectively, , , is a random number vector with n rows and 1 column, 0 to 1, used to balance exploration and development. After updating the solution through the above formula, PSA introduces a conditional factor , to prevent the solution from converging to a local optimum too early:

[0030] (2)

[0031] in is the maximum number of iterations, is the attenuation factor, which decreases with the number of iterations, gradually shifting from global search to local convergence. is a random number vector with n rows and 1 column from 0 to 1, is the Levy flight function, which is used to enhance the global search capability. The final update formula of the solution is:

[0032] (3)

[0033] in is a weight factor used to balance the contribution of different regulation strategies.

[0034] The PSA algorithm optimizes the amplitude and phase parameters to provide the best initial values ​​for NLS. This process consists of three main stages: initialization, objective function definition, and iterative optimization. Initialization involves setting the population size (50), the number of iterations (100), and the amplitude bounds (90-110% of the signal peak). The objective function minimizes the L2 norm of the error between the modeled signal and the observed signal. Finally, the optimal amplitude and phase are determined through iterative optimization and used as the input to the NLS for accurate fitting.

[0035] In order to observe the characteristics of the NLS solution and the feasibility of PSA optimization, we designed a set of simulation experiments to evaluate the impact of different initial values ​​on the stability and accuracy of the NLS solution results:

[0036] First, the phase range of the signal is set to a complete cycle, increasing in steps of 0.05, to generate a set of phase values ​​to form a phase set, and each phase value corresponds to the initial phase of a simulated signal. Next, 100 sine values ​​are uniformly generated within the elevation angle range of 5° to 20° as discrete sampling points of the signal within this range. For each phase value, combined with the true amplitude 15, frequency 17 and damping factor 40, a noise-free gradually attenuated signal is generated point by point, and random noise with a standard deviation of 2 is superimposed to obtain a simulated reflection signal. Each simulated signal is composed of the above 100 sampling points, and its morphology is affected by the initial phase and random noise. The signal finally generated has the characteristic that the amplitude gradually decays with increasing satellite elevation angle, which is close to the reflection signal in the real environment.

[0037] Figure 5This is a comparison chart of optimization simulation. In the figure, we can intuitively find the following fact: when using NLS solution, its accuracy is highly dependent on the closeness of the initial value to the true value. When the given initial value is close to the true value, NLS can accurately fit the phase and amplitude within the "correct range". However, when the initial value is far from the true value, the NLS solution value will produce a significant jump, accompanied by a drastic change in the fitted amplitude, including the occurrence of negative amplitude, and this jump is irregular. When the set initial value of the amplitude is close to the true amplitude, the correct range of the NLS phase is larger. PSA optimization provides NLS with a more reliable initial value, thereby expanding the correct range of the phase and restoring the amplitude to normal. In actual applications, PSA optimization also suppresses phase jumps, Figure 6 The comparison of the absolute values ​​of the correlation between the phase sequences solved by the two methods and soil moisture in two time periods is shown, where the red numbers represent the correlation calculation after the missing phases were assigned to 0.

[0038] Step 5: The time of each satellite phase is unified. The initial phase of each satellite has a corresponding SNR data sequence. The middle moment of the sequence period is used as the satellite phase acquisition time and recorded. Finally, together with all the amplitude and phase values ​​of a single satellite within the specified time range obtained by PSA, they are divided into corresponding annual days according to the acquisition time of the single satellite phase, and the phase value within a single day is processed by the equal-weighted average method as the phase value corresponding to the annual day. Finally, we can obtain the phase data of all single satellites of GPS L1, L2 and GLONASS G1, G2 on each day.

[0039] Step 6: Select appropriate satellites. First, preliminarily sort out the satellite phase data and remove satellites with serious phase loss, because they are difficult to provide valuable information in subsequent fusion and may introduce large uncertainties. Then, for the remaining satellites, deeply explore their phase sequence characteristics and amplitude information for processing.

[0040] For each remaining satellite, the DTW distance is used to measure the similarity of its phase sequence with other satellites. Assume that the phase sequences of the two satellites are and , its DTW distance is:

[0041] (4)

[0042] in, is the sum of the distances of all points on the path, In order to eliminate the scale difference of the original DTW distance between different satellites, the DTW distance is normalized:

[0043] (5)

[0044] Among them, is the normalized distance value, and are the minimum and maximum values ​​in the distance matrix respectively. Next, set the threshold and scale parameters For each satellite, the normalized DTW distance between it and the other satellites is calculated, and it is determined whether the distance is less than the threshold. In this way, the "similar satellite ratio" of each satellite is determined. :

[0045] (6)

[0046] in, is the total number of satellites, Defined as:

[0047] (7)

[0048] when When , it indicates that the normalized DTW distance between the satellite and most satellites in the frequency band is small, and it has a high similarity with other satellites, and can participate in the subsequent fusion. In this way, the satellite phase sequence of each frequency band is screened to eliminate the poor satellites. th and P should be adjusted according to the characteristics of the frequency band phase sequence to ensure that a sufficient number of high-quality satellites can be retained.

[0049] Step 7: Single-band normalized amplitude weighted fusion. After completing the satellite selection, enter the weighted fusion stage. Assume that the total number of satellites involved in the fusion is n. For the i-th satellite (i=1,2,...,n), its amplitude is recorded as , the weight of the satellite for:

[0050] (8)

[0051] Here j is the sum index, which is used to traverse all satellites involved in the fusion. Based on this weight, the humidity reflected on that day in this frequency band is for:

[0052] (9)

[0053] in represents the i-th satellite. The humidity value is calculated according to the algorithm proposed by Chew et al. (2014):

[0054] (10)

[0055] S is the slope between phase and SM. To avoid the influence of abnormal minimum value in phase sequence on estimation, the mean of the minimum 15% in phase sequence is taken as . In the actual data collection process, if the phase data of a certain satellite is missing on a certain day, its amplitude is 0. According to the weight formula, the weight of the satellite is 0, which does not affect the result during the summation calculation, thus avoiding errors caused by missing data or abnormal signals.

[0056] Step 8: Multi-frequency fusion. The moisture estimation value sequence of each frequency band obtained in step 7 is equally weighted to reduce possible systematic deviations, and the multi-frequency band soil moisture estimation value is obtained, which is compared with the measured soil moisture accuracy evaluation. Figure 7 This is the effect diagram of four-band fusion in two time periods.

Claims

1. A multi-mode and multi-frequency GNSS-IR soil moisture inversion method based on PID search and dynamic time warping, characterized in that: The following steps are involved: S1: GNSS data and soil moisture acquisition: randomly select soil moisture monitoring points, receive GNSS raw observation data through geodetic GNSS receivers, and deploy soil moisture collection probes around the measuring stations to monitor soil moisture at a depth of 0-5 cm every day in real time; S2: GNSS data preprocessing: extract information such as signal to noise ratio (SNR), altitude angle and azimuth angle of all satellites based on the original observation data obtained by the GNSS receiver; at the same time, determine the effective monitoring area of ​​soil moisture based on the Huygens-Fresnel principle; S3: GNSS satellite reflection signal separation: Use a third-order polynomial to fit the SNR data, remove the trend term, and retain the reflection signal component; S4: Use the PID search algorithm (PSA) and find the best parameters; repeatedly drive the algorithm through a preset number of cycles to explore the parameter space, and obtain and record the optimal objective function value (i.e., the minimum fitness value) and corresponding parameters of the current iteration after each run; obtain the appropriate initial amplitude and phase values ​​through multiple iterative optimizations of the PSA algorithm, and pass them as input to the nonlinear least squares (NLS) method for precise fitting, and finally solve the phase and amplitude; the fitting formula is as follows: ,in is the signal-to-noise ratio of the reflection component, is the amplitude of the reflected component, is the carrier wavelength, is the phase, is the vertical distance between the antenna phase center and the ground reflection surface; S5: Dynamic Time Warping (DTW) satellite selection and amplitude weighted fusion; use normalized DTW distance to measure the similarity of phase sequences between satellites to eliminate satellites with poor quality and make full use of satellites with partially missing phase data; then use normalized amplitude weighted fusion of single-band satellites; S6: Multi-frequency satellite fusion: The soil moisture sequences obtained from multiple frequency bands are equally weighted to reduce possible systematic deviations and further improve the stability and accuracy of the inversion results.

2. According to claim 1, a multi-mode multi-frequency GNSS-IR soil moisture inversion method based on PID search and dynamic time warping is characterized in that In S4, for the phase in the fitting formula, a more accurate NLS initial value is found using PSA without initial value and passed to NLS for more precise solution to obtain accurate amplitude and phase without jumps.

3. The multi-mode multi-frequency GNSS-IR soil moisture inversion method based on PID search and dynamic time warping according to claim 1 is characterized in that In S5, the normalized DTW distance between satellite phases is used to judge the similarity between satellites to screen high-quality satellites, and the satellite phase sequence does not need to be complete. After obtaining high-quality satellites, normalized amplitude weighted fusion is used to obtain a single-band soil moisture inversion sequence.

Images