CN109521444B - Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement - Google Patents
Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement Download PDFInfo
- Publication number
- CN109521444B CN109521444B CN201811230539.0A CN201811230539A CN109521444B CN 109521444 B CN109521444 B CN 109521444B CN 201811230539 A CN201811230539 A CN 201811230539A CN 109521444 B CN109521444 B CN 109521444B
- Authority
- CN
- China
- Prior art keywords
- signal
- adaptive
- gps
- observation
- velocity field
- 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
Images
Classifications
-
- 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/14—Receivers specially adapted for specific applications
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a self-adaptive least square fitting estimation algorithm for a GPS horizontal velocity field of the earth movement, which considers unreasonable factors of the relation between the range scale of a research domain and different covariance matrixes, firstly introduces a distance scale factor to construct a Gaussian empirical covariance function, solves the problem of data utilization rate, and effectively avoids the influence of a cross covariance negative value on the construction of an empirical covariance model; on the basis, influence unbalance of observation information and prior information on model parameters is adjusted by further introducing a self-adaptive factor, and coordination consistency of an observed value weight matrix and an empirical covariance matrix variance factor of a signal value is solved; compared with the conventional least square fitting estimation method which only considers a single factor for improvement, the method can more accurately fit and estimate the regional GPS velocity field, thereby more truly reflecting the motion characteristics of the regional crustal structure and providing more valuable basic data and reference for further researching the dynamic characteristics of the regional structure.
Description
Technical Field
The invention relates to a new algorithm for effective fitting estimation aiming at a GPS horizontal velocity field of crustal movement, in particular to a least square fitting estimation algorithm considering the fusion of distance scale factors and self-adaptive adjustment.
Background
The earth crust movement is a phenomenon that earth crust structure changes and earth crust internal substances are displaced due to changes of earth crust internal stress, so that earth surface displacement is generated, and the earth crust movement can be divided into horizontal movement and vertical movement according to the movement direction of the earth crust. The high-precision monitoring of the earth crust movement can obtain an earth crust movement velocity field, and the velocity field can provide valuable reference for deeply researching the migration and deformation characteristics of the earth crust deep structure. With the rapid development of modern space geodetic techniques, especially the rapid modernization construction of space monitoring techniques represented by Global Navigation Satellite systems (GPS, beiDou, navigation Satellite systems, etc.), the modern space geodetic techniques can be used to realize high-precision monitoring of horizontal earth motion with centimeter-level or even millimeter-level precision within a distance range of hundreds of kilometers or even thousands of kilometers. However, due to the influence of factors such as natural environment, construction position and economic cost, large-area and uniformly-covered GPS monitoring site construction is difficult to develop in some active crustal movement areas, so that regional multi-period and high-precision repeated GPS observation site data is difficult to acquire, and actual movement and deformation states of the crustal in a research area cannot be accurately mastered. Therefore, a reasonable fitting estimation model (a function model and a random model) needs to be constructed by combining mathematical knowledge based on the existing GPS monitoring data, and a regional continuous crustal motion deformation field is acquired to further research regional crustal motion and deformation characteristics.
Common function models include quadric surface models, moving surface models, multi-surface function models and the like; common stochastic models are: the weighted mean method, kriging method, etc.
Although the function model can accurately reflect trend changes of a research domain, the randomness of an approximation field is not considered; the stochastic model is mostly modeled by adopting a statistical method, only the randomness of the approaching field is considered, but the overall trend of the approaching field is ignored.
Disclosure of Invention
The invention aims to provide a self-adaptive least square fitting estimation algorithm for a GPS horizontal velocity field of crustal motion, which is used for fitting estimation to reflect regional crustal structure motion characteristics more accurately and provides valuable basic reference for further researching regional structure dynamic characteristics.
In order to realize the task, the invention adopts the following technical scheme:
an adaptive least square fitting estimation algorithm for a GPS horizontal velocity field of the earth movement comprises the following steps:
and step 6, determining a signal covariance matrix between the estimation point and the observation point according to the determined optimal signal covariance matrix and the determined optimal estimated Euler rotation parameter after introducing the distance scale factor and the adaptive factor, thereby estimating the signal value of any point in the research domain.
Further, the solution equation of the prior global motion trend is as follows:
in the above formula, V o A speed matrix, namely an observed value, formed by horizontal movement speed fields of all GPS monitoring stations in a research domain; n is observation error noise, omega is an Euler rotation parameter, and A is an omega coefficient array;
the method for obtaining the first signal value comprises the following steps:
according to saidThe prior overall motion trend is calculated by using an Euler equation to obtain an overall motion velocity field of a research domain, and V is calculated o Deducting the whole motion velocity field of the research domain, and obtaining a crust motion observation signal of the research domain, namely, a first signal value.
Further, the undetermined parameters required for fitting the empirical gaussian function are constructed according to the priori signal values, and the formula is as follows:
wherein C (0) is a pending parameter, q is the number of monitoring sites participating in modeling,
and eliminating a trend item for the observation speed of the GPS monitoring station i to obtain a first-check signal value.
Further, the distance scale factor is calculated as follows:
in the above formula, R is a distance scale factor, k is another undetermined parameter in empirical Gaussian parameters, d max Representing the distance between two GPS monitoring stations which are farthest away in the research domain;
the empirical gaussian function was constructed as follows:
F(d)＝C(0)exp(-k 2 d 2 )
wherein d is the distance between the GPS monitoring stations.
Further, the construction formula of the observation point signal covariance matrix is as follows:
C ij ＝F(d ij )
in the above formula, d ij Is the spherical distance of observation point i and observation point j;
the initial signal weight matrix P can be obtained by inverting the constructed covariance matrix s (0) (ii) a Contract Observation of noise and SignalAre all 1, i.e. the initial a priori unit weight variance of
Then there are:
in the formula, C oo Is a velocity matrix V o Covariance matrix of C nn Is a V o Observation error auto-covariance matrix, C uo A covariance matrix between the point signal and the observed point signal is estimated.
Further, introducing an adaptive factor, and performing adaptive iteration to determine an optimal estimation value of the euler rotation parameter and an optimal estimation value of the posterior signal, wherein the method comprises the following steps:
establishing a self-adaptive fitting estimation function model based on a Helmert variance component, and then establishing a method equation according to a parameter adjustment principle; constructing a self-adaptive factor according to a normal equation, and then adjusting the signal weight matrix by using the self-adaptive factor to obtain a new signal weight matrix;
and (3) performing adjustment calculation again, solving new posterior Euler vector parameters and posterior signal values at the same time, performing iterative calculation, and after lambda times of iteration, when the error in the fitting point residual error reaches the minimum, considering that the observed noise and the signal variance component are coordinated and consistent, thereby finally obtaining the optimal estimation of the Euler rotation parameter posterior signal.
Compared with the prior art, the invention has the following technical characteristics:
1. aiming at the negative definite problem of establishing the distance and the covariance, the distance scale factor is introduced to establish the Gaussian function, so that the data utilization rate is improved, and the influence of a cross covariance negative value on the establishment of an empirical covariance model is effectively avoided.
2. Aiming at the problem of unbalanced influence of observation information and prior information on model parameters, the invention introduces the adaptive factor for adjustment, improves the harmony of the signal covariance matrix established after considering the distance scale factor and the observed value signal covariance matrix compared with the conventional algorithm, but adds the adaptive adjustment process on the basis, has more definite physical significance of the adaptive factor in the process, and can effectively avoid signal distortion caused by calculation of a pure mathematical algorithm.
3. The least square fitting estimation algorithm which considers the fusion of the distance scale factor and the self-adaptive adjustment can more accurately fit and estimate the motion characteristics of the regional crustal structure, and can provide valuable basic data and reference for further researching the dynamic characteristics of the regional structure.
Drawings
FIG. 1 is a flow chart of the method of the present invention.
FIG. 2 is a high-precision GPS horizontal velocity field of the skyscraper plateau crustal sports obtained by the Chinese crustal sports monitoring network and the Chinese terrestrial observation network.
FIG. 3 is a statistical histogram of fitting residuals of a conventional least square fitting estimation method, a single consideration distance scale factor least square fitting estimation method, a single adaptive least square fitting estimation method, and a distance scale factor-consideration adaptive fusion algorithm to the Tibet plateau crustal motion high-precision GPS horizontal velocity field, respectively.
Detailed Description
In the research of fitting and estimating a horizontal earth crust velocity field of a GPS, a least square fitting estimation method which can simultaneously consider trend and randomness is more and more widely applied in recent years. The least square fitting estimation method can simultaneously take the overall trend of the earth crust motion of the research domain and the randomness of signals into consideration, and guarantees the randomness and the analyticity of the preparation of unknown quantities. The method can consider prior signal information, is optimal theoretically and has reliable fitting estimation values, but the core of the algorithm is how to accurately determine a covariance function for constructing a covariance matrix (strong in martial arts, 2009; jiangson, 2010). For the key core problem of the least square fitting estimation method, relevant scholars also propose a corresponding strategy: if the complexity of researching the domain scale range and the earth crust activity is considered, fitting an anisotropic covariance function to ensure the anisotropy of a random signal and the like; and (4) considering the complexity of observation data and the coordination consistency problem of various types of data, and adjusting the prior weight ratio of observation information and signals by using a variance component estimation method. However, the above strategies all involve or are improved for a certain factor, and especially, the comprehensive influence of unreasonable factors of research domain range scale and different covariance matrix relations on the construction of a reasonable empirical covariance function is not considered at the same time, so that the speed field fitting estimation result cannot reach the optimum.
Therefore, the invention considers the unreasonable factor of the relation between the research domain range scale and different covariance matrixes, introduces the distance scale factor to construct more reasonable covariance function and covariance matrix, and further utilizes the variance component estimation to achieve the coordination of the observed value and the signal prior weight ratio on the basis.
When the invention carries out least square fitting estimation on the earth crust movement GPS horizontal velocity field, two key problems exist in the least square fitting estimation:
(1) When the empirical Gaussian function is fitted to construct the covariance matrix, the general method mostly adopts a reduction avoidance measure for negative covariance, so that the fitting accuracy is low due to the reduction of the data utilization rate.
(2) The calculation of model parameters is affected in an unbalanced manner by the observation information and the prior information, namely, variance factors of an observation value weight matrix and an empirical covariance matrix of signal values cannot be coordinated and consistent, so that the calculation result is inconsistent with the actual result.
Aiming at the problems, firstly, a distance scale factor is introduced to construct a Gaussian empirical covariance function, so that the problem of data utilization rate is solved, and the influence of a cross covariance negative value on the construction of an empirical covariance model is effectively avoided; on the basis, influence unbalance of observation information and prior information on model parameters is adjusted by further introducing adaptive factors, and coordination and consistency of an observation value weight matrix and an empirical covariance matrix variance factor of a signal value are ensured. Compared with the conventional least square fitting estimation method which only considers a single factor for improvement, the novel algorithm obviously improves the fitting precision of the GPS horizontal velocity field, can more accurately fit and estimate the regional GPS velocity field, can more accurately fit and estimate the regional crustal structure motion characteristics, and can provide important and valuable basic reference for deeply researching the regional current crustal structure motion state and the internal mechanism characteristics. The specific technical scheme of the invention is as follows:
an adaptive least square fitting estimation algorithm for a GPS horizontal velocity field of the crustal movement comprises the following steps:
calculating the integral motion velocity field of the research domain according to an Euler equation to further obtain a crustal motion observation signal of the research domain, namely a first-check signal value;
wherein, the solution equation of the prior overall motion trend is as follows:
in the above formula, V o A speed matrix, namely an observed value, formed by horizontal movement speed fields of all GPS monitoring stations in a research domain; n is observation error noise, omega is an Euler rotation parameter, A is an omega coefficient array, and can be determined by GPS monitoring site coordinates, so that A omega represents the prior overall motion trend of a research domain (random signal values are not considered here or are reduced to observation noise). In the above formula, q × 1, q × 3, 3 × 1, and q × 1 represent the matrix V, respectively o A, Ω and n dimensions.
The method for obtaining the first signal value comprises the following steps:
according to the prior overall motion trend, an Euler equation is utilized to calculate and obtain an overall motion velocity field of a research domain, and the horizontal motion velocity of all GPS monitoring stations in the research domainVelocity matrix V of field composition o And deducting the integral motion velocity field of the research domain to obtain a crustal motion observation signal of the research domain, namely, a first signal value, and using the first signal value to subsequently construct an empirical Gaussian function.
wherein C (0) is the autocovariance of the GPS monitoring station, namely an undetermined parameter in the empirical Gaussian parameters; q is the number of monitoring stations participating in modeling, and the meaning is the same as q in the step 1;
the calculation expression of the first-check signal value after the trend item is removed from the observation speed of the GPS monitoring station i and j is as follows: v s ＝V o -AΩ，
Are each V s The ith, jth line in (1); q. q of di Is shown at d i Number of GPS monitoring site pairs contained within range segment, C (d) i ) To divide into distance segments d i The cross covariance corresponding to the distance segment.
The following step starts with the fitting of an empirical gaussian function, from which the distance scale factors and the adaptation factors proposed by the present invention are also introduced step-wise.
wherein, the calculation formula of R is as follows:
in the above formula, k is another undetermined parameter in the empirical Gaussian parameter, d max Indicating the distance between the two GPS monitoring sites that are furthest apart within the study. D is taken as the distance parameter d between the GPS monitoring stations max Determining corresponding k by taking different values for R under the criterion that the time R approaches infinity, and further determining an empirical Gaussian function according to an empirical Gaussian function expression taking distance scale factors into consideration, wherein the expression is as follows:
F(d)＝C(0)exp(-k 2 d 2 )
wherein d is the spherical distance between the GPS monitoring stations. And (3) analyzing errors in the fitting result residual error, selecting a value k with a smaller error in the residual error, considering the range size of a research domain, properly adjusting R, controlling the Gaussian function through k to be more matched with the research domain scale so as to obtain a better fitting result, and determining a specific mathematical expression of the empirical Gaussian function after introducing a distance scale factor.
When the Gaussian function is fitted in this step by conventional least squares fitting, it is usually calculated from C (0), C (d) in step 2 i ) D, mixing C (d) i ) The negative values in the data are removed to ensure that the data used for fitting are consistent with the positive character of a Gaussian function, then each distance section is taken as a horizontal axis, the cross covariance between sites of each distance section is taken as a vertical axis (the cross covariance between sites corresponding to the 0 distance section is the self covariance of the sites), and the empirical Gaussian function F (d) = aexp (-k) according to the distance scale of an unconscious research domain is taken into consideration 2 d 2 ) (the empirical Gaussian function coefficient a estimated by the conventional least square fitting is obtained by the least square, and the coefficient C (0) of the empirical Gaussian function considering the distance scale factor is directly calculated in the step 2), the Gaussian functions in the east and west directions are respectively fitted by the least squares to respectively obtain two undetermined parameters a and k in the east direction and two undetermined parameters a and k in the west direction of the empirical Gaussian function estimated by the conventional least square fitting, so far, in the conventional least square fitting estimation, the specific mathematical expression form of the empirical Gaussian function is considered to be completely determined and is used for constructing each subsequent covariance matrix.
However, in the above conventional least square fitting estimation algorithm, the reduction and avoidance measures taken for negative values will cause the data utilization rate to be greatly reduced, and the reliability of the fitting model will be reduced accordingly. In order to improve the utilization rate of data, the invention introduces a distance scale factor R considering the research domain range. After the distance scale factor R is introduced, the research domain scale range is considered, the utilization rate of data is obviously improved, the influence of a cross covariance negative value on the construction of an empirical covariance model is avoided, and the empirical Gaussian function is more accurate and reliable.
Next, an empirical covariance matrix is constructed through the empirical gaussian function constructed in step 3 for subsequent calculation.
the covariance matrix is constructed by the following formula:
C ij ＝F(d ij )
in the above formula, d ij Is the spherical distance of observation point i and observation point j;
the initial signal weight matrix P can be obtained by inverting the constructed covariance matrix s (0) (ii) a The initial a priori unity weight variance of the observation noise and signal are both agreed to be 1, i.e.
Then there are:
in the above formula, the first and second carbon atoms are,C oo is a velocity matrix V o The covariance matrix of (4) is used for describing the correlation and the dispersion of the spatial distribution among the GPS stations; c nn Is a V o The observation error auto-covariance matrix is used for describing the discreteness of the GPS station speed and reflecting the observation precision; c uo To estimate the covariance matrix between the point signals and the observed point signals. The above three covariance matrices, C nn Can be obtained according to the result of the GPS velocity field solution; c oo 、C uo And determining according to an empirical Gaussian function constructed after introducing the distance scale factor.
The posterior Euler rotation parameter can be obtained by solving through the adjustment calculation of the above three formulas
And the posterior signal value of the earth crust motion observation point of the research domain estimated
Solved at this time
Before adaptive iteration, in order to facilitate the parameter difference from the subsequent calculation, the value added with 0 will be
Is written as
Will be provided with
Writing into
So far, the posterior euler rotation parameters and the posterior signal values of the observation points have been solved, and the conventional least square fitting estimation method and the least square fitting estimation method considering only the distance scale factors have ended in this step. However, in the invention, only after the distance scale factor is introduced by considering the research domain scale, the constructed covariance matrix can not be directly optimized, and an unreasonable covariance function can cause the incoordination of the observed value noise and the signal variance factor to generate a system error, so that the self-adaptive adjustment strategy is further integrated and introduced on the basis of the conventional least square fitting estimation for estimating the distance scale factor, the observed value noise and the signal variance factor are coordinated and consistent, and the fitting estimation precision of the GPS horizontal velocity field is further improved; adaptive adjustment is performed as follows.
establishing a function model min based on self-adaptive fitting estimation of Helmert variance component:
min＝V T P e V+αS T P s S
in the above formula, α is an adaptive factor, and V is an observed value V 0 S is the signal value, P e For observation of vector weight arrays, P S Is a signal weight matrix. The model takes signals as virtual observation values, forms two types of observation quantities together with the observation value signals, and can establish the following matrix:
further, according to the parameter adjustment principle, a method equation is established:
several of the aboveIn the formula, A has the same meaning as step 1 and is a coefficient array of Euler rigid motion model parameter vectors omega, B is a coefficient of a signal value, and the shape of the coefficient is [ I0 ]] T Suppose that the estimate point has n 3 I in B is n 1 ×n 1 Order being a unit matrix, 0 being n 3 ×n 1 A zero matrix of order; l has the same meaning as V in step 1 o But is represented by L because it is in the form of an error equation here; v s Is a modified number of signal values, where V s I.e., equal to S;
is the observed signal unit weight variance; sigma e 2 For observation noise unit weight variance; for S and N matrices, the index i is one of the values 1,2, N 1 ＝tr(P s *P s -1 )，n 2 ＝tr(P e *P e -1 ) J is the same as i, i.e. S ii S ij Corresponding to S in the above equation 11 S 12 S 21 S 22 According to the above formulas, σ can be solved e 2 And
constructing an adaptive factor alpha, and then utilizing the adaptive factor alpha to carry out signal weight matrix P S Performing rising or falling adjustment to obtain a new signal weight array P s (1) ：
P s (1) ＝αP s
According to the formula
For newly obtained signal value weight array P s (1) Inverse operation is carried out to obtain a new signal covariance matrix
Then, the adjustment calculation is carried out again, and new posterior Euler vector parameters are solved
And a posteriori signal value
Solving the signal value by using a least square fitting estimation solution formula in the step 4
And the post-test global trend term
And step 5, the primary power difference component estimation re-determines the weight matrix
And σ e 2 As a first iteration calculation, after lambda iterations, when the error in the fitting point residual reaches the minimum, the observation noise and the signal variance component are considered to be coordinated and consistent, and the optimal estimated value of the Euler rotation parameter is finally obtained
Optimum estimation of a posteriori signals
After the adaptive factors are introduced for adjustment, the harmony between the signal covariance matrix established after considering the distance scale factors and the observed value signal covariance matrix is further improved compared with the conventional algorithm, the physical significance for researching the area scale range is realized in the adaptive adjustment process, and the problem of signal distortion caused by only pure mathematical solution can be effectively avoided. The distance scale factor and the adaptive factor are fused with a new algorithm, so that the regional crustal structure motion characteristics can be more accurately fitted and estimated.
in the formula (I), the compound is shown in the specification,
is composed of
Inverse matrix of (C) oo Is a velocity matrix V o The covariance matrix is used for describing correlation and dispersion of spatial distribution among the GPS stations; c nn Is a V o The observation error auto-covariance matrix is used for describing the discreteness of the GPS station speed and reflecting the observation precision; c uo A covariance matrix between the point signal and the observed point signal is estimated. The above three covariance matrices, C nn Can be obtained according to the result of the GPS velocity field solution; c oo 、C uo And determining according to an empirical Gaussian function constructed after introducing the distance scale factor.
Representing the signal value of the point to be estimated, where
The meaning is the same as in step 4
V o And A is the same as V in step 1 o 、A，
Then represents the best estimate of the euler rotation parameter
Represents the best estimate of the euler trend term for the entire study area.
Using the estimated signal value at any point of the domain
And the research on the movement characteristics of the regional crust structure can be carried out by combining the existing geological survey data of the research domain.
Example (c):
in order to verify the optimality and effectiveness of the algorithm provided by the invention, the Qinghai-Tibet plateau area of the strong activity area of the typical crustal structure of China is selected, and the least square fitting estimation new algorithm provided by the invention is tested by utilizing the long-term high-precision GPS horizontal velocity field of the crustal motion of the area. The GPS observation data comes from China crustal motion observation network and China land state observation network, the time span is 672 station data in 1999-2013, 127 continuous operation reference stations and 545 mobile observation stations are contained in the GPS observation data, and the GPS station speed field information is shown in FIG. 2. And performing corresponding fitting estimation calculation on the GPS horizontal velocity field of the Tibet plateau hull motion by respectively utilizing four algorithms of a routine algorithm, a self-adaption algorithm, a distance scale factor considered algorithm, a least square fitting estimation algorithm fusing the distance scale factor considered algorithm and the self-adaption algorithm. The method for analyzing and processing the GPS horizontal velocity field of the Tibet plateau hull motion by the least square fitting estimation method integrating the distance scale factor and the adaptive adjustment comprises the following steps:
(1) Based on the GPS horizontal velocity field of the Tibet plateau hull motion, firstly, the Euler rotation parameters before the experiment are solved by applying a least square method.
(2) The prior Euler rotation parameter of the Qinghai-Tibet plateau obtained in the last step is calculated according to a formula V s ＝V o And calculating to obtain a pre-test signal value by using-A omega, and calculating to obtain a parameter C (0) in the empirical Gaussian function of the research domain by using a formula of calculating C (0) by using the pre-test signal value in the specific implementation step 2.
(3) Introducing a distance scale factor R, and selecting the distance d of two stations with the farthest distance in the research domain max At the determination of d max And selecting different R in a later experiment to determine a parameter k in the empirical Gaussian function of the research area, and further determining a specific mathematical expression of the empirical Gaussian function for subsequently constructing an empirical covariance matrix.
(4) From C ij ＝F(d ij ) Construction of an observation velocity vector V between GPS observation sites in a research domain o Covariance of C oo And estimating a covariance matrix C between the point signal and the observed point signal uo And further, calculating a posterior integral trend term and a posterior signal value according to the conventional least square fitting, estimating and resolving formula mentioned in the step 4 to obtain the integral trend term and the displacement signal value of the station under consideration of the distance scale factor in the research domain.
(5) And introducing a self-adaptive factor to perform self-adaptive iterative calculation of observed value noise of a research domain and a signal after the experiment is solved. Solving sigma according to Helmert variance component estimation formula in step 5 in each iteration e 2 And
constructing an adaptive factor, and then adjusting the signal weight array by using the adaptive factor to obtain a new signal weight array P s (1) Then, adjustment is performed again, and new posterior Euler rotation parameters are solved
And a posteriori signal value
After k iterations, when the error in the fit point residual error is minimum, the observation noise and the signal variance component are considered to be consistent, and the optimal estimated value of the Euler rotation parameter of the research domain is obtained
Sum signal optimum estimation
(6) C of the above construction uo Substituting into formula
In this way, the signal value of any estimated point in the research domain can be calculatedAnd then, the earth crust movement GPS horizontal velocity field in any area of the research domain can be obtained by combining the block body integral movement trend item solved by posterior Euler rotation parameters.
The four algorithms are respectively utilized to carry out corresponding fitting estimation calculation on the GPS horizontal velocity field of the Tibet plateau shell movement, and errors in fitting residual errors of the four algorithms are obtained, as shown in table 1. Table 1 shows that the fitting accuracy of the improved algorithm is significantly higher than the other three classes of algorithms. Fig. 3 fits the residual statistical histogram, also visually showing that the improved algorithm proposed herein has a better fitting effect.
In order to further check the effectiveness of the algorithm, 20 GPS monitoring stations which are high in point location accuracy, good in reliability and stability and distributed uniformly are selected in a research domain to serve as check points, the positions of the check points are shown in FIG. 2, the rest 652 GPS monitoring stations serve as fitting points, and the result of the accuracy of the fitting residual errors of the check points is shown in Table 1. The results shown in table 1 also indicate that the improved algorithm proposed herein is the optimal algorithm among the four schemes, further verifying the effectiveness of the improved algorithm proposed herein.
TABLE 1 check points fitting residual accuracy results
Claims (4)
1. An adaptive least square fitting estimation algorithm for a GPS horizontal velocity field of the crustal movement is characterized by comprising the following steps of:
step 1, solving the prior overall motion trend of a research domain according to the least square principle based on the horizontal motion velocity fields of all GPS monitoring stations in the research domain obtained by observation; calculating the overall motion velocity field of the research domain to further obtain a crust motion observation signal of the research domain, namely a first signal value;
step 2, according to the first signal value, establishing undetermined parameters required by fitting an empirical Gaussian function, wherein the formula is as follows:
wherein C (0) is an undetermined parameter in empirical Gaussian parameters,qthe number of monitored sites participating in the modeling,
for GPS monitoring stationiThe observation speed of (2) eliminates the first signal value after the trend item;
step 3, introducing a distance scale factor considering the research domain range, determining another undetermined parameter in the empirical Gaussian function by selecting different distance scale factors, and then constructing the empirical Gaussian function;
the calculation formula of the distance scale factor is as follows:
in the above-mentioned formula, the compound has the following structure,Ris a function of the distance scale factor,kfor another undetermined parameter in the empirical gaussian parameter,d max the distance between two GPS monitoring stations which are farthest away in the research domain is represented;
the empirical gaussian function was constructed as follows:
whereindMonitoring the spherical distance between the stations by the GPS;
step 4, constructing an observation point signal covariance matrix according to an empirical Gaussian function, and solving posterior Euler rotation parameters and posterior signal values through adjustment calculation;
step 5, introducing a self-adaptive factor, and performing self-adaptive iteration to determine the optimal estimated value of the Euler rotation parameter and the optimal estimated value of the posterior signal;
and step 6, determining a signal covariance matrix between the estimation point and the observation point according to the determined optimal signal covariance matrix and the determined optimal estimated Euler rotation parameter after introducing the distance scale factor and the adaptive factor, thereby estimating the signal value of any point in the research domain.
2. The earth-movement GPS horizontal velocity field adaptive least squares fitting estimation algorithm of claim 1, characterized in that the solution equation of the prior global movement trend is:
in the above formula, the first and second carbon atoms are,V o a speed matrix, namely an observed value, formed by horizontal movement speed fields of all GPS monitoring stations in a research domain;nin order to observe the error noise, it is,
in order to be able to determine the euler rotation parameters,Ais composed of
A coefficient array;
the method for obtaining the first signal value comprises the following steps:
according to the prior overall motion trend, the Euler equation is utilized to calculate and obtain the overall motion velocity field of the research domainV o And deducting the integral motion velocity field of the research domain to obtain a crustal motion observation signal of the research domain, namely, a signal value of the first test.
3. The earth-motion GPS horizontal velocity field adaptive least squares fitting estimation algorithm of claim 2, wherein the construction formula of the observation point signal covariance matrix is as follows:
C ij = F(d ij )
in the above formula, the first and second carbon atoms are,d ij is an observation pointiAnd observation pointjThe spherical distance of (a);
the initial signal weight matrix can be obtained by inverting the constructed covariance matrix
(ii) a The initial a priori unity weight variance of the observation noise and signal are both agreed to be 1, i.e.
、
Then, there are:
in the formula (I), the compound is shown in the specification,C oo is a velocity matrixV o The covariance matrix of (a) is determined,C nn is composed ofV o An error auto-covariance matrix is observed,C uo to estimate the covariance matrix between the point signals and the observed point signals.
4. The earth motion GPS horizontal velocity field adaptive least squares fitting estimation algorithm of claim 1, wherein introducing adaptive factors, performing adaptive iteration to determine the optimal estimate of euler rotation parameters and the optimal estimate of a posteriori signals, comprises:
establishing a self-adaptive fitting estimation function model based on Helmert variance components, and then establishing a method equation according to a parameter adjustment principle; constructing a self-adaptive factor according to a normal equation, and then adjusting a signal weight matrix by using the self-adaptive factor to obtain a new signal weight matrix;
and (3) performing adjustment calculation again, solving new posterior Euler vector parameters and posterior signal values at the same time, performing iterative calculation, and after lambda iterations, when the error in the fitting point residual reaches the minimum, considering that the observed noise and the signal variance component are coordinated and consistent, thereby finally obtaining the optimal estimated value of the Euler rotation parameters and the optimal estimated value of the posterior signals.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811230539.0A CN109521444B (en) | 2018-10-22 | 2018-10-22 | Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811230539.0A CN109521444B (en) | 2018-10-22 | 2018-10-22 | Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109521444A CN109521444A (en) | 2019-03-26 |
| CN109521444B true CN109521444B (en) | 2023-03-14 |
Family
ID=65772865
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811230539.0A Active CN109521444B (en) | 2018-10-22 | 2018-10-22 | Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109521444B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110298013B (en) * | 2019-06-28 | 2023-08-04 | 安徽建筑大学 | Multi-period observation data processing method based on nonlinear Gauss-Helmert model |
| CN113254858B (en) * | 2021-04-26 | 2022-07-29 | 山东科技大学 | Resolving method for ill-conditioned separable nonlinear least square problem and application thereof |
| CN113268869B (en) * | 2021-05-19 | 2022-02-01 | 南方科技大学 | Method and system for monitoring change of earth surface quality |
| CN116955976B (en) * | 2023-09-20 | 2023-11-28 | 航天宏图信息技术股份有限公司 | Earthquake ground surface deformation analysis method and device based on deep learning and Beidou positioning |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6057800A (en) * | 1996-06-28 | 2000-05-02 | State University Of New York | RDOP surface for GPS relative positioning |
| CN102508279A (en) * | 2011-11-18 | 2012-06-20 | 中国测绘科学研究院 | Method for processing GNSS (global navigation satellite system) positioning posture measuring value of satellite navigation system and GNSS positioning posture measuring instrument |
| CN106066901A (en) * | 2016-04-22 | 2016-11-02 | 中南大学 | A kind of datum mark method for analyzing stability of GNSS automatization deformation monitoring |
| CN106441174A (en) * | 2016-09-09 | 2017-02-22 | 桂林电子科技大学 | High slope deformation monitoring method and system |
| CN107102342A (en) * | 2017-04-28 | 2017-08-29 | 武汉大学 | Gps coordinate time series discontinuity based on common-mode error supplies method |
| CN108345009A (en) * | 2018-02-08 | 2018-07-31 | 中国石油大学(华东) | The GPS three-dimensional steam chromatography methods and device of no prior information constraint |
-
2018
- 2018-10-22 CN CN201811230539.0A patent/CN109521444B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6057800A (en) * | 1996-06-28 | 2000-05-02 | State University Of New York | RDOP surface for GPS relative positioning |
| CN102508279A (en) * | 2011-11-18 | 2012-06-20 | 中国测绘科学研究院 | Method for processing GNSS (global navigation satellite system) positioning posture measuring value of satellite navigation system and GNSS positioning posture measuring instrument |
| CN106066901A (en) * | 2016-04-22 | 2016-11-02 | 中南大学 | A kind of datum mark method for analyzing stability of GNSS automatization deformation monitoring |
| CN106441174A (en) * | 2016-09-09 | 2017-02-22 | 桂林电子科技大学 | High slope deformation monitoring method and system |
| CN107102342A (en) * | 2017-04-28 | 2017-08-29 | 武汉大学 | Gps coordinate time series discontinuity based on common-mode error supplies method |
| CN108345009A (en) * | 2018-02-08 | 2018-07-31 | 中国石油大学(华东) | The GPS three-dimensional steam chromatography methods and device of no prior information constraint |
Non-Patent Citations (4)
| Title |
|---|
| Crustal strain fields in the surrounding areas of the Ordos Block,central China, estimated by the least-squares collocation technique;Wei Qu et al.;《Journal of Geodynamics》;20170127;第1-11页 * |
| 基于欧拉矢量的中国大陆地壳水平运动自适应拟合推估模型;杨元喜 等;《中国科学:地球科学》;20111231;第41卷(第08期);第1116-1125页 * |
| 最小二乘配置与REHSM求解GPS应变场的方法;管守奎 等;《大地测量与地球动力学》;20150831;第35卷(第04期);第604-607页 * |
| 自适应拟合推估在地壳形变分析中的应用;赵丽华 等;《大地测量与地球动力学》;20090228;第29卷(第01期);第132-134、139页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109521444A (en) | 2019-03-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109521444B (en) | Self-adaptive least square fitting estimation algorithm for GPS horizontal velocity field of crustal movement | |
| CN108981559B (en) | Real-time deformation monitoring method and system based on Beidou foundation enhancement system | |
| CN110058236B (en) | InSAR and GNSS weighting method oriented to three-dimensional surface deformation estimation | |
| CN105136054B (en) | The fine deformation monitoring method of structures and system based on Three Dimensional Ground laser scanning | |
| CN106912105B (en) | Three-dimensional positioning method based on PSO _ BP neural network | |
| CN105301601B (en) | A kind of GNSS ionosphere delay three-dimensional modeling methods suitable for Global Regional | |
| CN107247259B (en) | K distribution sea clutter shape parameter estimation method based on neural network | |
| CN108802770B (en) | High-precision dynamic positioning verification reference for INS enhanced GNSS | |
| CN108668358B (en) | Arrival time-based cooperative positioning method applied to wireless sensor network | |
| CN104729486A (en) | Bathymetric surveying method without tide observation based on quasigeoid refinement | |
| CN109283562A (en) | Three-dimensional vehicle localization method and device in a kind of car networking | |
| CN102841385A (en) | Local geomagnetic chart constructing method based on multi-fractal Krigin method | |
| US9213100B1 (en) | Bearing-only tracking for horizontal linear arrays with rapid, accurate initiation and a robust track accuracy threshold | |
| CN105759311A (en) | Near-real time earthquake source position positioning method | |
| CN109581281A (en) | Moving objects location method based on reaching time-difference and arrival rate difference | |
| CN109490855A (en) | A kind of trailer-mounted radar scaling method, device and vehicle | |
| Yang et al. | Horizontal crustal movement in China fitted by adaptive collocation with embedded Euler vector | |
| CN106644868B (en) | A kind of measuring method of two dimension not convex Random Aggregate ambient interfaces concentration | |
| CN106556877B (en) | A kind of earth magnetism Tonghua method and device | |
| CN110018501A (en) | A kind of multimode accurate one-point positioning method adjusted based on stochastic model On-line Estimation between system | |
| CN114004104A (en) | CORS site selection method based on checkerboard test | |
| CN107391794B (en) | Typhoon continuous three-dimensional wind field inversion method | |
| CN108447126B (en) | Laser point cloud precision evaluation method of mobile measurement system based on reference plane | |
| CN102043156A (en) | Adjustment processing method for measuring two-dimensional baseline vector network by GPS (Global Position System) | |
| CN111505575B (en) | Sensor selection method aiming at TDOA (time difference of arrival) location based on conversion TOA (time of arrival) model |
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 |