CN114839354B - Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy - Google Patents
Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy Download PDFInfo
- Publication number
- CN114839354B CN114839354B CN202210770775.1A CN202210770775A CN114839354B CN 114839354 B CN114839354 B CN 114839354B CN 202210770775 A CN202210770775 A CN 202210770775A CN 114839354 B CN114839354 B CN 114839354B
- Authority
- CN
- China
- Prior art keywords
- signal
- noise ratio
- soil humidity
- noise
- beidou
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000002689 soil Substances 0.000 title claims abstract description 115
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 52
- 238000000691 measurement method Methods 0.000 title claims abstract description 15
- 238000005259 measurement Methods 0.000 claims abstract description 77
- 238000000034 method Methods 0.000 claims abstract description 55
- 238000012935 Averaging Methods 0.000 claims abstract description 5
- 230000009977 dual effect Effects 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims description 33
- 238000001914 filtration Methods 0.000 claims description 22
- 238000012545 processing Methods 0.000 claims description 19
- 230000008569 process Effects 0.000 claims description 18
- 238000005070 sampling Methods 0.000 claims description 12
- 230000003044 adaptive effect Effects 0.000 claims description 10
- 230000009467 reduction Effects 0.000 claims description 9
- 239000013598 vector Substances 0.000 claims description 8
- 238000010183 spectrum analysis Methods 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 6
- 238000012886 linear function Methods 0.000 claims description 5
- 238000001228 spectrum Methods 0.000 claims description 4
- 108010076504 Protein Sorting Signals Proteins 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000013527 convolutional neural network Methods 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 230000000737 periodic effect Effects 0.000 claims description 3
- 238000010187 selection method Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 description 16
- 238000009795 derivation Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000007796 conventional method Methods 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000000523 sample Substances 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012272 crop production Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- DMBHHRLKUKUOEG-UHFFFAOYSA-N diphenylamine Chemical compound C=1C=CC=CC=1NC1=CC=CC=C1 DMBHHRLKUKUOEG-UHFFFAOYSA-N 0.000 description 1
- 238000001035 drying Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000003673 groundwater Substances 0.000 description 1
- 230000002262 irrigation Effects 0.000 description 1
- 238000003973 irrigation Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
- G01N33/24—Earth materials
- G01N33/246—Earth materials for water content
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/33—Multimode operation in different systems which transmit time stamped messages, e.g. GPS/GLONASS
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
- Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Health & Medical Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Mathematical Analysis (AREA)
- Medicinal Chemistry (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Algebra (AREA)
- Computational Mathematics (AREA)
- Geology (AREA)
- Food Science & Technology (AREA)
- Mathematical Optimization (AREA)
- Environmental & Geological Engineering (AREA)
- Pure & Applied Mathematics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Geophysics And Detection Of Objects (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention discloses a Beidou and GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy, which comprises the following steps: s1, establishing a reflecting signal-to-noise ratio phase reference value database; s2, calculating the soil humidity measured by the Beidou system by using a soil humidity measurement model and combining a signal-to-noise ratio phase reference value in a signal-to-noise ratio phase reference value database of the Beidou; s3, calculating the soil humidity measured by the GPS system by using the soil humidity measurement model and combining the signal-to-noise ratio phase reference value in the signal-to-noise ratio phase reference value database of the GPS; and S4, combining the soil humidity measured by the Beidou and GPS dual system through multi-frequency measurement, and averaging the results. The method considers the influence of the signal-to-noise ratio selection length, the effective reflection high fluctuation of the satellite and the frequency characteristic difference of the satellite on the soil humidity measurement precision, effectively solves the problems of low single-system measurement precision, unstable performance and the like, can meet the requirements of soil humidity measurement time and spatial resolution, and can meet the requirements of soil humidity measurement on high precision, high real-time performance and high stability.
Description
Technical Field
The invention relates to the technical field of navigation satellite remote sensing inversion, in particular to a Beidou and GPS soil humidity measuring method based on a sliding algorithm and a weighting strategy.
Background
The soil humidity, also called soil water content, can be used for the research of global atmospheric water circulation, is closely related to the material transfer and energy exchange among the global water circulation, the atmospheric circulation and the biological ecological circulation, and is an important component of the global ecological system. Meanwhile, the soil humidity not only influences regional rainfall, groundwater supply and weather change, but also can be used as an important index for agricultural drought monitoring, and information is provided for irrigation of crops. Therefore, the method and the device for real-time and effective soil humidity measurement are provided and established, which can provide effective information for water vapor ecological research, crop production, management decision and the like, and are one of the problems to be solved in the current soil humidity measurement.
Currently, there are three main types of soil moisture measurement methods. The first method is a traditional measurement method, mainly including a drying method and a probe measurement. Although the traditional method can obtain higher measurement precision, the method needs field operation and processing, has lower time and space convenience rate, cannot meet the requirement of large-range real-time measurement, and has higher labor cost. The second type is a soil humidity measuring method of a professional remote sensing satellite, the method utilizes an environment remote sensing satellite to measure, the measuring range of a large area can be measured, the time and the space resolution are high, but an optical lens and a sensor which are depended by the remote sensing satellite are easily influenced by the environment, and the poor weather such as cloud and haze can cause the measuring failure or the accuracy reduction based on the remote sensing satellite method. The third type is soil humidity measurement based on a GNSS (Global Navigation Satellite System) Navigation Satellite, and the method utilizes phase fluctuation of signal-to-noise ratio of a reflection signal of the GNSS Satellite to invert the soil humidity, so that a large-range soil humidity measurement result is obtained in real time.
However, the existing method does not consider the influence of the signal-to-noise ratio selection length on the effective reflection of the computed satellite and the subsequent phase detection quantity accuracy based on the signal-to-noise ratio, but adopts a fixed signal-to-noise ratio selection length of 5-30 degrees, so that the final soil humidity measurement accuracy is reduced. Meanwhile, the influence of the satellite effective reflection height fluctuating along with time is not considered in the conventional method, so that the accuracy of the signal-to-noise ratio phase of subsequent calculation is reduced. In addition, the existing soil humidity inversion usually adopts single frequency for inversion, and the influence of the difference of signal-to-noise ratios of different satellite frequency bands is ignored, so that the measurement precision is reduced.
Disclosure of Invention
The invention provides a Beidou and GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy, which is used for researching the influence of the selection length of a signal-to-noise ratio, the effective reflection high fluctuation of a satellite and the frequency characteristic difference of the satellite on the soil humidity measurement precision, effectively solving the problems of low measurement precision, unstable performance and the like of a single system, meeting the requirements of time and spatial resolution of soil humidity measurement, meeting the requirements of measurement precision, instantaneity, stability and the like and providing powerful support for soil humidity measurement.
In order to solve the technical problems, the technical scheme of the invention is as follows:
a Beidou and GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy comprises the following steps:
s1, establishing a phase reference value database of signal-to-noise ratio of a reflection signal
Continuously acquiring original observation values of Beidou and GPS without rainy days, wherein the original observation values comprise multi-frequency signal-to-noise ratios, pseudo ranges and carrier phases, calculating the altitude angle of a corresponding epoch moment, processing data of the original observation values, calculating to obtain signal-to-noise ratio phase reference values of the Beidou and the GPS, and storing the signal-to-noise ratio phase reference values as a database;
s2, calculating the current soil humidity of each signal-to-noise phase reference value by using a soil humidity measurement model and combining the signal-to-noise phase reference values in a signal-to-noise phase reference value database of the Beidou, and finally calculating the average value of the soil humidity obtained by calculating all frequencies of all the Beidou satellites to obtain the soil humidity measured by the Beidou system;
s3, calculating the current soil humidity of each signal-to-noise phase reference value by using a soil humidity measurement model and combining the signal-to-noise phase reference values in a signal-to-noise phase reference value database of the GPS, and finally, calculating the average value of the soil humidity obtained by calculating all frequencies of the GPS dual-frequency satellite to obtain the soil humidity measured by the GPS system;
and S4, combining the soil humidity measured by the Beidou and GPS dual system in a multi-frequency mode, averaging the results, and outputting the measurement result.
Preferably, the method for acquiring the signal-to-noise ratio phase reference value in step S1 includes:
s1-1, processing original data
S1-1-1, filtering the influence of a signal-to-noise ratio direct signal of an original observation value through a high-order Chebyshev polynomial fitting algorithm, and reserving a signal-to-noise ratio signal only containing a reflected signal and random noise;
s1-1-2, carrying out noise reduction processing on the signal-to-noise ratio signal processed in the step S1-1-1 through a Gaussian low-pass filtering noise reduction algorithm to obtain a signal-to-noise ratio signal only retaining a reflection signal;
s1-1-3, calculating the effective reflection height of the satellite through a self-adaptive sliding window and a weighting constraint strategy;
s1-2, calculating a signal-to-noise phase reference value through an extended Kalman filtering algorithm according to the satellite effective reflection height obtained in the step S1-1-3 and by combining the signal-to-noise ratio signal of the reflection signal obtained in the step S1-1-2.
Preferably, the principle of the high-order Chebyshev polynomial fitting algorithm in the step S1-1-1 is as follows:
at a time period t 0 ,t 0 +Δt]The signal-to-noise ratio signal in (c) is fitted with a Chebyshev polynomial, where t 0 Is the starting epoch of the signal-to-noise ratio, and at is the signal-to-noise ratio selected length,
to normalize the signal-to-noise ratio, the time is converted as follows:
at this time, the parameter τ ∈ -1,1], and thus, the signal-to-noise ratio can be expressed as a chebyshev polynomial:
in the formula, S d Expressing the signal-to-noise ratio of the fitted direct signal, n is the order of the Chebyshev polynomial, C i Chebyshev polynomial coefficient, chebyshev polynomial T, being the signal-to-noise ratio component i The recurrence formula is:
after fitting, the direct signal can be eliminated to obtain a signal only retaining the signal-to-noise ratio of the reflected signal and the high-frequency random noise, and the direct signal-to-noise ratio elimination formula is expressed as follows:
S r_n ＝S o -S d
in the formula: s r_n For signals containing only the signal-to-noise ratio and noise of the reflected signal, S o Representing the signal-to-noise ratio, S, just received d The signal-to-noise ratio of the direct signal of the corresponding point is obtained after Chebyshev polynomial fitting.
Preferably, the gaussian low-pass filtering noise reduction algorithm processing in the step S1-1-2 is expressed as follows:
suppose that the signal to be processed containing the SNR and noise of the reflected signal is denoted S r_n (x) And the signal after filtering which only contains the signal-to-noise ratio of the reflected signal is expressed as S r The calculation procedure is as follows:
S r (x)＝S r_n (x)*G(x)
in the formula, G (x) is a Gaussian function, the symbol represents convolution operation, x represents a signal-to-noise ratio sequence, high-frequency random noise can be effectively filtered through the processing of the process, a signal only retaining the signal-to-noise ratio of a reflection signal is obtained,
in the formula, G (x) is a Gaussian kernel function, the width of the Gaussian function is determined by sigma, the sigma can be calculated by the standard deviation of signal-to-noise ratio, the sigma is used as the distribution parameter of the Gaussian function, and e and pi are natural constants respectively.
Preferably, in step S1-1-3, the method for acquiring the satellite effective reflection is as follows:
(1) Determination of the signal-to-noise ratio signal S by means of an adaptive sliding algorithm r The data set is selected from a set of data,
the selection method comprises the following steps: the starting angle of the altitude angle is still 5 degrees, but the ending angle is selected in a sliding way between 25 and 35 degrees, a fixed value is not adopted, but the data groups are selected in a sliding way in a self-adaptive way to be respectively calculated in the middle, the selection rule is determined by the data sampling rate,
if the data sampling rate is 30s, the signal-to-noise ratio between 25 degrees and 35 degrees needs to be calculated in groups, when 20 epochs exist between 25 degrees and 35 degrees, the signal is divided into 20 groups, the first group of data consists of epochs between 5 degrees and 25 degrees of altitude plus the first epoch between 25 degrees and 35 degrees, the second group of data consists of epochs between 5 degrees and 25 degrees of altitude plus the first two epochs between 25 degrees and 35 degrees, and the like, so that 20 groups of data to be processed are formed;
if the data sampling rate is 15s, the step length is determined to be 2, namely, a group of data is selected for calculation every other epoch;
if the data sampling rate is 5s, the step length is 3, that is, a group of data is selected for calculation every three epochs. If the time is 1s, the step length is 5, a group of data is selected for calculation every five epochs,
after the adaptive sliding algorithm processing, the signal-to-noise ratio signals to be processed are divided into m groups, which can be expressed as:
that is, m groups of signals with different signal-to-noise ratio lengths can be used for calculating the effective reflection height of the satellite;
(2) And obtaining the effective reflection height of the satellite by adopting a Lomb-Scargle spectrum analysis method for each group of signal-to-noise ratio signals.
Wherein, lomb-Sacragle spectrum analysis method finds the power spectrum of the signal, and the expression is as follows:
in the formula, P x (f) Is the power of a periodic signal of frequency f; s. the r (x) For signal-to-noise ratio including the reflected signal, N is the sequence length, and τ is the time shift invariant, which can be calculated by the following equation:
can determine P x (f) F value when reaching peak value, and then calculating the effective reflection height h of the corresponding satellite according to the relationship between the frequency and the effective reflection height e The conversion relation is:
in the formula, λ is the wavelength of the signal and is a known quantity, and therefore, the effective reflection height h of the satellite can be obtained by the above formula e ，
Repeating the above calculation process to obtain m groups of signal-to-noise ratiosThe ratio signals are respectively resolved to obtain m effective reflection high results, which are expressed as:
averaging the calculated satellite effective reflection heights to obtain the satellite effective reflection height of the day
Wherein i represents the day
(3) Repeating the steps (1) and (2), resolving data of ten days before the measurement day, and obtaining the satellite effective reflection height of nearly ten days, which is expressed as
Weighting the satellite effective reflection height of the ten days by adopting a weighting constraint strategy, wherein a specific weighting model is expressed as follows:
in the formula, w i Representing a weighting coefficient, i represents a number of days, data for a total of ten days, data for the day of the measurement day is i =10, and so on, with the weight decreasing gradually.
Preferably, in step S1-2, the method for acquiring the snr phase reference value includes:
effective reflection height h by satellite e Combining the signal-to-noise ratio signal S containing the reflected signal r The phase value of the signal-to-noise ratio can be obtained by the relationship between the two, and the specific formula is as follows:
where A and φ represent the amplitude and phase of the signal-to-noise ratio, and are unknowns in the formula. E is the satellite altitude, a known quantity.
The solution is performed by using an extended kalman filter algorithm, as follows:
first, for a nonlinear system, the state equation and observation equation of its state space can be expressed as:
in the formula, the first term is a state equation, the second term is an observation equation, and x k And z k The actual state vector and the observation vector, respectively. Wherein x is k Is an unknown state quantity consisting of A and phi unknowns, which can be expressed as x k (A，φ)，z k Is formed by S r The observed state quantity, i.e. the output quantity of the system, of known quantity composition is denoted z k (S r )，u k As a function of the system, as an input to the equation of state, from a known quantity h e And E, k denotes the signal sequence, f is a nonlinear function of the system state at the time k-1 and k, h is the state quantity x k And the observed quantity z k Related non-linear function, w k And v k Respectively process noise and observation noise, and is Gaussian white noise with mutually independent mean value of zero,
in practical applications, w cannot be determined k And v k The specific value at each time instant, but the noise value can be ignored to obtain the estimated values of the state vector and the observed quantity at the k time instant, and therefore, the above equation can be approximately written as:
the extended kalman filter algorithm transforms a nonlinear system into a linear system by performing taylor first-order expansion of the nonlinear system near the optimal estimation point of its state. In this case, the above formula can be expressed as:
therefore, the snr phase reference value Φ and the amplitude value a can be obtained by solving. In this case only the phase value phi is used.
Preferably, the flow of the taylor first-order expansion algorithm is as follows:
1) Inputting initial conditions, inputting
And P k-1 ，
Is an a posteriori state estimate, P, of the system at time k-1 k-1 Is the error covariance of the system at time k-1;
2) Parameter prediction, the state prediction equation is:
wherein,
the posterior state estimation of the k-time system is carried out, and the error covariance prediction equation is as follows:
is the state at time k
And (3) predicting covariance, wherein A is a Jacobian matrix obtained by the partial derivation of x by the f function, W is a Jacobian matrix obtained by the partial derivation of W by f, and Q is a process noise covariance matrix.
A correction module, the system gain equation expressed as:
wherein, H is a Jacobian matrix after H is used for solving the partial derivative of x, V is a Jacobian matrix after H is used for solving the partial derivative of V, R is an observation noise covariance matrix, and the filtering equation is expressed as:
the meaning of the parameters in the formula is introduced before, and the error covariance update equation is expressed as:
wherein, I is a unit matrix and belongs to a linear algebra common sense matrix.
Preferably, the specific method of step S2 is as follows:
s2-1, reading a signal-to-noise ratio phase reference value of the Beidou system from a signal-to-noise ratio phase reference value database, and carrying out combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
wherein, delta phi is the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by preprocessing an original observation value obtained by a Beidou system on the day of measurement through the step S1;
s2-2, respectively calculating corresponding soil humidity by using data of three frequencies of all Beidou observation satellites, and then calculating an average value to obtain the final soil humidity measured by the Beidou system, wherein the final soil humidity is expressed as VWC B 。
Further, three frequencies of the Beidou satellite are respectively as follows: b1 is 1575.42MHz, B2 is 1176.45MHz, B3 is 1268.52MHz, and belongs to the common knowledge in the field.
Preferably, the specific method of step S3 is as follows:
s3-1, reading a signal-to-noise ratio phase reference value of a GPS system from a signal-to-noise ratio phase reference value database, and performing combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
in the formula, delta phi is the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by preprocessing an original observation value acquired by a GPS on the day of measurement through the step S1;
s3-2, respectively calculating corresponding soil humidity by using the dual-frequency data of all the GPS observation satellites, and then calculating the mean value to obtain the final soil humidity measured by the GPS system, wherein the final soil humidity is expressed as VWC G 。
It should be noted that, currently, all satellites of the GPS can transmit dual-frequency signals, and only some satellites transmit triple-frequency signals. Therefore, in order to maintain uniformity, dual frequency signals are used for measurement. The dual-frequency signal frequencies of the GPS satellite are respectively as follows: l1 is 1575.42MHz, and L2 is 1227.60MHz, which belongs to the common knowledge in the field.
Preferably, the implementation method of step S4 is:
and (5) solving the soil humidity of the final measurement day by utilizing the Beidou and GPS soil humidity obtained by calculation in the step (S3) and the step (S4), namely taking the mean value of the Beidou and GPS soil humidity, and expressing as follows:
in the formula, VWC B Is the soil humidity, VWC, measured by the Beidou system G Is the soil humidity, VWC, measured by the GPS system Final (a Chinese character of 'gan') I.e. the final measured soil moisture.
The invention has the following characteristics and beneficial effects:
by adopting the technical scheme, the existing Beidou/GPS navigation system satellite is utilized to measure the soil humidity, and compared with the existing method based on a remote sensing satellite measurement method and a probe measurement method, the problems of high measurement cost, high environmental influence and the like can be solved, and the problems of low measurement time and spatial resolution, incapability of meeting large-scale measurement and the like can be solved. Compared with the existing method based on the signal-to-noise ratio of the GPS signal, the method can solve the problems that the high precision of effective reflection of the satellite is influenced by the selection length of the signal-to-noise ratio, and the like, and can also solve the problems that the effective reflection of the satellite is high in fluctuation, the measurement precision is low and the measurement is unstable due to neglecting the frequency characteristic difference. Meanwhile, the Beidou and the GPS are combined to carry out dual-system multi-frequency combined measurement, so that the accuracy and the stability of the measurement can be guaranteed, and the soil humidity measurement with high precision, high real-time performance and high stability is met.
In addition, the signal-to-noise ratio of the direct signal is eliminated by using a high-order Chebyshev polynomial fitting method, denoising is carried out by adopting a Gaussian low-pass filtering denoising algorithm, and high-frequency random noise is filtered out, so that a signal only containing the signal-to-noise ratio of the reflected signal is obtained. Because the signal-to-noise ratio selection length directly influences the precision of the effective reflection height of the follow-up satellite, the effective reflection height of each frequency of each satellite is independently calculated by adopting a self-adaptive sliding algorithm, and the effective reflection heights of the satellite for continuous ten days are integrated by adopting a weighting constraint strategy, so that the precision of the effective reflection height of the satellite is effectively improved. And the problem of nonlinear solution is solved by adopting an extended Kalman filtering algorithm, and the solution accuracy of the phase value of the signal-to-noise ratio of the reflection signal is improved. The difference between different frequency characteristics of different satellites is considered, and the accuracy of soil humidity measurement is guaranteed by adopting a method of independently calculating all frequencies of all satellites and then calculating the mean value. The Beidou and GPS navigation system are combined to carry out dual-system multi-frequency combined detection, so that the precision, the real-time performance and the stability of soil humidity measurement are further ensured.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the embodiments or the description of the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the drawings without creative efforts.
Fig. 1 is a flow chart for establishing a snr-phase reference value database based on a sliding algorithm and a weighting strategy according to the present embodiment.
Fig. 2 is a flowchart of the Beidou/GPS soil humidity measurement method based on the sliding algorithm and the weighting strategy in the embodiment.
Detailed Description
It should be noted that the embodiments and features of the embodiments may be combined with each other without conflict.
In the description of the present invention, it is to be understood that the terms "central," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," "outer," and the like are used in the orientations and positional relationships indicated in the drawings, which are based on the orientations and positional relationships indicated in the drawings, and are used for convenience in describing the present invention and for simplicity in description, but do not indicate or imply that the device or element so referred to must have a particular orientation, be constructed in a particular orientation, and be operated, and thus should not be construed as limiting the present invention. Furthermore, the terms "first", "second", etc. are used for descriptive purposes only and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defined as "first," "second," etc. may explicitly or implicitly include one or more of that feature. In the description of the present invention, "a plurality" means two or more unless otherwise specified.
In the description of the present invention, it should be noted that, unless otherwise explicitly specified or limited, the terms "mounted," "connected," and "connected" are to be construed broadly, e.g., as meaning either a fixed connection, a removable connection, or an integral connection; can be mechanically or electrically connected; they may be connected directly or indirectly through intervening media, or they may be interconnected between two elements. The specific meanings of the above terms in the present invention can be understood by those of ordinary skill in the art through specific situations.
The invention provides a Beidou and GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy, which comprises the following steps as shown in figure 1:
s1, establishing a phase reference value database of signal-to-noise ratio of a reflection signal
Continuously acquiring original observation values of Beidou and GPS without rainy days, wherein the original observation values comprise multi-frequency signal-to-noise ratios, pseudo ranges and carrier phases, calculating altitude angles of corresponding epoch moments, processing data of the original observation values, calculating to obtain signal-to-noise ratio phase reference values of the Beidou and the GPS, and storing the signal-to-noise ratio phase reference values as a database.
It should be noted that the signal-to-noise ratio may be directly extracted from the observation file, and the altitude information corresponding to the epoch may be calculated by the pseudorange, the carrier phase, and the satellite ephemeris. This process is common knowledge in the field of positioning and the calculation process is not listed in detail here. For convenience of subsequent presentation, the raw SNR extracted here is used as S in the subsequent step o The calculated elevation angle is denoted by E.
Specifically, as shown in fig. 2, the method includes the following sub-steps:
s1-1, processing original data
S1-1-1, filtering out the influence of a signal-to-noise ratio direct signal of an original observation value through a high-order Chebyshev polynomial fitting algorithm, and reserving the signal-to-noise ratio signal only containing a reflected signal and random noise. The chebyshev polynomial fit is made by minimizing the sum of variances between the function value at a given point and the known function value at that point using the principle of the best approximation of the function, and the function is formed by fitting a chebyshev polynomial to the basis function, whose chebyshev polynomial fit algorithm principle is expressed as follows:
at a time period t 0 ,t 0 +Δt]The signal-to-noise ratio signal is fitted with a Chebyshev polynomial, where t 0 Is the starting epoch of the snr, and at is the snr selected length,
to normalize the signal-to-noise ratio, the time is converted as follows:
at this time, the parameter τ ∈ -1,1], and thus, the signal-to-noise ratio can be expressed as a chebyshev polynomial:
in the formula, S d The fitted direct signal to noise ratio is represented, n is the order of the Chebyshev polynomial, in order to improve the fitting accuracy of the signal to noise ratio, in the embodiment, a 5-order polynomial is selected for fitting, namely n =5, so that the fitting accuracy of the signal to noise ratio is improved. C i Chebyshev polynomial coefficient, chebyshev polynomial T, being a signal-to-noise ratio component i The recurrence formula is:
after fitting, the direct signal can be eliminated to obtain a signal only retaining the signal-to-noise ratio of the reflected signal and the high-frequency random noise, and the direct signal-to-noise ratio elimination formula is expressed as follows:
S r_n ＝S o -S d
in the formula: s r_n For signals containing only the signal-to-noise ratio and noise of the reflected signal, S o Representing the signal-to-noise ratio, S, just received d The signal-to-noise ratio of the direct signal of the corresponding point is obtained after Chebyshev polynomial fitting.
S1-1-2, noise reduction processing is carried out on the signal-to-noise ratio signal processed in the step S1-1-1 through a Gaussian low-pass filtering noise reduction algorithm, and a signal-to-noise ratio signal only retaining a reflection signal is obtained. Since the high-frequency random noise contained in the signal-to-noise ratio obeys normal distribution, a gaussian low-pass filter is adopted in the embodiment for denoising, so that the denoising effect is optimal. The gaussian low-pass filtering algorithm is essentially a linear smooth filtering algorithm, and mainly selects a weight value by using the shape of a gaussian function, so that not only can noise be accurately filtered, but also useful signals can be completely reserved. Since the snr signal in this embodiment is one-dimensional, a one-dimensional gaussian function with zero mean is used for filtering, which is expressed as follows:
suppose that the signal to be processed containing the signal-to-noise ratio and noise of the reflected signal is denoted S r_n (x) Only including inverse after filteringThe signal of the signal-to-noise ratio of the transmitted signal is denoted S r The calculation procedure is as follows:
S r (x)＝S r_n (x)*G(x)
in the formula, G (x) is a Gaussian function, the symbol represents convolution operation, x represents a signal-to-noise ratio sequence, high-frequency random noise can be effectively filtered through the processing of the process, a signal only retaining the signal-to-noise ratio of a reflection signal is obtained,
in the formula, G (x) is a Gaussian kernel function, the width of the Gaussian function is determined by sigma, the sigma can be calculated by the standard deviation of signal-to-noise ratio, the sigma is used as the distribution parameter of the Gaussian function, and e and pi are natural constants respectively.
By the technical scheme, high-frequency random noise can be effectively filtered to obtain a signal only retaining the signal-to-noise ratio of the reflected signal, and S is used in the subsequent steps r And (4) showing.
S1-1-3, calculating the effective reflection height of the satellite through an adaptive sliding window and a weighting constraint strategy.
The process can be divided into three steps, wherein the first step is to determine the signal-to-noise ratio selection length, and the second step is to solve the satellite effective reflection height by a Lomb-Scargle spectrum analysis method. And thirdly, obtaining the final effective reflection height of the satellite by adopting a weighting constraint strategy. The specific process is as follows:
(1) Determination of the signal-to-noise ratio signal S by means of an adaptive sliding algorithm r The selected data set.
It should be noted that, in the conventional method, the snr is selected to solve at an altitude angle of 5 to 30 degrees. However, the snr characteristics vary greatly from satellite to satellite and from frequency to frequency, resulting in a difference in effective reflection for satellites calculated at different snr lengths. Therefore, an adaptive sliding algorithm is used for the calculation. The specific process is as follows: the starting angle of the elevation angle is still 5 degrees, but the ending angle is between 25 and 35 degrees, and is not calculated with a fixed value, but with an adaptively chosen value between them. The selection rule is determined by the data sampling rate.
The specific selection method comprises the following steps: the starting angle of the elevation angle is still 5 degrees, but the ending angle is chosen in a sliding manner between 25 and 35 degrees, and is not calculated by using a fixed value, but by using adaptive sliding data sets. The selection rule is determined by the data sampling rate. Such as: if the data sampling rate is 30s, the snr of 25-35 degrees is calculated in groups (i.e. assuming 20 epochs between 25-35 altitudes, the signal is divided into 20 groups, the first group of data is composed of epochs at 5-25 degrees in altitude plus the first epoch at 25-35 degrees, the second group of data is composed of epochs at 5-25 degrees in altitude plus the first two epochs at 25-35 degrees, and so on to form 20 groups of data to be processed). If the time is 15s, the step length is determined to be 2, i.e. a group of data is selected for calculation every other epoch. If it is 5s, the step length is 3, i.e. a set of data is selected for calculation every three epochs. If the time is 1s, the step length is 5, and a group of data is selected for calculation every five epochs. After the adaptive sliding algorithm processing, it is assumed that the snr signals to be processed are divided into m groups, which can be expressed as:
i.e. m sets of signals of different snr lengths can be used to calculate the effective reflection height of the satellite.
(2) And (3) solving the satellite effective reflection height of each group of signal-to-noise ratio signals by adopting a Lomb-Scargle spectrum analysis method, firstly calculating the power spectrum of the signals, and finding out the corresponding frequency f when the energy is highest. And solving the corresponding satellite effective reflection height according to the relationship between the f value and the satellite effective reflection height.
Wherein, lomb-Sacragle spectrum analysis method is used for solving the power spectrum of the signal, and the expression is as follows:
in the formula, P x (f) Is the power of a periodic signal of frequency f; s. the r (x) For the signal-to-noise ratio of the included reflected signal, N is the sequence length, and τ isThe time shift invariant can be calculated by the following formula:
can determine P x (f) F value when reaching peak value, and then calculating the effective reflection height h of the corresponding satellite according to the relationship between the frequency and the effective reflection height e The conversion relation is:
where λ is the wavelength of the signal and is a known quantity. Therefore, the effective reflection height h of the satellite can be obtained by the above formula e ，
It should be noted that λ is a known quantity, and can be directly obtained, and the wavelength of the GPS L1 frequency band is: 19.0cm, the wavelength of big dipper B1 frequency channel is 19.2cm, belongs to common general knowledge in the field. Therefore, the effective reflection height of the satellite can be obtained by the above equation.
Repeating the calculation process, and respectively resolving m groups of signal-to-noise ratio signals to obtain m effective reflection high results, which are expressed as:
calculating the average value of the calculated satellite effective reflection heights to finally obtain the satellite effective reflection heights of the day
Wherein i represents the day;
(3) Repeating the steps (1) and (2), resolving data of ten days before the measurement day, and obtaining the satellite effective reflection height of nearly ten days, which is expressed as
Weighting the satellite effective reflection height of the ten days by adopting a weighting constraint strategy, wherein a specific weighting model is expressed as follows:
in the formula, w i Representing a weighting coefficient, i represents a number of days, data of ten days in total, data of the day of the measurement day is i =10, and the like, and the weight is gradually reduced.
It will be appreciated that the weighting patterns are weighted according to the proximity of time, i.e. the closer to the time of the measurement day, the higher the weight, the further away the time the lower the weight.
S1-2, calculating a signal-to-noise phase reference value through an extended Kalman filtering algorithm according to the satellite effective reflection height obtained in the step S1-1-3 and by combining the signal-to-noise ratio signal of the reflection signal obtained in the step S1-1-2.
Effective reflection height h by satellite e Combining the signal-to-noise ratio signal S containing the reflected signal r The phase value of the signal-to-noise ratio can be obtained by the relationship between the two, and the specific formula is as follows:
where A and φ represent the amplitude and phase of the signal-to-noise ratio, and are unknowns in the formula.
Furthermore, considering that the signal-to-noise ratio and the altitude angle are in a nonlinear relation at the moment, in order to improve the precision of parameter solution, an extended kalman filter algorithm is adopted for solution. The extended kalman filter algorithm has ready-made function modules in Matlab platform to call, and we only make a simple algorithm description here, as shown below:
first, for a nonlinear system, the state equation and observation equation of its state space can be expressed as:
in the formula, the first term is a state equation, the second term is an observation equation, and x k And z k Are respectively in actual stateVectors and observation vectors. Wherein x is k Is an unknown state quantity composed of A and phi unknowns, which can be expressed as x k (A，φ)；z k Is formed by S r The observed state quantity, i.e. the output quantity of the system, of known quantity composition is denoted z k (S r )；u k As a system function, as an input to the equation of state, from a known quantity h e And E. k denotes the signal sequence, f is a non-linear function of the system state at the time k-1 and k, h is the state quantity x k And the observed quantity z k Related non-linear function, w k And v k The noise is process noise and observation noise respectively, and is Gaussian white noise with mutually independent mean value of zero.
In practical applications, w cannot be determined k And v k The specific value at each time instant, but the noise value can be ignored to obtain the estimated values of the state vector and the observed quantity at the k time instant, and therefore, the above equation can be approximately written as:
the extended kalman filter algorithm transforms the nonlinear system into a linear system by performing a taylor first order expansion of the nonlinear system near the optimal estimation point of its state. In this case, the above formula can be expressed as:
therefore, the signal-to-noise ratio phase reference value phi can be obtained by solving.
Wherein, the flow of the Taylor first-order expansion algorithm is as follows:
1) Inputting initial conditions, inputting
And P k-1 ，
Is an a posteriori state estimate, P, of the system at time k-1 k-1 Is the error covariance of the system at time k-1;
2) Parameter prediction, the state prediction equation is:
wherein,
the posterior state estimation of the system at the moment k is carried out, and an error covariance prediction equation is as follows:
is the state at time k
And (3) predicting covariance, wherein A is a Jacobian matrix obtained by the partial derivation of x by the f function, W is a Jacobian matrix obtained by the partial derivation of W by f, and Q is a process noise covariance matrix.
3) A correction module, the system gain equation expressed as:
wherein, H is a Jacobian matrix after H is used for solving the partial derivative of x, V is a Jacobian matrix after H is used for solving the partial derivative of V, R is an observation noise covariance matrix, and the filtering equation is expressed as:
the meaning of the parameters in the formula is introduced before, and the error covariance update equation is expressed as:
wherein, I is a unit matrix and belongs to a linear algebra common sense matrix.
It should be noted that, in this embodiment, the signal-to-noise ratio phase reference value is obtained through extended kalman filtering, and may be obtained through a least square algorithm.
It can be understood that, in the above technical solution, the snr phase reference value database is established. By continuously acquiring navigation data of the Beidou and GPS systems in non-rainy days and by utilizing the process of the step S1, the sub-system calculates corresponding signal-to-noise ratio phase values according to the sub-frequency of the satellite, establishes a database according to the system type, the satellite number and the frequency number and stores the database on a local computer in a text form. For example, the phase value of the B1 frequency band of the Beidou B01 satellite is phi B01_1 And the phase value of the B2 frequency band is phi B01_2 . By analogy, a signal-to-noise ratio phase reference value database of multi-frequency of the Beidou and GPS dual system can be established.
S2, calculating the current soil humidity of each signal-to-noise phase reference value by using a soil humidity measurement model and combining the signal-to-noise phase reference value database of the Beidou, and finally calculating the average value of the soil humidity obtained by calculating all frequencies of all satellites of the Beidou to obtain the soil humidity measured by the Beidou system.
S2-1, reading a signal-to-noise ratio phase reference value of the Beidou system from a signal-to-noise ratio phase reference value database, and performing combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
in the formula, delta phi is the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by processing an original observation value obtained by a Beidou system on the day of measurement through the step S1;
it should be noted that the formula is an empirical formula for solving soil moisture by using signal-to-noise ratio, and belongs to the common knowledge in the field. Therefore, this embodiment will not be described in detail.
S2-2, respectively calculating corresponding soil humidity by using data of three frequencies of all Beidou observation satellites, and then calculating an average value to obtain the final soil humidity measured by the Beidou system, wherein the final soil humidity is expressed as VWC B 。
Further, the three frequencies of the Beidou satellite are respectively as follows: b1 is 1575.42MHz, B2 is 1176.45MHz, B3 is 1268.52MHz, and belongs to the common knowledge in the field.
S3, calculating the current soil humidity of each signal-to-noise ratio phase reference value by utilizing a soil humidity measurement model and combining with a signal-to-noise ratio phase reference value database of the GPS, and finally, calculating the average value of the soil humidity obtained by calculating all frequencies of the GPS dual-frequency satellite to obtain the soil humidity measured by the GPS system;
s3-1, reading a signal-to-noise ratio phase reference value of a GPS system from a signal-to-noise ratio phase reference value database, and performing combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
in the formula, Δ Φ is a difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by preprocessing an original observation value acquired by the GPS on the day of measurement through step S1.
S3-2, respectively calculating corresponding soil humidity by utilizing dual-frequency data of all observation satellites of the GPS, and then calculating an average value to obtain the final soil humidity measured by the GPS, wherein the final soil humidity is expressed as VWC G 。
It should be noted that, currently, all satellites of the GPS transmit dual-frequency signals, and only some satellites transmit triple-frequency signals. Therefore, to maintain uniformity, measurements are taken here using both frequency signals. The dual-frequency signal frequencies of the GPS satellite are respectively as follows: l1 is 1575.42MHz, and L2 is 1227.60MHz, which belongs to the common knowledge in the field.
And S4, combining the soil humidity measured by the Beidou and GPS dual-system multi-frequency measurement, averaging the results, and outputting the measurement result.
And (3) solving the soil humidity of the final measurement day by utilizing the Beidou obtained by calculation in the step (S3) and the step (S4) and the GPS soil humidity, namely taking the mean value of the Beidou and the GPS soil humidity, and expressing as follows:
in the formula, VWC B Is the soil humidity, VWC, measured by the Beidou system G Is the soil humidity, VWC, measured by the GPS system Terminal I.e. the final measured soil moisture.
The embodiments of the present invention have been described in detail with reference to the accompanying drawings, but the present invention is not limited to the described embodiments. It will be apparent to those skilled in the art that various changes, modifications, substitutions and alterations can be made in these embodiments, including the components, without departing from the principles and spirit of the invention, and still fall within the scope of the invention.
Claims (4)
1. A Beidou and GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy is characterized by comprising the following steps:
s1, establishing a phase reference value database of signal-to-noise ratio of a reflection signal
Continuously acquiring original observation values of Beidou and GPS without rainy days, wherein the original observation values comprise multi-frequency signal-to-noise ratios, pseudo ranges and carrier phases, calculating elevation angles of corresponding epoch moments, processing data of the original observation values, calculating to obtain signal-to-noise ratio phase reference values of the Beidou and the GPS, and storing the signal-to-noise ratio phase reference values as a database;
s1-1, processing original data
S1-1-1, filtering out the influence of a signal-to-noise ratio direct signal in an original observation value through a high-order Chebyshev polynomial fitting algorithm, and reserving the signal-to-noise ratio signal only containing a reflected signal and random noise;
the principle of the high-order Chebyshev polynomial fitting algorithm is expressed as follows:
in a time period t 0 ,t 0 +Δt]The signal-to-noise ratio signal is fitted with a Chebyshev polynomial, where t 0 Is the starting epoch of the snr, and at is the snr selected length,
to normalize the signal-to-noise ratio, the time is converted as follows:
at this time, the parameter τ ∈ [ -1,1], and therefore, the signal-to-noise ratio can be expressed as a Chebyshev polynomial:
in the formula, S d Representing the signal-to-noise ratio of the fitted direct signal, n being the order of the Chebyshev polynomial, C i Chebyshev polynomial coefficient, chebyshev polynomial T, being the signal-to-noise ratio component i The recurrence formula is:
after fitting, the direct signal can be eliminated to obtain a signal only retaining the signal-to-noise ratio of the reflected signal and the high-frequency random noise, and the direct signal-to-noise ratio elimination formula is expressed as follows:
S r_n ＝S o -S d
in the formula: s. the r_n For signals containing only the signal-to-noise ratio and noise of the reflected signal, S o Representing the original just-received signal-to-noise ratio, S d The signal-to-noise ratio of the direct signal of the corresponding point is obtained after Chebyshev polynomial fitting;
s1-1-2, noise reduction processing is carried out on the signal-to-noise ratio signal processed in the step S1-1-1 through a Gaussian low-pass filtering noise reduction algorithm to obtain a signal-to-noise ratio signal only retaining a reflection signal,
gaussian low-pass filtering noise reduction algorithm processing, expressed as follows:
suppose that the signal to be processed containing the signal-to-noise ratio and noise of the reflected signal is denoted S r_n (x) The signal after filtering, which contains only the signal-to-noise ratio of the reflected signal, is denoted as S r The calculation process is shown as follows：
S r (x)＝S r_n (x)*G(x)
In the formula, G (x) is a Gaussian function, the symbol x represents convolution operation, x represents a signal-to-noise ratio sequence, high-frequency random noise can be effectively filtered through the processing of the process, a signal only keeping the signal-to-noise ratio of the reflection signal is obtained,
in the formula, G (x) is a Gaussian kernel function, the width of the Gaussian function is determined by sigma, the sigma can be calculated by the standard deviation of signal-to-noise ratio and is used as the distribution parameter of the Gaussian function, and e and pi are natural constants respectively;
s1-1-3, calculating the effective reflection height of the satellite through an adaptive sliding window algorithm and a weighting constraint strategy,
the method for acquiring the effective reflection height of the satellite comprises the following steps:
(1) Determination of the signal-to-noise ratio signal S by means of an adaptive sliding algorithm r The data set is selected from a set of data,
the selection method comprises the following steps: the starting angle of the altitude angle is still 5 degrees, but the ending angle is selected in a sliding way between 25 and 35 degrees, a fixed value is not adopted, but the data groups are selected in a sliding way in a self-adaptive way to be respectively calculated in the middle, the selection rule is determined by the data sampling rate,
if the data sampling rate is 30s, the signal-to-noise ratio between 25 degrees and 35 degrees needs to be calculated in groups, when 20 epochs exist between 25 degrees and 35 degrees, the signal is divided into 20 groups, the first group of data consists of epochs between 5 degrees and 25 degrees of altitude plus the first epoch between 25 degrees and 35 degrees, the second group of data consists of epochs between 5 degrees and 25 degrees of altitude plus the first two epochs between 25 degrees and 35 degrees, and the like, so that 20 groups of data to be processed are formed;
if the data sampling rate is 15s, the step length is determined to be 2, namely, a group of data is selected for calculation every other epoch;
if the data sampling rate is 5s, the step length is 3, namely, one group of data is selected for calculation every three epochs, if the data sampling rate is 1s, the step length is 5, one group of data is selected for calculation every five epochs,
after the adaptive sliding algorithm processing, the signal-to-noise ratio signals to be processed are divided into m groups, which can be expressed as:
that is, m groups of signals with different signal-to-noise ratio lengths can be used for calculating the effective reflection height of the satellite;
(2) The effective reflection height of the satellite is obtained by adopting a Lomb-Scargle spectrum analysis method for each group of signal-to-noise ratio signals,
wherein, lomb-Sacragle spectrum analysis method is used for solving the power spectrum of the signal, and the expression is as follows:
in the formula, P x (f) Is the power of a periodic signal of frequency f; s r (x) The signal-to-noise ratio of the reflected signal, N the sequence length, and τ the time shift invariant, can be calculated by the following equation:
can determine P x (f) F value when reaching peak value, and then calculating the effective reflection height h of the corresponding satellite according to the relationship between the frequency and the effective reflection height e The conversion relation is:
in the formula, λ is the wavelength of the signal and is a known quantity, and therefore, the effective reflection height h of the satellite can be obtained by the above formula e ，
Repeating the calculation process, and respectively resolving m groups of signal-to-noise ratio signals to obtain m effective reflection high results, which are expressed as:
calculating the average value of the calculated satellite effective reflection heights to finally obtain the satellite effective reflection heights of the day
Wherein i represents the day;
(3) Repeating the steps (1) and (2), resolving data of ten days before the measurement day, and obtaining the effective reflection height of the satellite of nearly ten days, which is expressed as
And weighting the satellite effective reflection heights of the ten days by adopting a weighting constraint strategy, wherein a specific weighting model is represented as follows:
in the formula, w i Representing a weighting coefficient, i represents the number of days, the data of ten days in total are measured, the data of the day of the measurement day is i =10, and the like, and the weight is gradually reduced;
s1-2, calculating a signal-to-noise phase reference value through an extended Kalman filtering algorithm according to the satellite effective reflection height obtained in the step S1-1-3 and by combining the signal-to-noise ratio signal of the reflection signal obtained in the step S1-1-2,
the method for acquiring the signal-to-noise ratio phase reference value comprises the following steps:
effective reflection height h by satellite e Combining the signal-to-noise ratio signal S containing the reflected signal r The phase value of the signal-to-noise ratio can be obtained by the relationship between the two, and the specific formula is as follows:
where A and φ represent the amplitude and phase of the SNR, and are unknowns in the formula, and E is the height of the epoch timeAngle, a known quantity, obtained by preprocessing in step S1, satellite effective reflection height h e Obtained in step S1-1-2, including the signal-to-noise ratio signal S of the reflected signal r And is also a known quantity obtained by the step S1-1-3, therefore, the unknown quantity A and phi can be solved by an extended Kalman filtering algorithm,
the solution is performed by using an extended kalman filter algorithm, as follows:
first, for a nonlinear system, the state equation and observation equation of its state space can be expressed as:
in the formula, the first term is a state equation, the second term is an observation equation, and x k And z k Respectively, an actual state vector and an observation vector, where x k Is an unknown state quantity consisting of A and phi unknowns, which can be expressed as x k (A，φ)；z k Is formed by S r The observed state quantity, i.e. the output quantity of the system, of known quantity composition is denoted z k (S r )；u k As a function of the system, as an input to the equation of state, from a known quantity h e And E, k represents a signal sequence,
f is a non-linear function of the state of the system at time k-1 and k, and h is a state quantity x k And the observed quantity z k Related non-linear function, w k And v k Respectively process noise and observation noise, and is Gaussian white noise with mutually independent mean value of zero,
in practice, w cannot be determined k And v k The specific values at each time instant, but the noise values can be ignored to obtain the estimated values of the state vector and the observed quantity at the time instant k, so the above equation can be approximately written as:
the extended kalman filter algorithm transforms the nonlinear system into a linear system by performing taylor first order expansions of the nonlinear system near the optimal estimation point of its state, which can be expressed as:
therefore, the signal-to-noise ratio phase reference value phi can be obtained by solving;
s2, calculating the current soil humidity of each signal-to-noise phase reference value by using a soil humidity measurement model and combining the signal-to-noise phase reference values in a signal-to-noise phase reference value database of the Beidou, and finally, calculating the average value of the soil humidity obtained by calculating all frequencies of all the Beidou satellites to obtain the soil humidity measured by a Beidou system;
s3, calculating the current soil humidity of each signal-to-noise phase reference value by using a soil humidity measurement model and combining the signal-to-noise phase reference values in a signal-to-noise phase reference value database of the GPS, and finally, calculating the average value of the soil humidity obtained by calculating all frequencies of the GPS dual-frequency satellite to obtain the soil humidity measured by the GPS system;
and S4, combining the soil humidity measured by the Beidou and GPS dual system in a multi-frequency mode, averaging the results, and outputting the measurement result.
2. The Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy as claimed in claim 1 is characterized in that the specific method of the step S2 is as follows:
s2-1, reading a signal-to-noise ratio phase reference value of the Beidou system from a signal-to-noise ratio phase reference value database, and performing combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
in the formula, delta phi is the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by preprocessing an original observation value obtained by a Beidou system on the day of measurement through the step S1;
s2-2, respectively calculating corresponding soil humidity by using data of three frequencies of all Beidou observation satellites, and then solving an average value to obtain the final soil humidity measured by the Beidou system, wherein the final soil humidity is expressed as VWC B The three frequencies of the Beidou satellite are respectively as follows: b1 is 1575.42MHz, B2 is 1176.45MHz, and B3 is 1268.52MHz.
3. The Beidou and GPS soil moisture measurement method based on sliding algorithm and weighting strategy according to claim 2, characterized in that the specific method of the step S3 is as follows:
s3-1, reading a signal-to-noise ratio phase reference value of a GPS system from a signal-to-noise ratio phase reference value database, and performing combined calculation by combining a signal-to-noise ratio phase value of the day to obtain soil humidity of the measurement day, wherein the calculation formula is as follows:
in the formula, delta phi is the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the day of measurement, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the day of measurement is a value obtained by preprocessing an original observation value acquired by a GPS on the day of measurement through the step S1;
s3-2, respectively calculating corresponding soil humidity by using the dual-frequency data of all the GPS observation satellites, and then calculating the mean value to obtain the final soil humidity measured by the GPS system, wherein the final soil humidity is expressed as VWC G The dual-frequency signal frequencies of the GPS satellite are respectively as follows: l1 is 1575.42MHz and L2 is 1227.60MHz.
4. The Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy according to claim 3, characterized in that the implementation method of the step S4 is as follows:
and (5) solving the soil humidity of the final measurement day by utilizing the Beidou and GPS soil humidity obtained by calculation in the step (S3) and the step (S4), namely taking the mean value of the Beidou and GPS soil humidity, and expressing as follows:
in the formula, VWC B Is the soil humidity, VWC, measured by the Beidou system G Is the soil humidity, VWC, measured by the GPS system Final (a Chinese character of 'gan') I.e. the final measured soil moisture.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210770775.1A CN114839354B (en) | 2022-07-02 | 2022-07-02 | Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210770775.1A CN114839354B (en) | 2022-07-02 | 2022-07-02 | Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114839354A CN114839354A (en) | 2022-08-02 |
| CN114839354B true CN114839354B (en) | 2022-11-18 |
Family
ID=82573543
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210770775.1A Active CN114839354B (en) | 2022-07-02 | 2022-07-02 | Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114839354B (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107796484A (en) * | 2017-01-11 | 2018-03-13 | 中南大学 | One kind is based on BDStar navigation system signal-to-noise ratio data observed stage changing method |
Family Cites Families (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101900692B (en) * | 2010-06-18 | 2012-05-30 | 武汉大学 | Method for measuring large-area soil humidity |
| CN102968552B (en) * | 2012-10-26 | 2016-01-13 | 郑州威科姆科技股份有限公司 | A kind of satellite orbit data estimation and modification method |
| CN104020180B (en) * | 2014-06-19 | 2016-02-10 | 武汉大学 | Based on the soil moisture retrieval method of the low elevation signals that Big Dipper base station receives |
| CN104224181B (en) * | 2014-09-26 | 2016-04-20 | 中国科学院生物物理研究所 | A kind of SAR real-time monitoring system of multi-channel magnetic resonance imaging equipment and method |
| CN105352979B (en) * | 2015-12-08 | 2017-11-28 | 武汉大学 | Soil moisture method of estimation based on Big Dipper GEO satellite signal |
| CN105787474A (en) * | 2016-03-29 | 2016-07-20 | 武汉瑞能电力设备有限责任公司 | Processing method of bridge vibration monitoring data |
| CN106290408B (en) * | 2016-07-21 | 2019-01-04 | 清华大学 | Soil moisture measurement method based on the station continuous operation GNSS signal-to-noise ratio data |
| CN106125106A (en) * | 2016-08-10 | 2016-11-16 | 清华大学 | The method measuring soil moisture based on the ground Big Dipper/GPS dual-mode survey station |
| CN109932735B (en) * | 2019-03-25 | 2023-04-07 | 中国铁路设计集团有限公司 | Positioning method for resolving Beidou short baseline single-frequency single epoch |
| CN111474209A (en) * | 2020-03-13 | 2020-07-31 | 山东航向电子科技有限公司 | Real-time soil humidity measuring method based on intelligent terminal equipment |
| CN111399015B (en) * | 2020-04-07 | 2022-05-31 | 中船重工鹏力（南京）大气海洋信息***有限公司 | Dual-mode satellite fusion positioning method suitable for ship traffic management system |
| CN111896984B (en) * | 2020-07-10 | 2022-10-28 | 自然资源部第一海洋研究所 | GNSS-based real-time high-precision wave measurement method and device |
| CN212692782U (en) * | 2020-09-10 | 2021-03-12 | 安徽阡陌网络科技有限公司 | High-precision lightweight handheld land area measuring device |
| CN112505068B (en) * | 2020-11-03 | 2023-08-11 | 桂林理工大学 | GNSS-IR-based earth surface soil humidity multi-star combination inversion method |
| CN215953863U (en) * | 2020-12-25 | 2022-03-04 | 武汉大学 | Beidou/GNSS water vapor real-time inversion device |
| CN112782689A (en) * | 2020-12-29 | 2021-05-11 | 西南交通大学 | Multi-satellite data fusion GNSS-IR soil humidity monitoring method |
| CN113049777B (en) * | 2021-03-12 | 2022-02-08 | 北京航空航天大学 | Device for measuring soil humidity through GNSS direct reflection signal carrier interference |
| CN113959329A (en) * | 2021-12-17 | 2022-01-21 | 西南交通大学 | Snow depth inversion method based on multi-satellite data fusion |
| CN114355411B (en) * | 2021-12-22 | 2023-07-25 | 杭州电子科技大学 | Flood detection method based on Beidou or GPS carrier-to-noise ratio observation value |
| CN114397425B (en) * | 2021-12-22 | 2024-04-19 | 杭州电子科技大学 | GNSS-IR soil humidity inversion method based on generalized extension approximation |
-
2022
- 2022-07-02 CN CN202210770775.1A patent/CN114839354B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107796484A (en) * | 2017-01-11 | 2018-03-13 | 中南大学 | One kind is based on BDStar navigation system signal-to-noise ratio data observed stage changing method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114839354A (en) | 2022-08-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114518586B (en) | GNSS precise single-point positioning method based on spherical harmonic expansion | |
| CN111583214B (en) | Sea surface wind speed inversion method based on RBF neural network and based on marine radar image | |
| Vousdoukas et al. | Performance of intertidal topography video monitoring of a meso-tidal reflective beach in South Portugal | |
| CN109001845B (en) | rainfall forecasting method | |
| CN109522516B (en) | Soil humidity detection method and device based on random forest regression algorithm and electronic equipment | |
| CN111829957A (en) | System and method for inverting moisture content of winter wheat plants based on multispectral remote sensing of unmanned aerial vehicle | |
| CN109597040A (en) | A kind of field-free geometric calibration method of satellite-borne SAR image | |
| CN116519913B (en) | GNSS-R data soil moisture monitoring method based on fusion of satellite-borne and foundation platform | |
| CN114297938B (en) | Inversion method of optical shallow water depth based on neural network | |
| CN114924241A (en) | Frequency correction method and system for satellite-borne rainfall measurement radar and ground-based weather radar | |
| CN112782701A (en) | Visibility perception method, system and equipment based on radar | |
| CN111144350B (en) | Remote sensing image positioning accuracy evaluation method based on reference base map | |
| Levin et al. | Observation impacts on the Mid-Atlantic Bight front and cross-shelf transport in 4D-Var ocean state estimates: Part I—Multiplatform analysis | |
| CN114839354B (en) | Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy | |
| Yu et al. | Generating daily 100 m resolution land surface temperature estimates continentally using an unbiased spatiotemporal fusion approach | |
| CN115980317B (en) | Foundation GNSS-R data soil moisture estimation method based on corrected phase | |
| CN116029162B (en) | Flood disaster inundation range monitoring method and system by using satellite-borne GNSS-R data | |
| CN114880883B (en) | Mountain land surface soil moisture remote sensing estimation method and device and electronic equipment | |
| AU2021105536A4 (en) | A High Spatial-Temporal Resolution Method for Near-Surface Air Temperature Reconstruction | |
| CN115205698A (en) | Lake water volume change estimation method combining height measurement data and water level-volume curve | |
| CN111044489B (en) | Method for obtaining atmosphere refractive index height distribution profile based on multi-wavelength measurement | |
| CN113625307A (en) | Landslide monitoring system and method based on GNSS | |
| Danylenko et al. | Monitoring of soil moisture in the south of Ukraine using active and passive remote sensing data | |
| CN114814173B (en) | Dielectric constant-based satellite-borne GNSS-R soil humidity inversion method and system | |
| Mironov et al. | A Method for Continues Calibration of a Rotating Antenna Scatterometer in Application to CFOSAT Measurements |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |