CN114839354A - Beidou/GPS soil humidity measurement method based on sliding algorithm and weighting strategy - Google Patents

Beidou/GPS soil humidity measurement method based on sliding algorithm and weighting strategy Download PDF

Info

Publication number
CN114839354A
CN114839354A CN202210770775.1A CN202210770775A CN114839354A CN 114839354 A CN114839354 A CN 114839354A CN 202210770775 A CN202210770775 A CN 202210770775A CN 114839354 A CN114839354 A CN 114839354A
Authority
CN
China
Prior art keywords
signal
noise ratio
soil humidity
beidou
gps
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.)
Granted
Application number
CN202210770775.1A
Other languages
Chinese (zh)
Other versions
CN114839354B (en
Inventor
苏明坤
杨艳茜
乔磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN202210770775.1A priority Critical patent/CN114839354B/en
Publication of CN114839354A publication Critical patent/CN114839354A/en
Application granted granted Critical
Publication of CN114839354B publication Critical patent/CN114839354B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/24Earth materials
    • G01N33/246Earth materials for water content
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13Receivers
    • G01S19/33Multimode operation in different systems which transmit time stamped messages, e.g. GPS/GLONASS
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02Preprocessing
    • G06F2218/04Denoising
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A40/00Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
    • Y02A40/10Adaptation 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/GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy, which comprises the following steps: s1, establishing a reflected signal-to-noise ratio phase reference value database; s2, calculating the soil humidity measured by the Beidou 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 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 in a multi-frequency mode, 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 measurement precision, unstable performance and the like of a single system, can meet the requirements of soil humidity measurement time and spatial resolution, and can meet the requirements of soil humidity measurement with high precision, high real-time performance and high stability.

Description

Beidou/GPS soil humidity measurement method based on sliding algorithm and weighting strategy
Technical Field
The invention relates to the technical field of navigation satellite remote sensing inversion, in particular to a Beidou/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, a large area range can be measured, time and spatial resolution are high, an optical lens and a sensor which are depended by the remote sensing satellite are easily influenced by the environment, and measurement failure or accuracy reduction based on the remote sensing satellite method can be caused by severe weather such as cloud, fog, haze and the like. The third type is soil humidity measurement based on a GNSS (Global Navigation Satellite System) Navigation Satellite, and the method utilizes the phase fluctuation of the signal-to-noise ratio of the reflection signal of the GNSS Satellite to invert the soil humidity, thereby acquiring the soil humidity measurement result in a large range 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 amount accuracy based on the signal-to-noise ratio, but adopts the fixed signal-to-noise ratio of 5 to 30 degrees to select the length, so that the final soil humidity measurement accuracy is reduced. Meanwhile, the influence of the high fluctuation of the effective reflection of the satellite along with the 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/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 space 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/GPS soil humidity measurement method based on a sliding algorithm and a weighting strategy comprises the following steps:
s1, establishing a database of the phase reference value of the signal-to-noise ratio of the reflected signal
Connecting and collecting original observation values of the Beidou and the GPS on 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;
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 the soil humidity measurement model and combining the signal-to-noise phase reference values in the 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 multi-frequency measurement of the Beidou and GPS dual systems, averaging the results, and outputting the measurement result.
Preferably, the method for obtaining the signal-to-noise ratio phase reference value in step S1 includes:
s1-1, processing the original data
S1-1-1, filtering the influence of the signal-to-noise ratio direct signal of the original observation value through a high-order Chebyshev polynomial fitting algorithm, and reserving the signal-to-noise ratio signal only containing the 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 an 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:
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:
Figure 72905DEST_PATH_IMAGE001
at this time, the parameter τ ϵ [ -1, 1], and thus, the signal-to-noise ratio can be expressed as a Chebyshev polynomial:
Figure 865281DEST_PATH_IMAGE002
in the formula,
Figure 471842DEST_PATH_IMAGE003
the signal-to-noise ratio of the fitted direct signal is represented,nfor the order of the chebyshev polynomial,C i chebyshev polynomial coefficients, chebyshev polynomials, being signal-to-noise ratio componentsT i The recurrence formula is:
Figure 327191DEST_PATH_IMAGE004
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:
Figure 550361DEST_PATH_IMAGE005
in the formula:S r_n for signals that contain only the signal-to-noise ratio and noise of the reflected signal,S o representing the signal-to-noise ratio that was originally just received,S 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:
assuming that the signal to be processed containing the signal-to-noise ratio and noise of the reflected signal is represented asS r_n (x) And the signal after filtering which only contains the signal-to-noise ratio of the reflected signal is expressed asS r The calculation procedure is as follows:
Figure 552953DEST_PATH_IMAGE006
in the formula,G(x) Is a gaussian function, the symbol denotes the convolution operation,xrepresenting the SNR sequence, the high-frequency random noise can be effectively filtered through the processing of the process, the signal only retaining the SNR of the reflected signal is obtained,
Figure 404234DEST_PATH_IMAGE007
wherein G (x) is a Gaussian kernel function having a width of
Figure 99657DEST_PATH_IMAGE008
It is decided that,
Figure 544545DEST_PATH_IMAGE008
can be calculated from the standard deviation of the signal-to-noise ratio, as a distribution parameter of the gaussian function,eandπare natural constants respectively.
Preferably, in step S1-1-3, the method for acquiring the satellite effective reflection is as follows:
(1) determination of signal-to-noise ratio signal by adaptive sliding algorithmS 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 altitudes, the signals are 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 every other epoch for calculation;
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 signal to be processed is divided intomGroup, can be represented as:
Figure 413144DEST_PATH_IMAGE009
That is to say havemThe signals with different signal-to-noise ratio lengths can be used for calculating the effective reflection height of the satellite;
(2) the Lomb-Scargle spectrum analysis method is adopted to solve the satellite effective reflection height 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:
Figure 322194DEST_PATH_IMAGE010
in the formula,P x ( f ) Is at a frequency off The power of the periodic signal of (a);S r (x) In order to contain the signal-to-noise ratio of the reflected signal,Nin order to be the length of the sequence,τthe time shift invariant can be calculated by the following formula:
Figure 860623DEST_PATH_IMAGE011
can find outP x ( f ) At the time of peakf The value is calculated according to the relation between the frequency and the effective reflection heighth e The conversion relation is:
Figure 183020DEST_PATH_IMAGE012
in the formula,λis the wavelength of the signal, and is a known quantity, therefore, the effective reflection height of the satellite can be obtained by the above formulah e
Repeating the above calculation processmRespectively resolving the group signal-to-noise ratio signals to obtainmThe effective reflection is high, expressed as:
Figure 527414DEST_PATH_IMAGE013
the satellite to be calculated isThe effective reflection height is averaged to finally obtain the effective reflection height of the satellite in the day
Figure 228653DEST_PATH_IMAGE014
WhereiniTo show 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
Figure 300DEST_PATH_IMAGE015
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:
Figure 809993DEST_PATH_IMAGE016
in the formula,w i which represents the weighting coefficient(s) of the,idata representing the number of days, ten days in total, measured on the day of the dayi=10, and so on, the weights decrease gradually.
Preferably, in step S1-2, the signal-to-noise ratio phase reference value obtaining method includes:
high effective reflection through satellite
Figure 895761DEST_PATH_IMAGE017
Combining the signal-to-noise ratio signal containing the reflected signal
Figure 513824DEST_PATH_IMAGE018
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:
Figure 784268DEST_PATH_IMAGE019
in the formula,Aand
Figure 691045DEST_PATH_IMAGE020
the amplitude and phase representing the signal-to-noise ratio are of the formulaAn unknown quantity.EIs 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:
Figure 377241DEST_PATH_IMAGE021
wherein the first term is a state equation, the second term is an observation equation,x k andz k the actual state vector and the observation vector, respectively. Wherein,x k is formed byAAnd
Figure 177707DEST_PATH_IMAGE020
the unknown state quantity of the unknown quantity composition can be expressed asx k A
Figure 291156DEST_PATH_IMAGE020
),z k Is formed by
Figure 685228DEST_PATH_IMAGE022
The observed state quantity, i.e. the output quantity of the system, of known quantity composition is expressed asz k
Figure 237432DEST_PATH_IMAGE022
),u k As a function of the system, as input quantities to the equation of state, from known quantities
Figure 830088DEST_PATH_IMAGE023
AndEthe components of the composition are as follows,kwhich is indicative of a sequence of signals,f to relate tok-1 andka non-linear function of the state of the system at the moment,his a state quantityx k And observed quantityz k The non-linear function of interest is,w k andv k respectively process noise and observationsNoise, and is white Gaussian noise with mutually independent mean value of zero,
in practical applications, it cannot be determinedw k Andv k the specific value at each time instant, but the estimate of the state vector and the observed quantity for which the noise value comes at time k can be ignored, and therefore the above equation can be approximately written as:
Figure 317701DEST_PATH_IMAGE024
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:
Figure 326633DEST_PATH_IMAGE025
therefore, the signal-to-noise ratio phase reference value can be obtained by solving
Figure 354632DEST_PATH_IMAGE020
Sum amplitude valueA. In this case only phase values
Figure 739477DEST_PATH_IMAGE020
Preferably, the flow of the taylor first-order expansion algorithm is as follows:
1) inputting initial conditions, inputting
Figure 788204DEST_PATH_IMAGE026
And
Figure 156869DEST_PATH_IMAGE027
Figure 988558DEST_PATH_IMAGE028
is thatk-1 a posteriori state estimates of the system at time instance,
Figure 352544DEST_PATH_IMAGE029
is thatk-1 time instant systematic error covariance;
2) parameter prediction, the state prediction equation is:
Figure 447538DEST_PATH_IMAGE030
wherein
Figure 100237DEST_PATH_IMAGE031
is thatkThe posterior state estimation of the time system, the error covariance prediction equation is:
Figure 532355DEST_PATH_IMAGE032
Figure 157371DEST_PATH_IMAGE033
is thatkState of the moment
Figure 220005DEST_PATH_IMAGE034
The covariance of the prediction is determined by the prediction,Ais thatfFunction pairxThe jacobian matrix after the partial derivation is solved,Wis thatfTo pairwThe jacobian matrix after the partial derivation is solved,Qis the process noise covariance matrix.
A correction module, the system gain equation expressed as:
Figure 156737DEST_PATH_IMAGE035
whereinHis thathTo pairxThe jacobian matrix after the partial derivation is solved,Vis thathTo pairvThe jacobian matrix after the partial derivation is solved,Ris an observed noise covariance matrix, where the filter equation is expressed as:
Figure 595809DEST_PATH_IMAGE036
in the formula, the meaning of the parameters is introduced before, and the error covariance update equation is expressed as:
Figure 809753DEST_PATH_IMAGE037
whereinIthe unit matrix belongs to a linear algebra common sense matrix.
Preferably, the specific method of step S2 is as follows:
s2-1, reading the signal-to-noise ratio phase reference value of the Beidou system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 371184DEST_PATH_IMAGE038
in the formula,
Figure 732895DEST_PATH_IMAGE039
the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the measurement day is obtained by preprocessing the original observation value obtained by the Beidou system on the measurement day through VWC, VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the measurement day is obtained through 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, and B3 is 1268.52MHz, which belongs to the common knowledge in the field.
Preferably, the specific method of step S3 is as follows:
s3-1, reading the signal-to-noise ratio phase reference value of the GPS system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 382182DEST_PATH_IMAGE040
in the formula,
Figure 106425DEST_PATH_IMAGE039
the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the measurement day is obtained by preprocessing the original observation value acquired by the GPS on the measurement day through the step S1, wherein the VWC represents the soil humidity required by the measurement day;
s3-2, respectively calculating corresponding soil humidity by using 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 can 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, L2 is 1227.60MHz, and is common knowledge in the field.
Preferably, the implementation method of step S4 is:
and (4) solving the soil humidity of the final measurement day by using the Beidou and GPS soil humidity obtained by calculation in the step S3 and the step S4, namely, taking the average value of the Beidou and GPS soil humidity, and expressing as follows:
Figure 510861DEST_PATH_IMAGE041
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 method can solve the problems of high measurement cost, high environmental influence and the like, and can solve the problems of low measurement time and spatial resolution, incapability of meeting large-scale measurement and the like. 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 using 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. And the Beidou and GPS navigation systems 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 description of the embodiments or the prior art will be briefly described below, and 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 these 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 "center", "longitudinal", "lateral", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", and the like, indicate orientations or positional relationships based on those shown in the drawings, and are used only for convenience in describing the present invention and for simplicity in description, and do not indicate or imply that the referenced devices or elements must have a particular orientation, be constructed and operated in a particular orientation, and thus, are not to 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 meaning 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/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 the signal-to-noise ratio of the reflected signal
The method comprises the steps of connecting and collecting original observation values of the Beidou and the GPS in 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.
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 here in detail. For the convenience of subsequent presentation, the raw SNR extracted here is used in subsequent stepsS o Indicating, for calculated altitude angleEAnd (4) showing.
Specifically, as shown in fig. 2, the method includes the following sub-steps:
s1-1, processing the original data
S1-1-1, filtering the influence of the signal-to-noise ratio direct signal of the original observation value through a high-order Chebyshev polynomial fitting algorithm, and reserving the signal-to-noise ratio signal only containing the reflected signal and random noise. The chebyshev polynomial fitting is based on the principle of the best approximation of a function, which is formed by fitting a chebyshev polynomial as a basis function, using the function as a minimum of the sum of the variances between the function value at a given point and the known function value at that point, and the principle of the chebyshev polynomial fitting algorithm is expressed as follows:
in the time periodt 0 , t 0t]The signal-to-noise ratio signal is fitted with a Chebyshev polynomial, wheret 0 Is the starting epoch of the signal-to-noise ratio, and ΔtThe length is chosen for the signal-to-noise ratio,
to normalize the signal-to-noise ratio, the time is converted as follows:
Figure 297552DEST_PATH_IMAGE001
at this time, the parameter τ ϵ [ -1, 1], and thus, the signal-to-noise ratio can be expressed as a Chebyshev polynomial:
Figure 812847DEST_PATH_IMAGE042
in the formula,
Figure 391596DEST_PATH_IMAGE003
the signal-to-noise ratio of the fitted direct signal is represented,nin order to improve the fitting accuracy of the signal-to-noise ratio, in this embodiment, a polynomial of order 5 is selected for fitting, that is, n =5, so as to improve the fitting accuracy of the signal-to-noise ratio.C i Chebyshev polynomial coefficients, chebyshev polynomials, being signal-to-noise ratio componentsT i The recurrence formula is:
Figure 904616DEST_PATH_IMAGE043
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:
Figure 240920DEST_PATH_IMAGE044
in the formula:S r_n for signals that contain only the signal-to-noise ratio and noise of the reflected signal,S o representing the signal-to-noise ratio that was originally just received,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, and a signal-to-noise ratio signal only retaining the reflection signal is obtained. Since the high-frequency random noise included in the signal-to-noise ratio obeys normal distribution, the gaussian low-pass filter is used for denoising in the embodiment, 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:
assuming that the signal to be processed containing the signal-to-noise ratio and noise of the reflected signal is represented asS r_n (x) The signal after filtering which only contains the signal-to-noise ratio of the reflected signal is represented asS r The calculation procedure is as follows:
Figure 890732DEST_PATH_IMAGE006
in the formula,G(x) Is a gaussian function, the symbol denotes the convolution operation,xrepresenting the SNR sequence, the high-frequency random noise can be effectively filtered through the processing of the process, the signal only retaining the SNR of the reflected signal is obtained,
Figure 933774DEST_PATH_IMAGE007
in the formula,G(x) Is a Gaussian kernel function, the width of the Gaussian function is
Figure 680013DEST_PATH_IMAGE008
It is decided that,
Figure 300350DEST_PATH_IMAGE008
can be calculated from the standard deviation of the signal-to-noise ratio, as a distribution parameter of the gaussian function,eandπrespectively, are natural constants.
By the technical scheme, high-frequency random noise can be effectively filtered, a signal only retaining the signal-to-noise ratio of the reflection signal is obtained, and the signal is expressed in the subsequent steps.
And S1-1-3, calculating the effective reflection height of the satellite through an adaptive sliding window and a weighted 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 through 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 signal-to-noise ratio signal by adaptive sliding algorithmS 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 size is 3, i.e. a group 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 signal-to-noise ratio signal to be processed is divided intomGroups, which can be represented as:
Figure 423027DEST_PATH_IMAGE045
that is to say havemA set of signals of different signal-to-noise ratio lengths can be used to calculate the satellite effective reflection height.
(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, the power spectrum of the signals is firstly calculated, and the corresponding frequency when the energy is highest is found outf. Then pass throughfAnd solving the corresponding satellite effective reflection height according to the relation between the 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:
Figure 586155DEST_PATH_IMAGE046
in the formula,P x ( f ) Is at a frequency off The power of the periodic signal of (a);S r (x) In order to contain the signal-to-noise ratio of the reflected signal,Nin order to be the length of the sequence,τfor time shift invariants, this can be calculated by:
Figure 565612DEST_PATH_IMAGE047
can find P x (f) F value when reaching peak value, and then calculating the effective reflection height of the corresponding satellite according to the relationship between the frequency and the effective reflection height
Figure 79770DEST_PATH_IMAGE048
The conversion relation is:
Figure 68455DEST_PATH_IMAGE049
in the formula, λ is the wavelength of the signal and is a known quantity, and therefore, the effective reflection height of the satellite can be obtained by the above formula
Figure 86090DEST_PATH_IMAGE050
It should be noted that λ is known quantity and can be directly obtained, and the wavelength of the GPS L1 frequency band is: 19.0cm, and the wavelength of the Beidou B1 frequency band is 19.2 cm. Belonging to the common knowledge in the field. Therefore, the effective reflection height of the satellite can be obtained by the above equation.
Repeating the above calculation processmRespectively resolving the group signal-to-noise ratio signals to obtainmThe effective reflection is high, expressed as:
Figure 767607DEST_PATH_IMAGE051
averaging the calculated satellite effective reflection heights to obtain the satellite effective reflection height of the day
Figure 300219DEST_PATH_IMAGE052
WhereiniRepresents the day of 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
Figure 702382DEST_PATH_IMAGE053
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:
Figure 699156DEST_PATH_IMAGE054
in the formula,w i which represents the weighting coefficient(s) of the,idata representing the number of days, ten days in total, measured on the day of the dayi=10, and so on, the weights decrease gradually.
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.
High effective reflection through satellite
Figure 426941DEST_PATH_IMAGE055
Combining the signal-to-noise ratio signal containing the reflection signal, the phase value of the signal-to-noise ratio can be obtained through the relation between the signal-to-noise ratio signal and the reflection signal, and the specific formula is as follows:
Figure 509167DEST_PATH_IMAGE056
in the formula, A and
Figure 980599DEST_PATH_IMAGE020
the amplitude and phase, which represent the signal-to-noise ratio, are unknowns in the formula.
Furthermore, considering that the signal-to-noise ratio and the altitude angle are in a nonlinear relationship at the moment, in order to improve the precision of parameter calculation, an extended kalman filter algorithm is adopted for the calculation. The extended kalman filter algorithm is called by an existing function module in the Matlab platform, and here we only make a simple algorithm description, as shown below:
first, for a nonlinear system, the state equation and observation equation of its state space can be expressed as:
Figure 97460DEST_PATH_IMAGE057
wherein the first term is a state equation, the second term is an observation equation,x k andz k the actual state vector and the observation vector, respectively. Wherein,x k is formed byAAnd
Figure 527304DEST_PATH_IMAGE020
the unknown state quantity of the unknown quantity composition can be expressed asx k A
Figure 580316DEST_PATH_IMAGE020
);z k Is formed by
Figure 386598DEST_PATH_IMAGE022
The observed state quantity, i.e. the output quantity of the system, of known quantity composition is expressed asz k
Figure 561227DEST_PATH_IMAGE022
);u k As a function of the system, as input quantities to the equation of state, from known quantities
Figure 427552DEST_PATH_IMAGE048
AndEand (4) forming.kWhich is indicative of a sequence of signals,f to relate tok-1 andka non-linear function of the state of the system at the moment,his a state quantityx k And observed quantityz k The non-linear function of interest is,w k andv 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 estimate of the state vector and the observed quantity for which the noise value comes at time k can be ignored, and therefore the above equation can be approximately written as:
Figure 94157DEST_PATH_IMAGE058
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:
Figure 32026DEST_PATH_IMAGE059
therefore, the signal-to-noise ratio phase reference value can be obtained by solving
Figure 795582DEST_PATH_IMAGE020
Wherein, the flow of the Taylor first-order expansion algorithm is as follows:
1) inputting initial conditions, inputting
Figure 504912DEST_PATH_IMAGE060
And
Figure 314605DEST_PATH_IMAGE061
Figure 462690DEST_PATH_IMAGE062
is thatk-1 a posteriori state estimates of the system at time instance,
Figure 18436DEST_PATH_IMAGE063
is thatk-1 time instant systematic error covariance;
2) parameter prediction, the state prediction equation is:
Figure 288880DEST_PATH_IMAGE064
wherein, in the process,
Figure 992394DEST_PATH_IMAGE065
is thatkThe posterior state estimation of the time system, the error covariance prediction equation is:
Figure 881853DEST_PATH_IMAGE066
Figure 947898DEST_PATH_IMAGE067
is thatkState of the moment
Figure 264610DEST_PATH_IMAGE068
The covariance of the prediction is determined by the prediction,Ais thatfFunction pairxThe jacobian matrix after the partial derivation is solved,Wis thatfFor is towThe jacobian matrix after the partial derivation is solved,Qis the process noise covariance matrix;
a correction module, the system gain equation expressed as:
Figure 783316DEST_PATH_IMAGE069
whereinHis thathTo pairxAfter the derivation of the partial derivatives, the jacobian matrix is obtained,Vis thathTo pairvThe jacobian matrix after the partial derivation is solved,Ris an observed noise covariance matrix, where the filter equation is expressed as:
Figure 273203DEST_PATH_IMAGE070
in the formula, the meaning of the parameters is introduced before, and the error covariance update equation is expressed as:
Figure 397017DEST_PATH_IMAGE071
whereinIthe unit matrix 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 filter calculation, and may be obtained through a least square algorithm.
It can be understood that, in the above technical solution, the signal-to-noise ratio phase reference value database is established. By continuously acquiring the navigation data of the Beidou and GPS systems in non-rainy days, and by utilizing the process of the step S1, the sub-system respectively calculates the corresponding signal-to-noise ratio phase values according to the satellite sub-frequency, 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
Figure 946947DEST_PATH_IMAGE072
The phase value of B2 frequency band is
Figure 562736DEST_PATH_IMAGE073
. 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 ratio phase reference value by using the soil humidity measurement model and combining the signal-to-noise ratio 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 the signal-to-noise ratio phase reference value of the Beidou system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 187140DEST_PATH_IMAGE040
in the formula,
Figure 571985DEST_PATH_IMAGE039
the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the measurement day is obtained by processing an original observation value acquired through a Beidou system on the measurement day through step S1, wherein VWC represents the soil humidity required by the measurement day;
it should be noted that the formula is an empirical formula for solving the soil moisture by the 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 solving the mean 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, and B3 is 1268.52MHz, which belongs to the common knowledge in the field.
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 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 the signal-to-noise ratio phase reference value of the GPS system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 620712DEST_PATH_IMAGE040
in the formula,
Figure 786114DEST_PATH_IMAGE039
the VWC represents the soil humidity required, and the signal-to-noise ratio phase value on the measurement day is a value obtained by preprocessing the original observation value acquired by the GPS on the measurement day through step S1.
S3-2, respectively calculating corresponding soil humidity by using 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, L2 is 1227.60MHz, and is common knowledge in the field.
And S4, combining the multi-frequency measurement of the Beidou and GPS dual systems, averaging the results, and outputting the measurement result.
And (4) solving the soil humidity of the final measurement day by using the Beidou and GPS soil humidity obtained by calculation in the step S3 and the step S4, namely, taking the average value of the Beidou and GPS soil humidity, and expressing as follows:
Figure 821066DEST_PATH_IMAGE074
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 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 (9)

1. A Beidou/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 the signal-to-noise ratio of the reflected signal
Connecting and collecting original observation values of the Beidou and the GPS on 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;
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 the soil humidity measurement model and combining the signal-to-noise phase reference values in the 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 multi-frequency measurement of the Beidou and GPS dual systems, averaging the results, and outputting the measurement result.
2. The Beidou/GPS soil humidity measurement method based on sliding algorithm and weighting strategy as claimed in claim 1, wherein the acquisition method of signal-to-noise ratio phase reference value in step S1 is:
s1-1, processing the original data
S1-1-1, filtering 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;
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 an adaptive sliding window algorithm 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.
3. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy according to claim 2, characterized in that the principle of the high-order Chebyshev polynomial fitting algorithm in the step S1-1-1 is expressed as follows:
in the time periodt 0 , t 0t]The signal-to-noise ratio signal is fitted with a Chebyshev polynomial, wheret 0 Is the starting epoch of the signal-to-noise ratio, and ΔtThe length is chosen for the signal-to-noise ratio,
to normalize the signal-to-noise ratio, the time is converted as follows:
Figure 484625DEST_PATH_IMAGE001
at this time, the parameter τ ϵ [ -1, 1], and thus, the signal-to-noise ratio can be expressed as a Chebyshev polynomial:
Figure 237817DEST_PATH_IMAGE002
in the formula,
Figure 243820DEST_PATH_IMAGE003
the signal-to-noise ratio of the fitted direct signal is shown,nfor the order of the chebyshev polynomial,C i chebyshev polynomial coefficients, chebyshev polynomials, being signal-to-noise ratio componentsT i The recurrence formula is:
Figure 740660DEST_PATH_IMAGE004
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:
Figure 784227DEST_PATH_IMAGE005
in the formula:S r_n for signals that contain only the signal-to-noise ratio and noise of the reflected signal,S o representing the signal-to-noise ratio that was originally just received,S d the signal-to-noise ratio of the direct signal of the corresponding point is obtained after Chebyshev polynomial fitting.
4. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy according to claim 3, characterized in that the Gaussian low-pass filtering noise reduction algorithm processing in the step S1-1-2 is expressed as follows:
assuming that the signal to be processed containing the signal-to-noise ratio and noise of the reflected signal is represented asS r_n (x) And the signal after filtering which only contains the signal-to-noise ratio of the reflected signal is expressed asS r The calculation procedure is as follows:
Figure 177162DEST_PATH_IMAGE006
in the formula,G(x) Is a gaussian function, the symbol denotes the convolution operation,xrepresenting the SNR sequence, the high-frequency random noise can be effectively filtered through the processing of the process, the signal only retaining the SNR of the reflected signal is obtained,
Figure 670460DEST_PATH_IMAGE007
in the formula,G(x) Is a Gaussian kernel function, the width of the Gaussian function is
Figure 502150DEST_PATH_IMAGE008
It is decided that,
Figure 272660DEST_PATH_IMAGE008
can be calculated from the standard deviation of the signal-to-noise ratio, as a distribution parameter of the gaussian function,eandπrespectively, are natural constants.
5. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy as claimed in claim 4, wherein in the step S1-1-3, the satellite effective reflection high obtaining method comprises:
(1) determination of signal-to-noise ratio signal by adaptive sliding algorithmS 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 altitudes, the signals are 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, 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 signal to be processed is divided intomGroups, which can be represented as:
Figure 226709DEST_PATH_IMAGE009
that is to say havemThe 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:
Figure 817091DEST_PATH_IMAGE010
in the formula,P x ( f ) Is at a frequency off The power of the periodic signal of (a);S r (x) In order to contain the signal-to-noise ratio of the reflected signal,Nin order to be the length of the sequence,τthe time shift invariant can be calculated by the following formula:
Figure 45947DEST_PATH_IMAGE011
can find outP x ( f ) At the time of peakfValue ofThen the corresponding satellite effective reflection height is obtained through the relation between the frequency and the effective reflection height
Figure 139805DEST_PATH_IMAGE012
Has a conversion relation of
Figure 530335DEST_PATH_IMAGE013
In the formula,λis the wavelength of the signal, and is a known quantity, therefore, the effective reflection height of the satellite can be obtained by the above formulah e
Repeating the above calculation processmRespectively resolving the group signal-to-noise ratio signals to obtainmThe effective reflection is high, expressed as:
Figure 608012DEST_PATH_IMAGE014
averaging the calculated satellite effective reflection heights to obtain the satellite effective reflection height of the day
Figure 374980DEST_PATH_IMAGE015
WhereiniRepresents the day of 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
Figure 854503DEST_PATH_IMAGE016
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:
Figure 88038DEST_PATH_IMAGE017
in the formula,w i which represents the weighting coefficient(s) of the,idata representing the number of days, ten days in total, measured on the day of the dayi=10, and so on, the weights decrease gradually.
6. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy as claimed in claim 5, wherein in the step S1-2, the method for obtaining the SNR phase reference value comprises:
high effective reflection through satellite
Figure 777645DEST_PATH_IMAGE018
Combining signal-to-noise ratio signals containing the reflected signal
Figure 958091DEST_PATH_IMAGE019
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:
Figure 416754DEST_PATH_IMAGE020
in the formula, A and
Figure 758874DEST_PATH_IMAGE021
the amplitude and phase representing the signal-to-noise ratio are unknown quantities in the formula, and E is the altitude angle at the epoch time and is a known quantity, which is obtained by preprocessing in step S1
Figure 607881DEST_PATH_IMAGE022
Signal-to-noise ratio signal obtained in step S1-1-2 and including the reflected signal
Figure 657264DEST_PATH_IMAGE023
Also known quantity, from step S1-1-3, and therefore, can be solved for A and A by the extended Kalman Filter Algorithm
Figure 908117DEST_PATH_IMAGE021
The amount of the unknown quantity is,
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:
Figure 811351DEST_PATH_IMAGE024
wherein the first term is a state equation and the second term is an observation equation,x k Andz k respectively, the actual state vector and the observation vector, wherein,x k is formed byAAnd
Figure 85337DEST_PATH_IMAGE021
the unknown state quantity of the unknown quantity composition can be expressed asx k A
Figure 732219DEST_PATH_IMAGE021
);z k Is formed by
Figure 40841DEST_PATH_IMAGE019
The observed state quantity, i.e. the output quantity of the system, of known quantity composition is expressed asz k
Figure 787080DEST_PATH_IMAGE019
);u k As a system function, as input quantities of a state equation, from a known quantity andEthe components of the composition are as follows,kwhich is indicative of a sequence of signals,
f is related tok-1 andka non-linear function of the state of the system at the moment,his a state quantityx k And observed quantityz k The non-linear function of interest is,w k andv k respectively process noise and observation noise, and is Gaussian white noise with mutually independent mean value of zero,
in practical applications, it cannot be determinedw k Andv k specific values at each moment, but negligible noise values comekThe state vector at that time and the estimate of the observed quantity, and therefore the above equation can be approximated by:
Figure 672997DEST_PATH_IMAGE025
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:
Figure 998936DEST_PATH_IMAGE026
therefore, the signal-to-noise ratio phase reference value can be obtained by solving
Figure 21118DEST_PATH_IMAGE021
7. The Beidou/GPS soil moisture measurement method based on the sliding algorithm and the weighting strategy as claimed in any one of claims 1 to 6, wherein the specific method of the step S2 is as follows:
s2-1, reading the signal-to-noise ratio phase reference value of the Beidou system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 141521DEST_PATH_IMAGE027
in the formula,
Figure 452417DEST_PATH_IMAGE028
the method comprises the steps that the difference between a signal-to-noise ratio phase value and a signal-to-noise ratio phase reference value on the day of measurement is obtained by preprocessing an original observation value obtained through a Beidou system on the day of measurement in step S1, wherein VWC represents the soil humidity required by the day of measurement;
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 Three frequencies of the Beidou satelliteRespectively as follows: b1 is 1575.42MHz, B2 is 1176.45MHz, and B3 is 1268.52MHz, which belongs to the common knowledge in the field.
8. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy as claimed in claim 7, wherein the specific method of the step S3 is as follows:
s3-1, reading the signal-to-noise ratio phase reference value of the GPS system from the signal-to-noise ratio phase reference value database, and performing combined calculation by combining the signal-to-noise ratio phase value of the day to obtain the soil humidity of the measurement day, wherein the calculation formula is as follows:
Figure 441101DEST_PATH_IMAGE029
in the formula,
Figure 458736DEST_PATH_IMAGE028
the difference between the signal-to-noise ratio phase value and the signal-to-noise ratio phase reference value on the measurement day is obtained by preprocessing the original observation value acquired by the GPS on the measurement day through the step S1, wherein the VWC represents the soil humidity required by the measurement day;
s3-2, respectively calculating corresponding soil humidity by using 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 The dual-frequency signal frequencies of the GPS satellite are respectively as follows: l1 is 1575.42MHz, L2 is 1227.60MHz, and is common knowledge in the field.
9. The Beidou/GPS soil moisture measurement method based on sliding algorithm and weighting strategy as claimed in claim 8, wherein the implementation method of the step S4 is:
and (4) solving the soil humidity of the final measurement day by using the Beidou and GPS soil humidity obtained by calculation in the step S3 and the step S4, namely, taking the average value of the Beidou and GPS soil humidity, and expressing as follows:
Figure 874674DEST_PATH_IMAGE030
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.
CN202210770775.1A 2022-07-02 2022-07-02 Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy Active CN114839354B (en)

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 true CN114839354A (en) 2022-08-02
CN114839354B 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 (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101900692A (en) * 2010-06-18 2010-12-01 武汉大学 Method for measuring large-area soil humidity
CN102968552A (en) * 2012-10-26 2013-03-13 郑州威科姆科技股份有限公司 Satellite orbit data estimation and correction method
CN104020180A (en) * 2014-06-19 2014-09-03 武汉大学 Soil humidity inversion method based on low elevation signal received by Beidou base station
CN104224181A (en) * 2014-09-26 2014-12-24 中国科学院生物物理研究所 SAR (Specific Absorption Rate) real-time monitoring system and method of multi-channel magnetic resonance imaging equipment
CN105352979A (en) * 2015-12-08 2016-02-24 武汉大学 Soil humidity estimation method based on Beidou GEO satellite signals
CN105787474A (en) * 2016-03-29 2016-07-20 武汉瑞能电力设备有限责任公司 Processing method of bridge vibration monitoring 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
CN106290408A (en) * 2016-07-21 2017-01-04 清华大学 Based on the soil moisture measurement method running GNSS station signal-to-noise ratio data continuously
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
CN109932735A (en) * 2019-03-25 2019-06-25 中国铁路设计集团有限公司 The localization method of the short baseline single-frequency simple epoch solution of Beidou
CN111399015A (en) * 2020-04-07 2020-07-10 中船重工鹏力(南京)大气海洋信息系统有限公司 Dual-mode satellite fusion positioning method suitable for ship traffic management system
CN111474209A (en) * 2020-03-13 2020-07-31 山东航向电子科技有限公司 Real-time soil humidity measuring method based on intelligent terminal equipment
CN212692782U (en) * 2020-09-10 2021-03-12 安徽阡陌网络科技有限公司 High-precision lightweight handheld land area measuring device
CN112505068A (en) * 2020-11-03 2021-03-16 桂林理工大学 Surface soil humidity multi-satellite combined inversion method based on GNSS-IR
CN112782689A (en) * 2020-12-29 2021-05-11 西南交通大学 Multi-satellite data fusion GNSS-IR soil humidity monitoring method
CN113049777A (en) * 2021-03-12 2021-06-29 北京航空航天大学 Device for measuring soil humidity through GNSS direct reflection signal carrier interference
WO2022007211A1 (en) * 2020-07-10 2022-01-13 自然资源部第一海洋研究所 Gnss-based real-time high-precision wave measurement method and apparatus
CN113959329A (en) * 2021-12-17 2022-01-21 西南交通大学 Snow depth inversion method based on multi-satellite data fusion
CN215953863U (en) * 2020-12-25 2022-03-04 武汉大学 Beidou/GNSS water vapor real-time inversion device
CN114355411A (en) * 2021-12-22 2022-04-15 杭州电子科技大学 Flood detection method based on Beidou or GPS carrier-to-noise ratio observation value
CN114397425A (en) * 2021-12-22 2022-04-26 杭州电子科技大学 GNSS-IR soil humidity inversion method based on generalized continuation approximation

Patent Citations (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101900692A (en) * 2010-06-18 2010-12-01 武汉大学 Method for measuring large-area soil humidity
CN102968552A (en) * 2012-10-26 2013-03-13 郑州威科姆科技股份有限公司 Satellite orbit data estimation and correction method
CN104020180A (en) * 2014-06-19 2014-09-03 武汉大学 Soil humidity inversion method based on low elevation signal received by Beidou base station
CN104224181A (en) * 2014-09-26 2014-12-24 中国科学院生物物理研究所 SAR (Specific Absorption Rate) real-time monitoring system and method of multi-channel magnetic resonance imaging equipment
CN105352979A (en) * 2015-12-08 2016-02-24 武汉大学 Soil humidity estimation method based on Beidou GEO satellite signals
CN105787474A (en) * 2016-03-29 2016-07-20 武汉瑞能电力设备有限责任公司 Processing method of bridge vibration monitoring data
CN106290408A (en) * 2016-07-21 2017-01-04 清华大学 Based on the soil moisture measurement method running GNSS station signal-to-noise ratio data continuously
CN106125106A (en) * 2016-08-10 2016-11-16 清华大学 The method measuring soil moisture based on the ground Big Dipper/GPS dual-mode survey station
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
CN109932735A (en) * 2019-03-25 2019-06-25 中国铁路设计集团有限公司 The localization method of the short baseline single-frequency simple epoch solution of Beidou
CN111474209A (en) * 2020-03-13 2020-07-31 山东航向电子科技有限公司 Real-time soil humidity measuring method based on intelligent terminal equipment
CN111399015A (en) * 2020-04-07 2020-07-10 中船重工鹏力(南京)大气海洋信息系统有限公司 Dual-mode satellite fusion positioning method suitable for ship traffic management system
WO2022007211A1 (en) * 2020-07-10 2022-01-13 自然资源部第一海洋研究所 Gnss-based real-time high-precision wave measurement method and apparatus
CN212692782U (en) * 2020-09-10 2021-03-12 安徽阡陌网络科技有限公司 High-precision lightweight handheld land area measuring device
CN112505068A (en) * 2020-11-03 2021-03-16 桂林理工大学 Surface soil humidity multi-satellite combined inversion method based on GNSS-IR
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
CN113049777A (en) * 2021-03-12 2021-06-29 北京航空航天大学 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
CN114355411A (en) * 2021-12-22 2022-04-15 杭州电子科技大学 Flood detection method based on Beidou or GPS carrier-to-noise ratio observation value
CN114397425A (en) * 2021-12-22 2022-04-26 杭州电子科技大学 GNSS-IR soil humidity inversion method based on generalized continuation approximation

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
MINGKUN SU等: "BeiDou system satellite-induced pseudorange multipath bias mitigation based on different orbital characteristic for static applications", 《IET RADAR SONAR NAVIG.》 *
SIBO ZHANG等: "Use of GNSS SNR data to retrieve soil moisture and vegetation variables over a wheat crop", 《HYDROL. EARTH SYST. SCI. DISCUSS.》 *
张双成等: "BDS/GPS多卫星解译土壤湿度变化研究", 《测绘科学》 *

Also Published As

Publication number Publication date
CN114839354B (en) 2022-11-18

Similar Documents

Publication Publication Date Title
CN114518586B (en) GNSS precise single-point positioning method based on spherical harmonic expansion
Vousdoukas et al. Performance of intertidal topography video monitoring of a meso-tidal reflective beach in South Portugal
CN109001845B (en) rainfall forecasting method
CN110530901B (en) Medium and small scale soil water monitoring system and method integrating cosmic ray neutron method and unmanned aerial vehicle remote sensing
CN111751286B (en) Soil moisture extraction method based on change detection algorithm
CN116519913B (en) GNSS-R data soil moisture monitoring method based on fusion of satellite-borne and foundation platform
CN115980751A (en) Power law model InSAR troposphere delay correction method
CN115047499A (en) Inversion method and system for satellite-borne GNSS-R soil temperature and humidity
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
Yu et al. Generating daily 100 m resolution land surface temperature estimates continentally using an unbiased spatiotemporal fusion approach
CN107273797B (en) Rice sub-pixel identification method based on water body index variation coefficient
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
Wang et al. Soil moisture retrieval from sentinel-1 and sentinel-2 data using ensemble learning over vegetated fields
CN114821360A (en) Method and device for intelligently extracting crop leaf area index and electronic equipment
CN114611699A (en) Soil moisture downscaling method and device, electronic equipment and storage medium
CN115980317B (en) Foundation GNSS-R data soil moisture estimation method based on corrected phase
CN114839354B (en) Beidou and GPS soil humidity measurement method based on sliding algorithm and weighting strategy
CN112733445A (en) Large-area-scale soil moisture inversion method based on evapotranspiration vegetation index spatial characteristics
CN116029162B (en) Flood disaster inundation range monitoring method and system by using satellite-borne GNSS-R data
CN114252875B (en) High-precision meshing method for imaging altitude data
CN113625307A (en) Landslide monitoring system and method based on GNSS
CN111044489B (en) Method for obtaining atmosphere refractive index height distribution profile based on multi-wavelength measurement
CN110543835A (en) Satellite sea surface salinity remote sensing product precision evaluation method based on triple matching theory
Farah Assessment of UNB3M neutral atmosphere model and EGNOS model for near-equatorial-tropospheric delay correction

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