CN108332649B - Landslide deformation comprehensive early warning method and system - Google Patents

Landslide deformation comprehensive early warning method and system Download PDF

Info

Publication number
CN108332649B
CN108332649B CN201810123208.0A CN201810123208A CN108332649B CN 108332649 B CN108332649 B CN 108332649B CN 201810123208 A CN201810123208 A CN 201810123208A CN 108332649 B CN108332649 B CN 108332649B
Authority
CN
China
Prior art keywords
deformation
monitoring
displacement
early warning
landslide
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
Application number
CN201810123208.0A
Other languages
Chinese (zh)
Other versions
CN108332649A (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.)
Guilin University of Electronic Technology
Original Assignee
Guilin University of Electronic Technology
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 Guilin University of Electronic Technology filed Critical Guilin University of Electronic Technology
Priority to CN201810123208.0A priority Critical patent/CN108332649B/en
Publication of CN108332649A publication Critical patent/CN108332649A/en
Application granted granted Critical
Publication of CN108332649B publication Critical patent/CN108332649B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/004—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring coordinates of points
    • 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
    • G—PHYSICS
    • G08—SIGNALLING
    • G08B—SIGNALLING OR CALLING SYSTEMS; ORDER TELEGRAPHS; ALARM SYSTEMS
    • G08B21/00—Alarms responsive to a single specified undesired or abnormal condition and not otherwise provided for
    • G08B21/02—Alarms for ensuring the safety of persons
    • G08B21/10—Alarms for ensuring the safety of persons responsive to calamitous events, e.g. tornados or earthquakes
    • G—PHYSICS
    • G08—SIGNALLING
    • G08B—SIGNALLING OR CALLING SYSTEMS; ORDER TELEGRAPHS; ALARM SYSTEMS
    • G08B29/00—Checking or monitoring of signalling or alarm systems; Prevention or correction of operating errors, e.g. preventing unauthorised operation
    • G08B29/18—Prevention or correction of operating errors
    • G08B29/185—Signal analysis techniques for reducing or preventing false alarms or for enhancing the reliability of the system

Abstract

The invention discloses a landslide early warning method and a system, comprising a landslide deformation monitoring network and a comprehensive early warning unit; the landslide deformation monitoring network comprises a satellite positioning reference station and more than 3 satellite positioning monitoring stations; more than 3 satellite positioning monitoring stations are distributed at different monitoring points on the earth surface of the monitoring target; each satellite positioning monitoring station is respectively connected with a satellite positioning reference station to realize RTK positioning resolving; and the positioning results of all the satellite positioning monitoring stations are transmitted to the comprehensive early warning unit. According to the method, the prediction of deformation displacement is realized by accurately calculating the overall attitude and displacement of the monitored target, and comprehensive early warning is realized by fusing multi-source monitoring information. The method has the characteristics of high early warning accuracy, simple algorithm, easy realization and stronger practicability, and can meet the actual requirement of the mountain landslide deformation early warning.

Description

Landslide deformation comprehensive early warning method and system
Technical Field
The invention relates to the technical field of geological early warning, in particular to a landslide deformation comprehensive early warning method and system.
Background
Due to the influence of factors such as the construction of key projects such as water and electricity, environmental climate change and the like, geological disasters such as landslide, ground deformation, earthquake, land degradation and the like in China are increasingly serious. Landslide, collapse, ground subsidence and the like become one of main geological disasters in areas such as southwest mountainous areas of China, and one of important measures for actively and effectively preventing and treating the landslide and other geological disasters is to observe deformation hidden trouble points such as the landslide for a long time. Once the landslide hidden danger point is deformed, if the landslide hidden danger point cannot be monitored and effectively early warned in time, huge economic loss can be brought to local residents, and even the life safety of the residents is damaged.
Disclosure of Invention
The invention provides a landslide deformation comprehensive early warning method and a landslide deformation comprehensive early warning system, which have the characteristics of simplicity in implementation and high early warning accuracy.
In order to solve the problems, the invention is realized by the following technical scheme:
a landslide deformation comprehensive early warning method specifically comprises the following steps:
step 1, solving three-dimensional position coordinates of monitoring points of monitoring targets of landslide deformation monitoring networks by each satellite positioning monitoring station, and solving displacement according to initial positions of the satellite positioning monitoring stations;
step 2, smoothing the displacement obtained by each satellite positioning monitoring station by using a Kalman filter to obtain an actual measurement value of the deformation displacement; according to the sampling time interval and the displacement difference of the measured values of the deformation displacement at the front and rear moments, the measured values of the deformation speed and the deformation acceleration of the monitoring target are solved;
step 3, extracting the measured values of the deformation displacement of each monitoring point obtained in the step 2, and solving the attitude of the monitoring target by using a least square method, namely the measured values of the course angle, the pitch angle and the roll angle of the monitoring target;
step 4, extracting the measured values of the deformation displacement, the deformation speed and the deformation acceleration of each monitoring point obtained in the step 2, and solving the measured value of the comprehensive displacement of the monitoring target by using an extended Kalman filter;
step 5, setting an initial epoch and a length of a prediction sliding window, sampling the comprehensive displacement of the monitored target at equal time intervals, utilizing an iterative gray model to realize a predicted value of the deformation displacement of the monitored target, and calculating predicted values of the deformation speed and the deformation acceleration according to the sampling time interval and the displacement difference of the predicted values of the deformation displacement at the previous moment and the next moment;
step 6, carrying out actual measurement landslide grade early warning on the current state of the monitoring target according to a preset actual measurement early warning rule by using the actual measurement values of the deformation displacement, the deformation speed and the deformation acceleration calculated in the steps 3 and 4; and meanwhile, carrying out prediction landslide grade early warning on the subsequent state of the monitored target according to a preset prediction early warning rule by using the predicted values of the deformation displacement, the deformation speed and the deformation acceleration predicted in the step 5.
The specific process of the step 5 is as follows:
step 5.1, setting a prediction sliding window initial epoch, and sampling a measured value of the comprehensive displacement of the monitoring target to obtain an observation sequence of the current prediction sliding window;
step 5.2, accumulating the observation sequence for one time to obtain a data generation sequence, and calculating an adjacent mean value generation sequence according to the data generation sequence;
step 5.3, constructing a Jacobian matrix by utilizing the adjacent mean value generation sequence, and constructing a least square estimation parameter coefficient, namely a development coefficient and a gray action amount according to the Jacobian matrix;
step 5.4, carrying out one-time subtraction on the observation number series by using the development coefficient and the gray effect quantity to obtain a prediction number series of the current prediction sliding window;
step 5.5, continuously pushing the initial observation data forward by using the prediction sliding window, and performing iterative computation by repeating the step 5.2-5.4 to obtain the prediction sequence of each prediction sliding window;
step 5.6, performing weighted accumulation on the prediction series of the prediction sliding windows to obtain a prediction value of the deformation displacement of the monitoring target;
and 5.7, calculating the predicted values of the deformation speed and the deformation acceleration according to the sampling time interval and the displacement difference of the predicted values of the deformation displacement at the previous moment and the next moment.
The specific process of the step 3 is as follows:
step 3.1, extracting the measured values of the deformation displacement of each monitoring point obtained in the step 2;
3.2, establishing a baseline vector by using the three-dimensional position coordinates of any more than 3 non-collinear satellite positioning monitoring stations in the landslide deformation monitoring network, and randomly selecting 2 baseline vectors to calculate the initial values of a pitch angle, a course angle and a roll angle;
3.3, calculating corresponding weight of a base line according to the constructed attitude coefficient matrix, and estimating an attitude correction number by using a least square method;
and 3.4, calculating the attitude of the monitoring target according to the attitude correction number and the initial values of the pitch angle, the course angle and the roll angle, namely the actual measured values of the course angle, the pitch angle and the roll angle of the monitoring target.
The specific process of the step 4 is as follows:
step 4.1, extracting measured values of deformation displacement, deformation speed and deformation acceleration of each monitoring point obtained in the step 2;
step 4.1, constructing a state quantity matrix, a global observation matrix and a measurement relation matrix of the extended Kalman filter by using the deformation displacement and the deformation speed, constructing a state transition matrix of the extended Kalman filter by using a sampling time interval, and resolving a process noise covariance matrix and a process noise covariance matrix of the extended Kalman filter according to the deformation acceleration;
and 4.1, updating time by using the extended Kalman filter, predicting the current system state, correcting the priori estimation value of the predicted state through measurement updating, and obtaining the global optimal estimation of the overall displacement of the monitoring target, namely the measured value of the comprehensive displacement of the monitoring target.
As an improvement, the landslide deformation comprehensive early warning method further comprises the following steps: and 7, collecting rainfall by a rainfall gauge in the landslide deformation monitoring network, collecting the rainfall and collecting the crack displacement by a crack sensor, performing weighted operation on the rainfall and the crack displacement and the measured values and the predicted values of the deformation displacement, the deformation speed and the deformation acceleration obtained in the steps 4 and 5, and performing comprehensive landslide grade early warning on the monitored target according to a preset comprehensive early warning rule.
A landslide deformation comprehensive early warning system for realizing the method comprises a landslide deformation monitoring network and a comprehensive early warning unit; the landslide deformation monitoring network comprises a satellite positioning reference station and more than 3 satellite positioning monitoring stations; more than 3 satellite positioning monitoring stations are distributed at different monitoring points on the earth surface of the monitoring target; each satellite positioning monitoring station is respectively connected with a satellite positioning reference station to realize RTK positioning resolving; and the positioning results of all the satellite positioning monitoring stations are transmitted to the comprehensive early warning unit.
In the above scheme, the satellite positioning monitoring station is a GPS satellite positioning receiver, a BDS satellite positioning receiver and/or a GNSS satellite positioning receiver.
As an improvement, the landslide deformation monitoring network of the landslide deformation comprehensive early warning system further comprises a rain gauge and more than 2 crack sensors; the rain gauges are distributed on the ground surface of the monitoring target; more than 2 crack sensors are distributed at different monitoring points on the ground surface of the monitoring target; and the monitoring data of the rain gauge and all the crack sensors are transmitted to the comprehensive early warning unit.
In the scheme, the number of the satellite positioning monitoring stations is the same as that of the crack sensors; at the moment, each monitoring point on the earth surface of the monitoring target is provided with 1 satellite positioning monitoring station and 1 crack sensor.
Compared with the prior art, the method and the device realize prediction of deformation displacement by accurately calculating the overall attitude and displacement of the monitored target, and realize comprehensive early warning by fusing multi-source monitoring information. The method has the characteristics of high early warning accuracy, simple algorithm, easy realization and stronger practicability, and can meet the actual requirement of the mountain landslide deformation early warning.
Drawings
Fig. 1 is a schematic block diagram of a landslide deformation comprehensive early warning system.
FIG. 2 is a flow chart of an algorithm for monitoring the attitude of a target.
FIG. 3 is a flow chart of an extended Kalman filter data fusion algorithm for a multi-satellite positioning monitoring station.
Fig. 4 is a flowchart of an iterative gray model deformation prediction algorithm.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail below with reference to the accompanying drawings in conjunction with specific examples.
A landslide deformation comprehensive early warning system is mainly composed of a landslide deformation monitoring network and a comprehensive early warning unit as shown in figure 1.
The landslide deformation monitoring network comprises a satellite positioning landslide deformation monitoring network and a landslide monitoring auxiliary network. The satellite positioning landslide deformation monitoring network consists of 1 satellite positioning reference station and N satellite positioning monitoring stations and is used for monitoring the ground surface displacement of a monitoring target. The satellite positioning reference station is arranged at the position where the foundation is stable, has no signal shielding and no high-power radio emission source, and the N satellite positioning monitoring stations are respectively arranged on different monitoring points of the potential deformation displacement direction of the monitored target. The distance between the satellite positioning reference station and each satellite positioning monitoring station generally cannot exceed 5000 meters. Each satellite positioning monitoring station is connected with a satellite positioning reference station, RTK precision positioning resolving is achieved respectively, and information such as three-dimensional position coordinates, deformation displacement, speed and acceleration of each satellite positioning monitoring station is obtained. In the preferred embodiment of the invention, the satellite positioning landslide deformation monitoring network is a BDS, GPS or GNSS landslide deformation monitoring network. The landslide deformation monitoring auxiliary network consists of a rain gauge and a plurality of crack sensors and is used for respectively collecting local rainfall information and measuring the size of cracks. The rain gauge is arranged in the monitoring target, and the plurality of crack sensors are arranged on different monitoring points of the monitoring target. The number and placement positions of the crack sensors do not necessarily correspond to the number and placement positions of the satellite positioning monitoring stations. However, in the preferred embodiment of the invention, in order to simplify the installation and maintenance of the system, the number of the crack sensors is consistent with that of the satellite positioning monitoring stations, and 1 satellite positioning monitoring station and 1 crack sensor are arranged at each monitoring point of the ground surface of the monitoring target.
The comprehensive early warning of landslide deformation is realized to above-mentioned comprehensive early warning unit, includes wherein: smoothing RTK positioning information data of each monitoring point by using a Kalman filter, eliminating outliers, improving monitoring precision, and extracting information such as effective displacement; resolving the real-time attitude of the monitoring target by utilizing an attitude resolving algorithm, and realizing attitude prediction according to deformation prediction information; the method comprises the steps that an extended Kalman filter is utilized to realize the data fusion of displacement deformation, speed and acceleration of multiple sites in a monitoring network, and the optimal estimation of comprehensive displacement is realized; predicting deformation displacement by using an MGM (1, 1) gray model algorithm; the landslide deformation comprehensive early warning module judges the landslide deformation grade through the fusion of attitude angle information, real-time deformation alarm grade information, prediction deformation alarm grade information, rainfall information, crack size information and the like, so that the landslide deformation comprehensive early warning is realized.
(1) Data fusion algorithm for landslide mass multi-satellite positioning monitoring station
In the actual landslide monitoring, for the same monitoring target, in order to early warn the whole deformation trend of the landslide monitoring target, a landslide monitoring network usually adopts a network arrangement mode that one landslide satellite positioning reference station is matched with a plurality of satellite positioning monitoring stations. During deformation early warning analysis, in order to avoid human errors and other uncertain factors, the displacement monitoring data of a plurality of satellite positioning monitoring stations need to be fused for comprehensive early warning, and the prediction is not only carried out according to a single satellite positioning monitoring station. Therefore, the extended Kalman filter is designed to fuse the positioning data of the multiple satellite positioning monitoring stations, and the global optimal estimation of the overall displacement of the monitored target is realized. The specific design is as follows:
in the stage of monitoring the deformation of the target, a state equation and an observation equation can be established as follows:
XE(k+1)＝FEXE(k)+GEWE(k) (1)
ZE(k+1)＝HEXE(k)+VE(k) (2)
in the formula:
Figure BDA0001572702440000041
k represents the time k, x (k) is the accumulated displacement of the sliding mass, and v (k) is the speed of the sliding mass;
Figure BDA0001572702440000042
t is the sampling time interval of the RTK positioning system;
Figure BDA0001572702440000043
WE(k) a (k) is the acceleration of the sliding mass at two adjacent moments; zE(k) Monitoring the deformation displacement of the target; hE＝[10]；VE(k) Is a Gaussian white noise sequence with variance of sigma2Has a mean value of zero, and is equal to WE(k) Is not relevant. The monitoring object motion model can be described as
Figure BDA0001572702440000051
Assuming that N satellite positioning monitoring stations need to be fused, Z in the global observation equation (2)E(k)、HE(k)、VE(k) Can be described as:
ZE(k)＝[Z1 T(k),Z2 T(k),…ZN T(k)]T(4)
HE(k)＝[H1 T(k),H2 T(k),…HN T(k)]T(5)
VE(k)＝[V1 T(k),V2 T(k),…VN T(k)]T(6)
the specific algorithm for applying the extended Kalman filter algorithm to the model is as follows:
Figure BDA0001572702440000052
Figure BDA0001572702440000053
PE(k+1|k)＝FE(k)PE(k|k)FE T(k)+GE(k)QE(k)GE T(k) (9)
KEi(k+1)＝PE(k+1)HEi T(k+1)REi -1(k+1) (10)
process noise is WECan be provided with WE(k) A (k), solving a process noise covariance matrix QE，QE＝Cov(WE)＝E(WEWT E)，QEIs a symmetric matrix of 2N × 2N. Measuring noise vector as VESolving a process noise covariance matrix RE，RE＝Cov(VE)＝E(VEVT E)，REIs a symmetric matrix of 2N × 2N.
Therefore, the comprehensive displacement information of the monitored target can be obtained by fusing the displacement information data of the multiple monitoring points.
(2) Solving for monitoring target attitude
During deformation monitoring, each attitude angle index of the monitored target can effectively reflect the change condition of the deformation, so that the attitude values of the monitored target, such as a pitch angle, a roll angle, a course angle and the like, can be solved by utilizing the position coordinate change information of the front epoch and the rear epoch of a multi-antenna (multi-monitoring point). Solving for the pose involves two coordinate systems: a geographical coordinate system, a carrier coordinate system. The geographic coordinate system (L system) is based on the phase center of the main antenna as the origin, the north local meridian is the Y axis, the X axis points to the east and is perpendicular to the Y axis, and the Z axis and the X, Y axis form a right-hand coordinate system. The carrier coordinate system (b system) takes the phase center of the main antenna as the origin; the Y axis points to the connecting line direction of the M antenna and the M +1 antenna, the X axis is perpendicular to the Y axis and points to the right side of the monitoring target, and the Z axis and the X, Y axis form a right-hand coordinate system.
As a certain monitoring station is selected as a main monitoring station in the monitoring network, the main monitoring station and other monitoring stations (which are not collinear) can respectively obtain N-1 baselines formed by a geographic coordinate system and a carrier coordinate system according to the basic principle of coordinate system conversion:
Figure BDA0001572702440000054
in the formula αiA base line vector S in a carrier coordinate system for the ith base lineiA baseline vector of the ith baseline in the geographic coordinate system, Cb LIs a coordinate rotation matrix.
Wherein
Figure BDA0001572702440000061
Psi is the heading angle, theta is the pitch angle, and gamma is the roll angle.
Two baseline vectors are arbitrarily chosen in the geographic coordinate system: s12＝[x12y12z12]T、S13＝[x13y13z13]T. The initial attitude angle, the initial values of the pitch angle theta DEG and the heading angle psi DEG can be obtained by defining the attitude angle as follows:
Figure BDA0001572702440000062
ψ°＝-arctan(x12/y12)(13)
the initial value γ ° of the roll angle can be obtained by the following transformation:
Figure BDA0001572702440000063
γ°＝-arctan(z1′3/x1′3) (15) writing equation (11) in error form gives:
Figure BDA0001572702440000064
in the formula, δ ψ, δ θ, δ γ are parameters to be estimated, that is, attitude angle correction numbers. Observed value weight matrix
Figure BDA0001572702440000065
Figure BDA0001572702440000066
Are respectively SiAnd αiCovariance matrix of (2).
The attitude correction number can be estimated according to the principle of the least squares method:
Figure BDA0001572702440000067
therefore, after the attitude angle is estimated by least squares, the course angle, the pitch angle and the roll angle of the monitored target are respectively
Figure BDA0001572702440000068
Figure BDA0001572702440000069
Figure BDA00015727024400000610
(3) Landslide deformation prediction gray MGM (1, 1) model
The optimal length of the data model of the gray prediction model is 40 groups of data, and if the number of the groups is too large, errors are accumulated, and the prediction accuracy is affected. In practical application of BDS/GPS landslide monitoring, one epoch is generally solved once, that is, each epoch obtains a positioning result, and if each positioning result is predictedThe analysis results in excessive calculation and increased data redundancy. The positioning data of each epoch is now sampled with a sampling interval set to Tc，TcCan be set according to actual requirements. Now TcSet to 4 hours. The sampling method is to take the average value of the displacement within a period of time to avoid the influence of factors such as wild value and the like. Observed value array X after sampling(0)The following were used:
X(0)＝{x(0)(1),x(0)(2),…,x(0)(n)} (21)
where n is the length of the observation sequence, i.e., the length of the data model, and is now set to 40; x is the number of(0)(k) For the amount of deformation displacement, in the actual landslide prediction model, since the deformation displacement is generally unidirectional, x(0)(k) Is a non-negative number. According to the basic theory of the gray GM (1, 1) model, the sequence X will now be observed(0)Accumulating once to obtain data generation sequence X(1)The specific formula is as follows:
X(1)＝{x(1)(1),x(1)(2)，…,x(1)(n)} (22)
wherein
Figure BDA0001572702440000071
Generating a sequence X from data(1)Can calculate the adjacent mean value generation sequence Z(1)The concrete formula is as follows:
Z(1)＝{z(1)(1),z(1)(2),z(1)(3),…,z(1)(n)} (23)
wherein z is(1)(k) For whitening background value, two generated data X are generated before and after the value is taken(1)Average value of (d); z is a radical of(1)(k)＝0.5x(1)(k)+0.5x(1)(k-1),(k＝2,3,…,n)；
To X(1)Constructing a simulated white differential equation:
Figure BDA0001572702440000072
the symbols are replaced and arranged to obtain a differential equation of equation (24):
x(0)(k)+aiz(1)(k)＝bi(25)
the formula (25) is the GM (1, 1) model. In the formula, ai、biAll are equations requiring undetermined parameters, aiCan reflect the deformation displacement x(0)The development tendency of (a), called the development coefficient; biThe magnitude of the system effect, which can reflect the grey information coverage, is called the grey effect amount.
The time response functions of the white differential equation and the gray differential equation of the formula (25) are constructed, and the specific formula is as follows:
Figure BDA0001572702440000073
one time of accumulation and subtraction can be carried out for restoring the original number sequence, and formula variables are replaced and sorted after accumulation and subtraction
Figure BDA0001572702440000074
To solve the parameter aiAnd biThe least square estimation can be used to solve the formula (25), which is specifically:
Figure BDA0001572702440000081
wherein the content of the first and second substances,
Figure BDA0001572702440000082
from this, the coefficient of development a can be obtainediAnd biThe amount of gray effect is substituted into the formula (27) to obtain a predicted value at each time.
The prediction is carried out in a sliding window mode, namely the data model length of the gray prediction model is fixed, but the initial observation data is continuously pushed forward, so that the iterative computation of the predicted value is realized. For each initial observation, a set of parameters a can be calculated by equations (22) to (28)iAnd biWhich corresponds to a prediction of
Figure BDA0001572702440000083
The specific formula is as follows.
Figure BDA0001572702440000084
For the prediction of the landslide deformation at the same moment, a plurality of predicted values can be obtained through iterative prediction. These predicted values
Figure BDA0001572702440000085
Fluctuation within a certain range, and predicted value fitting can be carried out in a weighting mode. If the number of the current predicted values is set as MnowThe number of the prediction points at the end of the sliding window of the prediction model is MendIf the number of the predicted values of the gray model is Mall：
Mall＝Mnow-Mend(30) When M isallGreater than a certain value EnumAfter the meeting, the user can use the device,
Figure BDA0001572702440000086
accuracy may be degraded continuously, so MallWithin a certain range
Figure BDA0001572702440000087
It is effective.
Figure BDA0001572702440000088
In the formula, wiIs a predicted value
Figure BDA0001572702440000089
The corresponding weight.
Specifically, the landslide deformation comprehensive early warning method realized by the device specifically comprises the following steps:
step 1, collecting the current three-dimensional position [ x ] obtained by resolving each satellite positioning monitoring station in the landslide deformation monitoring network by using an RTK positioning algorithmL(k),yL(k),zL(k)]NAnd storing the data in a specific data file for later use in comprehensive early warning analysis. Reading the position information of the satellite positioning monitoring station of the data file, firstly positioning the initial position [ x (0), y (0), z (0) of the monitoring station according to each satellite]NSolving the deformation displacement x corresponding to each timebin(k)N(ii) a Solving the deformation speed and the deformation acceleration of the landslide body according to the sampling time interval and the displacement difference of the deformation displacement amount at the front moment and the rear moment; the data information is mm-level information.
Step 2, according to the collected deformation displacement x of a plurality of satellite positioning monitoring stationsbin(k)NUsing Kalman filter pair xbin(k)NAnd carrying out smooth filtering, eliminating or repairing abnormal monitoring values and improving the precision of monitoring data.
Using Kalman filter pair xbin(k)NSmoothing is performed. For one of xbin(k) Comprises the following steps: quantity of state xp(k) And an observation vector xbin(k) The relationship (measurement equation) of (1) is: x is the number ofpin(k)＝Cpxbin(k)+vbin(k) (ii) a The state equation is: x is the number ofp(k)＝Apxp(k-1)+Bpu (k-1), wherein the relation matrix of the measurement equation is: cp＝[10](ii) a The state transition matrix and the input relation matrix of the state equation are respectively as follows:
Figure BDA0001572702440000091
Figure BDA0001572702440000092
implementation of x using the basic principles of KF time update and measurement updatebin(k) Wherein the noise Q is reduced by adjusting the processpAnd measuring the noise RpRealizing optimal filtering result to obtain deformation displacement xpin(k) And eliminating outliers.
Meanwhile, solving the deformation speed v of the landslide body according to the front-back displacement difference and the sampling time intervalpin(k)N＝[vpin(k),vpin(k),vpin(k)]NDeformation acceleration apin(k)N＝[apin(k),apin(k),apin(k)]NAnd (4) equivalence.
Step 3, establishing a baseline vector S according to deformation monitoring position information of any N '(N' is not less than 3) non-collinear satellite positioning monitoring stations in the monitoring networkiAnd calculating a monitoring target course angle psi, a pitch angle theta and a roll angle gamma by using a least square attitude calculation algorithm.
Extracting three-dimensional position coordinates x of each smoothed satellite positioning monitoring stationpin(k)＝(xk,yk,zk)pinAnd the like, and sampling is realized according to requirements so as to reduce the calculation amount of data. And performing time synchronization matching on the data after each monitoring sampling, judging whether the time synchronization is realized, if not, continuing the matching, and if the synchronization is realized, entering the next step. Selecting a certain monitoring station as a main monitoring station, and constructing a baseline vector S with the positions of other non-collinear monitoring stationsiThe transformation of the coordinate system is achieved by the following equation.
Figure BDA0001572702440000093
In the formula αiA base line vector S in a carrier coordinate system for the ith base lineiA baseline vector of the ith baseline in the geographic coordinate system, Cb LIs a coordinate rotation matrix.
Wherein the content of the first and second substances,
Figure BDA0001572702440000101
psi is the heading angle, theta is the pitch angle, and gamma is the roll angle.
Arbitrarily selecting two baseline vectors S in a geographic coordinate system12、S13And resolving initial values theta, psi and gamma of the pitch angle, the course angle and the roll angle. Constructing an attitude coefficient matrix Ge、LeThe matrix is equalized, and the corresponding weight P of the base line is calculatede。
Wherein the content of the first and second substances,
Figure BDA0001572702440000102
Figure BDA0001572702440000103
Figure BDA0001572702440000104
are respectively SiAnd αiCovariance matrix of (2).
The attitude correction numbers can be settled according to the principle of the least square method:
Figure BDA0001572702440000105
and resolving the attitude of the monitoring target according to the attitude correction number and initial values of the pitch angle, the course angle and the roll angle.
In the invention, the attitude angle of the monitoring target is calculated by utilizing a least square attitude calculation algorithm, and a specific algorithm flow chart is shown in FIG. 2. The specific implementation steps are as follows:
step 3.1: extracting three-dimensional position coordinates x of each smoothed satellite positioning monitoring stationpin(k)＝(xk,yk,zk)pinAnd the like, and sampling is realized according to requirements so as to reduce the calculation amount of data.
Step 3.2: and performing time synchronization matching on the data after each monitoring sampling in the previous step, judging whether the time synchronization is realized, if not, continuing the matching, and if the synchronization is realized, entering the next step.
Step 3.3: selecting a certain monitoring station as a main monitoring station, and constructing a baseline vector S with the positions of other non-collinear monitoring stationsi。
Step 3.4: the baseline vector of the previous step is in a geographic coordinate system, and the transformation of the coordinate system is realized by using the following formula.
Figure BDA0001572702440000106
In the formula αiA base line vector of the ith base line in the carrier coordinate system,SiA baseline vector of the ith baseline in the geographic coordinate system, Cb LIs a coordinate rotation matrix.
Wherein the content of the first and second substances,
Figure BDA0001572702440000111
psi is the heading angle, theta is the pitch angle, and gamma is the roll angle.
Step 3.5: arbitrarily selecting two baseline vectors S in a geographic coordinate system12、S13And resolving initial values theta, psi and gamma of the pitch angle, the course angle and the roll angle.
Assume a baseline vector of S12＝[x12y12z12]T、S13＝[x13y13z13]T. The initial attitude angle, the initial values of the pitch angle theta DEG and the heading angle psi DEG can be obtained by defining the attitude angle as follows:
Figure BDA0001572702440000112
ψ°＝-arctan(x12/y12)
the initial value γ ° of the roll angle can be obtained by the following transformation:
Figure BDA0001572702440000113
γ°＝-arctan(z1′3/x1′3)
step 3.6: constructing an attitude coefficient matrix Ge、LeMatrix and calculating a baseline correspondence weight Pe。
Figure BDA0001572702440000114
Figure BDA0001572702440000115
Figure BDA0001572702440000116
Are respectively SiAnd αiCovariance matrix of (2).
Step 3.7: the attitude correction number can be solved according to the principle of the least square method:
Figure BDA0001572702440000117
step 3.8: and resolving the attitude of the monitoring target according to the attitude correction number and initial values of the pitch angle, the course angle and the roll angle.
Step 4, extracting the position x of each satellite positioning monitoring station after smoothingpin(k) Velocity vpin(k) Acceleration apin(k) Performing time synchronization alignment of the satellite positioning monitoring station to construct an extended Kalman filter state transition matrix GEMeasuring relation matrix HEAnd the EKF is utilized to fuse data of a plurality of satellite positioning monitoring stations in the monitoring network, so that the optimal estimation of the comprehensive displacement information of the monitoring target is realized.
Extracting the smoothed position x of each satellite positioning monitoring stationpin(k) Velocity vpin(k) Acceleration apin(k) And performing time synchronization alignment of the satellite positioning monitoring stations, and judging whether the time of each satellite positioning monitoring station is aligned. Constructing a state quantity matrix XE(k)：
Figure BDA0001572702440000121
k denotes the time k, xpin(k) Integrating the displacement, v, for the sliding masspin(k) Is the speed of the landslide mass. XE(k) Expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations;
construction of a Global Observation matrix ZE(k) And measuring the relation matrix HE(k) The method comprises the following steps Wherein Z isE(k) For monitoring the amount of deformation displacement of the target, HE(k)＝[10]Expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations;
constructing a State transition matrix FE(k) And matrix GE(k)：
Figure BDA0001572702440000122
T is the sampling time interval of the RTK positioning system;
Figure BDA0001572702440000123
expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations;
process noise is WECan be provided with WE(k)＝a(k)pinSolving the process noise covariance matrix QE，QE＝Cov(WE)＝E(WEWT E)，QEIs a symmetric matrix of 2N × 2N. Measuring noise vector as VESolving a process noise covariance matrix RE，RE＝Cov(VE)＝E(VEVT E)，REIs a symmetric matrix of 2N × 2N.
Updating time, and predicting the current system state; and (5) measuring and updating, and correcting the priori estimation value of the prediction state to obtain the global optimal estimation of the overall displacement of the monitored target.
In the invention, the global optimal estimation of the whole displacement of the monitoring target is realized by utilizing a landslide mass multi-satellite positioning monitoring station data fusion algorithm, and a specific algorithm flow chart is shown in figure 3. The specific implementation steps are as follows:
step 4.1: extracting the smoothed position x of each satellite positioning monitoring stationpin(k) Velocity vpin(k) Acceleration apin(k) And performing time synchronization alignment of the satellite positioning monitoring stations, and judging whether the time of each satellite positioning monitoring station is aligned.
Step 4.2: constructing a state quantity matrix XE(k)：
Figure BDA0001572702440000124
k denotes the time k, xpin(k) Integrating the displacement, v, for the sliding masspin(k) Is the speed of the landslide mass. XE(k) Expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations;
step 4.3: construction of a Global Observation matrix ZE(k) And measuring the relation matrix HE(k) The method comprises the following steps Wherein Z isE(k) For monitoring the targetAmount of deformation displacement, HE(k)＝[10]Expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations;
step 4.4: constructing a State transition matrix FE(k) And matrix GE(k)：
Figure BDA0001572702440000125
Figure BDA0001572702440000126
Expanding the dimension according to the number N of the actual fusion satellite positioning monitoring stations; wherein T is the sampling time interval of the RTK positioning system;
step 4.5: constructing a process noise covariance matrix QEMeasuring the covariance matrix RE: process noise is WECan be provided with WE(k)＝a(k)pinSolving the process noise covariance matrix QE，QE＝Cov(WE)＝E(WEWT E)，QEIs a symmetric matrix of 2N × 2N. Measuring noise vector as VESolving a process noise covariance matrix RE，RE＝Cov(VE)＝E(VEVT E)，REIs a symmetric matrix of 2N × 2N.
Step 4.6: updating time, and predicting the current system state; measuring and updating, and correcting the priori estimated value of the prediction state to obtain the global optimal estimation x of the overall displacement of the monitored targetall(k)。
Step 5, synthesizing displacement information x of the monitored targetall(k) Sampling at equal time intervals, setting the sampling time interval as T, setting the initial epoch and the length n of a prediction sliding window, and monitoring the target landslide deformation by utilizing the corresponding principle of an iterative gray model
Figure BDA0001572702440000131
And (6) predicting.
Setting a prediction sliding window initial epoch, and carrying out comprehensive displacement information x on a monitored targetall(k) Sampling is carried out, wherein sampling intervals are 4 hours, 40 groups are sampled in total, and the sampling method is to take the average value of displacement in a period of time.The sampled data are:
X(0)＝{x(0)(1),x(0)(2),…,x(0)(n)}
will observe the sequence X(0)Accumulating once to obtain data generation sequence X(1)：
X(1)＝{x(1)(1),x(1)(2)，…,x(1)(n)}
Wherein
Figure BDA0001572702440000132
Generating a sequence X from data(1)Can calculate the adjacent mean value generation sequence Z(1)：
Z(1)＝{z(1)(1),z(1)(2),z(1)(3),…,z(1)(n)}
Wherein z is(1)(k) For whitening background value, two generated data X are generated before and after the value is taken(1)Average value of (a). Constructing a Jacobian matrix B:
Figure BDA0001572702440000133
estimating the parameter coefficients using least squares:
Figure BDA0001572702440000134
and (3) carrying out once subtraction to restore the original number sequence:
Figure BDA0001572702440000135
the initial observation data is continuously pushed forward, the iterative computation of the predicted value is realized, and N is obtainedwinGroup parameter aiAnd biCorresponding to a predicted value of
Figure BDA0001572702440000141
Deviation fromCalculating the ratio of the average values to each predicted value
Figure BDA0001572702440000142
Weight w ofiAnd is and
Figure BDA0001572702440000143
wherein E isnumThe number of effective predicted values.
Calculating a predicted value according to the weight:
Figure BDA0001572702440000144
in the invention, the prediction of the monitored target deformation displacement is realized by using a landslide deformation prediction gray MGM (1, 1) model, and a specific algorithm flow chart is shown in FIG. 4. The specific implementation steps are as follows:
step 5.1: setting a prediction sliding window initial epoch, and carrying out comprehensive displacement information x on a monitored targetall(k) Sampling is carried out, wherein sampling intervals are 4 hours, 40 groups are sampled in total, and the sampling method is to take the average value of displacement in a period of time. The sampled data are:
X(0)＝{x(0)(1),x(0)(2),…,x(0)(n)}
step 5.2: will observe the sequence X(0)Accumulating once to obtain data generation sequence X(1)：
X(1)＝{x(1)(1),x(1)(2)，…,x(1)(n)}
Wherein
Figure BDA0001572702440000145
Step 5.3: generating a sequence X from data(1)Can calculate the adjacent mean value generation sequence Z(1)：
Z(1)＝{z(1)(1),z(1)(2),z(1)(3),…,z(1)(n)}
Wherein z is(1)(k) For whitening background value, taking the value of anteroposteriorGenerating data X in two(1)Average value of (a).
Step 5.4: constructing a Jacobian matrix B:
Figure BDA0001572702440000146
step 5.5: estimating the parameter coefficients using least squares:
Figure BDA0001572702440000147
wherein:
Figure BDA0001572702440000148
aito develop the coefficient, biThe grey effect is indicated.
Step 5.6: and (5) carrying out once subtraction to restore the original sequence.
The first data point of the sequence is now taken as an example for explanation:
Figure BDA0001572702440000151
step 5.7: the initial observation data is continuously pushed forward to realize the iterative calculation of the predicted value to obtain a plurality of groups of parameters aiAnd biCorresponding to a predicted value of
Figure BDA0001572702440000152
Step 5.8: calculating each predicted value
Figure BDA0001572702440000153
Weight w ofiAnd calculating a predicted value according to the weight:
Figure BDA0001572702440000154
step 6, carrying out actual measurement early warning on the current state of the monitoring target according to a preset early warning rule by using the actual measurement values of the deformation displacement, the deformation speed and the deformation acceleration calculated in the steps 3 and 4; and meanwhile, predicting and early warning the subsequent state of the monitored target according to a preset early warning rule by using the predicted values of the deformation displacement, the deformation speed and the deformation acceleration predicted in the step 5.
The preset early warning rules are artificially set according to specific conditions, and if the preset early warning rules are combined in a mode of weighting the deformation displacement, the deformation speed and the deformation acceleration, early warning judgment is carried out, and the like can also be carried out after the preset early warning rules are combined with the sum or the sum of the deformation displacement, the deformation speed and the deformation acceleration.
If using the integrated displacement X of the monitored targetall(k) And monitoring target attitude for real-time early warning grade GlealdmJudging that the horizontal displacement of the earth surface exceeds 15mm and the vertical displacement exceeds 20mm, and judging the earth surface to be a blue grade; the horizontal displacement of the earth surface exceeds 30mm, the vertical displacement exceeds 35mm, and the earth surface can be judged to be yellow grade; the horizontal displacement of the earth surface exceeds 40mm, the vertical displacement exceeds 45mm, and the orange grade can be judged; the horizontal displacement of the earth surface exceeds 50mm, the vertical displacement exceeds 55mm, and the red grade can be judged. The threshold for actual grading must be set according to geological conditions.
E.g. using deformation prediction
Figure BDA0001572702440000155
Grade judgment is carried out on speed, acceleration and the like, and grade G of displacement and posture is predictedpremThe judgment is the same as the real-time displacement and posture grade judgment. Velocity tangent angle and acceleration early warning grade GprevaWhen the acceleration is less than 0 and the speed tangent angle β is less than 45 degrees, G can be judgedprevaIs a blue color grade; oscillating when the acceleration is equal to 0 or near 0 value, and the speed tangent angle is 45 °<β<At 80 deg.C, G can be judgedprevaYellow grade, when the acceleration is larger positive value and the speed tangent angle is between 80 and β<At 85 deg.C, G can be judgedprevaOrange grade, when the acceleration is larger positive value and increased, the speed tangent angle is between 85 degrees and β degrees, and the speed slides downwards at about 89 degrees, G can be judgedprevaIs a red grade;
if the subsequent state early warning judgment is carried out in a weighting mode, the early warning level weight of the predicted displacement is wpremPredicting the velocity tangent angle and the early warning level weight of the acceleration as wprevaThen monitoring the early warning level G of the subsequent state prediction of the targetpregradeComprises the following steps:
Gpregrade＝wpremGprem+wprevaGpreva
wherein G ispremFor predicting displacement and attitude warning levels, GprevaThe velocity tangent angle and the early warning level of the acceleration are obtained.
And 7, collecting rainfall and/or crack displacement by a crack sensor by a rain gauge in the landslide deformation monitoring network, and performing comprehensive prediction and early warning on a monitored target according to a preset comprehensive early warning rule by combining the deformation displacement, the deformation speed and the measured value and the predicted value of the deformation acceleration.
And (4) utilizing landslide factors such as rainfall, crack displacement and the like to judge the early warning grade, wherein a preset early warning rule is manually set according to specific conditions. In the invention, the rainfall is measured by a rain gauge, and the rainfall L is collectedrInformation, rainfall landslide warning factor wrAnd grade GrAnd days of rainfall drHave close relationship. Grade GrThe higher, wrThe larger. The specific calculation relation is determined according to the local rainfall and the landslide number statistical rate. Statistics of days of rainfall drAnalysis of landslide and days of rainfall drBetween landslide early warning factor wrd：
Figure BDA0001572702440000161
Wherein: setting the historical landslide number of the surrounding area of the monitoring area as x, xiFor the ith value of x, the value of x,
Figure BDA0001572702440000162
is the mean of x; the rainfall 7 days before the monitoring time is set as y, yiFor the ith value of y, the value of y,
Figure BDA0001572702440000163
is the mean value of y.
Early warning grade G of rainfallrainThe judgment is as follows: when w isrd<0.3 and LrG when the thickness is less than or equal to 20mmrainCan be judged as blue when the color is 0.3<wrd<0.5 and LrG when the thickness is less than or equal to 30mmrainCan be judged as yellow grade; when 0.5<wrd<0.8 and LrG when the thickness is less than or equal to 40mmrainCan be judged as orange grade; when w isrd>0.8 and LrWhen the thickness is less than or equal to 90mm, GrainA red color level may be determined.
And selecting a silicon micro-inclinometer which can be buried under the ground surface crack of the monitored target for a long time. Calculating the variation of the target inclination angle according to the reading of the silicon micro-inclinometer:
△θ＝K1Fi+K2F2 i+K3Fi 3-(K1F0+K2F0 2+K3F0 3)
wherein △ theta is the variation of the target inclination angle to be monitored, KiThe parameter is a cubic fitting parameter of the inclinometer, i is 1,2 and 3, and i is a power term; fiMeasuring a reading for the inclinometer; f0Is an inclinometer reference reading. Calculating the horizontal displacement:
△x(k)＝Lsin△θ
Figure BDA0001572702440000164
wherein L is the distance between the guide wheels, xlAnd d is the horizontal distance of the crack and the measuring times of the inclinometer. Calculating the vertical displacement:
△z(k)＝Lcos△θ
Figure BDA0001572702440000165
wherein z islIs a vertical displacement.
Early warning grade G of internal displacement of ground surface crackfissureThe judgment is as follows: when the internal displacement x of the ground surface cracklLess than or equal to 10mm and ZlWhen the thickness is less than or equal to 15mm, GfissureCan be judged as blue grade when 10<xlLess than or equal to 20mm and 15mm<ZlG when the thickness is less than or equal to 25mmfissureCan be judged as yellow grade; when 20 is turned on<xlLess than or equal to 40mm and 25<Zl≤45mm，GfissureCan be judged as orange grade; when x isl>40mm and Zl>At 45mm time, GfissureA red color level may be determined. The threshold for actual grading must be adjusted according to geological conditions.
And 8, carrying out grade judgment of real-time comprehensive early warning and subsequent state comprehensive early warning according to each early warning grade, and carrying out comprehensive early warning by using a set fault tolerance mechanism.
From the step 6, the comprehensive displacement X of the target is monitoredall(k) And monitoring target attitude real-time early warning grade GlealdmAnd step 7, rainfall landslide early warning grade GrainAnd early warning grade G of internal displacement of ground surface crackfissureWeighted resolving to obtain a real-time comprehensive early warning grade GlealgradeComprises the following steps:
Glealgrade＝wlealdmGlealdm+wrain1Grain+wfissure1Gfissure
wherein, wlealdmIs GlealdmWeight factor of, wrain1Weight factor, w, for rainfall landslide warning level in real timefissure1And weighting factors of the early warning grade of the internal displacement of the ground surface fracture in a real-time state.
Early warning level G of subsequent state prediction of the monitored target by step 6pregradeAnd step 7, rainfall landslide early warning grade GrainAnd early warning grade G of internal displacement of ground surface crackfissurePerforming weighted calculation to obtain a comprehensive early warning grade G of the subsequent state prediction of the real monitoring targetpgComprises the following steps:
Gpg＝wpregradeGpregrade+wrain2Grain+wfissure2Gfissure
wherein, wpregradeIs GpregradeWeight factor of, wrain2Weight factor, w, for the subsequent state prediction of rainfall landslide warning levelfissure2And the weight factor is used for predicting the subsequent state of the early warning level of the internal displacement of the earth surface fracture.
The fault tolerance mechanism is set manually. The fault tolerance mechanism is set as follows: when the early warning grade changes at a certain moment, storing the early warning result, inquiring whether the records of a plurality of moments before and after the early warning result meet the requirement of the early warning grade changes, if the records meet the requirement after continuous monitoring, carrying out early warning, and if the records do not meet the requirement, judging again.
The method applies an RTK multimode high-precision positioning technology to landslide deformation early warning, and carries out comprehensive early warning according to factors such as displacement, acceleration, speed, local rainfall and the like of landslide body (monitoring target) deformation. The threshold value of the early warning level is different due to different geology, rainfall, monitoring target postures and the like. And by utilizing a prediction algorithm, calculating the future deformation trend of the landslide monitoring target according to the displacement deformation quantity of the monitoring target, and combining data such as a crack meter and the like to realize early warning of landslide occurrence. The early warning grade can be divided into four grades of blue, yellow, orange and red, and the landslide forming degree is gradually increased.
It should be noted that, although the above-mentioned embodiments of the present invention are illustrative, the present invention is not limited thereto, and thus the present invention is not limited to the above-mentioned embodiments. Other embodiments, which can be made by those skilled in the art in light of the teachings of the present invention, are considered to be within the scope of the present invention without departing from its principles.

Claims (9)

1. A landslide deformation comprehensive early warning method is characterized by comprising the following steps:
step 1, solving three-dimensional position coordinates of monitoring points of monitoring targets of landslide deformation monitoring networks by each satellite positioning monitoring station, and solving displacement according to initial positions of the satellite positioning monitoring stations;
step 2, smoothing the displacement obtained by each satellite positioning monitoring station by using a Kalman filter to obtain an actual measurement value of the deformation displacement; according to the sampling time interval and the displacement difference of the measured values of the deformation displacement at the front and rear moments, the measured values of the deformation speed and the deformation acceleration of the monitoring target are solved;
step 3, extracting the measured values of the deformation displacement of each monitoring point obtained in the step 2, and solving the attitude of the monitoring target by using a least square method, namely the measured values of the course angle, the pitch angle and the roll angle of the monitoring target;
step 4, extracting the measured values of the deformation displacement, the deformation speed and the deformation acceleration of each monitoring point obtained in the step 2, and solving the measured value of the comprehensive displacement of the monitoring target by using an extended Kalman filter;
step 5, setting an initial epoch and a length of a prediction sliding window, sampling the comprehensive displacement of the monitored target at equal time intervals, utilizing an iterative gray model to realize a predicted value of the deformation displacement of the monitored target, and calculating predicted values of the deformation speed and the deformation acceleration according to the sampling time interval and the displacement difference of the predicted values of the deformation displacement at the previous moment and the next moment;
step 6, carrying out actual measurement landslide grade early warning on the current state of the monitoring target according to a preset actual measurement early warning rule by using the actual measurement values of the deformation displacement, the deformation speed and the deformation acceleration calculated in the step 2; and meanwhile, carrying out prediction landslide grade early warning on the subsequent state of the monitored target according to a preset prediction early warning rule by using the predicted values of the deformation displacement, the deformation speed and the deformation acceleration predicted in the step 5.
2. The landslide deformation comprehensive early warning method according to claim 1, wherein the specific process of the step 5 is as follows:
step 5.1, setting a prediction sliding window initial epoch, and sampling a measured value of the comprehensive displacement of the monitoring target to obtain an observation sequence of the current prediction sliding window;
step 5.2, accumulating the observation sequence for one time to obtain a data generation sequence, and calculating an adjacent mean value generation sequence according to the data generation sequence;
step 5.3, constructing a Jacobian matrix by utilizing the adjacent mean value generation sequence, and constructing a least square estimation parameter coefficient, namely a development coefficient and a gray action amount according to the Jacobian matrix;
step 5.4, carrying out one-time subtraction on the observation number series by using the development coefficient and the gray effect quantity to obtain a prediction number series of the current prediction sliding window;
step 5.5, continuously pushing the initial observation data forward by using the prediction sliding window, and performing iterative computation by repeating the step 5.2-5.4 to obtain the prediction sequence of each prediction sliding window;
step 5.6, performing weighted accumulation on the prediction series of the prediction sliding windows to obtain a prediction value of the deformation displacement of the monitoring target;
and 5.7, calculating the predicted values of the deformation speed and the deformation acceleration according to the sampling time interval and the displacement difference of the predicted values of the deformation displacement at the previous moment and the next moment.
3. The landslide deformation comprehensive early warning method according to claim 1, wherein the specific process of step 3 is as follows:
step 3.1, extracting the measured values of the deformation displacement of each monitoring point obtained in the step 2;
3.2, establishing a baseline vector by using the three-dimensional position coordinates of any more than 3 non-collinear satellite positioning monitoring stations in the landslide deformation monitoring network, and randomly selecting 2 baseline vectors to calculate the initial values of a pitch angle, a course angle and a roll angle;
3.3, calculating corresponding weight of a base line according to the constructed attitude coefficient matrix, and estimating an attitude correction number by using a least square method;
and 3.4, calculating the attitude of the monitoring target according to the attitude correction number and the initial values of the pitch angle, the course angle and the roll angle, namely the actual measured values of the course angle, the pitch angle and the roll angle of the monitoring target.
4. The landslide deformation comprehensive early warning method according to claim 1, wherein the specific process of the step 4 is as follows:
step 4.1, extracting measured values of deformation displacement, deformation speed and deformation acceleration of each monitoring point obtained in the step 2;
step 4.1, constructing a state quantity matrix, a global observation matrix and a measurement relation matrix of the extended Kalman filter by using the deformation displacement and the deformation speed, constructing a state transition matrix of the extended Kalman filter by using a sampling time interval, and resolving a process noise covariance matrix and a process noise covariance matrix of the extended Kalman filter according to the deformation acceleration;
and 4.1, updating time by using the extended Kalman filter, predicting the current system state, correcting the priori estimation value of the predicted state through measurement updating, and obtaining the global optimal estimation of the overall displacement of the monitoring target, namely the measured value of the comprehensive displacement of the monitoring target.
5. The landslide deformation comprehensive early warning method according to claim 1, further comprising the steps of:
and 7, collecting rainfall by a rainfall gauge in the landslide deformation monitoring network, collecting the rainfall and collecting the crack displacement by a crack sensor, performing weighted operation on the rainfall and the crack displacement and the measured values and the predicted values of the deformation displacement, the deformation speed and the deformation acceleration obtained in the steps 4 and 5, and performing comprehensive landslide grade early warning on the monitored target according to a preset comprehensive early warning rule.
6. A landslide deformation comprehensive early warning system for realizing the method of claim 1, which is characterized by comprising a landslide deformation monitoring network and a comprehensive early warning unit; the landslide deformation monitoring network comprises a satellite positioning reference station and more than 3 satellite positioning monitoring stations; more than 3 satellite positioning monitoring stations are distributed at different monitoring points on the earth surface of the monitoring target; each satellite positioning monitoring station is respectively connected with a satellite positioning reference station to realize RTK positioning resolving; and the positioning results of all the satellite positioning monitoring stations are transmitted to the comprehensive early warning unit.
7. The landslide deformation comprehensive early warning system according to claim 6, wherein the satellite positioning monitoring station is a GPS satellite positioning receiver, a BDS satellite positioning receiver and/or a GNSS satellite positioning receiver.
8. The landslide deformation comprehensive early warning system according to claim 6, wherein the landslide deformation monitoring net further comprises a rain gauge and more than 2 crack sensors; the rain gauges are distributed on the ground surface of the monitoring target; more than 2 crack sensors are distributed at different monitoring points on the ground surface of the monitoring target; and the monitoring data of the rain gauge and all the crack sensors are transmitted to the comprehensive early warning unit.
9. The landslide deformation comprehensive early warning system according to claim 8, wherein the number of satellite positioning monitoring stations is the same as the number of crack sensors; at the moment, each monitoring point on the earth surface of the monitoring target is provided with 1 satellite positioning monitoring station and 1 crack sensor.
CN201810123208.0A 2018-02-07 2018-02-07 Landslide deformation comprehensive early warning method and system Active CN108332649B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810123208.0A CN108332649B (en) 2018-02-07 2018-02-07 Landslide deformation comprehensive early warning method and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810123208.0A CN108332649B (en) 2018-02-07 2018-02-07 Landslide deformation comprehensive early warning method and system

Publications (2)

Publication Number Publication Date
CN108332649A CN108332649A (en) 2018-07-27
CN108332649B true CN108332649B (en) 2020-04-24

Family

ID=62928630

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810123208.0A Active CN108332649B (en) 2018-02-07 2018-02-07 Landslide deformation comprehensive early warning method and system

Country Status (1)

Country Link
CN (1) CN108332649B (en)

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109443188B (en) * 2018-09-29 2020-05-22 桂林电子科技大学 Double-layer multi-dimensional landslide monitoring method
CN109373911B (en) * 2018-11-02 2020-02-14 中国地质科学院地质力学研究所 Ground surface displacement gridding dynamic monitoring method
CN109556897A (en) * 2018-11-16 2019-04-02 王玉波 A kind of bridge construction system in science of bridge building field
CN109696153B (en) * 2018-12-25 2021-09-14 广州市中海达测绘仪器有限公司 RTK tilt measurement accuracy detection method and system
CN110660196A (en) * 2019-11-25 2020-01-07 中国地质调查局水文地质环境地质调查中心 High-steep karst mountain deformation monitoring device and method based on LoRa ad hoc network
CN110986747B (en) * 2019-12-20 2021-03-19 桂林电子科技大学 Landslide displacement combined prediction method and system
CN111022124B (en) * 2019-12-31 2021-04-30 山东交通学院 Advanced early warning method for short-term and long-term deformation of bridge and tunnel engineering
CN111561917B (en) * 2020-03-30 2021-10-26 同济大学 Road side slope monitoring system
CN111694021B (en) * 2020-05-26 2022-08-05 中国科学院国家授时中心 GNSS environment model-based single-station landslide deformation monitoring and early warning method
CN111931345B (en) * 2020-07-09 2021-11-02 西南交通大学 Monitoring data prediction method, device, equipment and readable storage medium
CN112085731A (en) * 2020-09-18 2020-12-15 深圳市易图资讯股份有限公司 Security early warning method, device and equipment based on satellite map and storage medium
CN112581725B (en) * 2020-12-08 2022-07-22 重庆邮电大学 Mountain landslide early warning monitoring system based on NBIOT and LoRa dual-mode communication
CN112924990B (en) * 2021-01-25 2024-03-22 智连空间测绘技术(苏州)有限公司 Landslide body monitoring method and system based on GNSS accelerometer fusion
CN113281742B (en) * 2021-06-02 2023-07-25 西南交通大学 SAR landslide early warning method based on landslide deformation information and meteorological data
CN113421404B (en) * 2021-06-16 2022-04-15 深圳防灾减灾技术研究院 Satellite-ground cooperative slope multi-risk factor combined real-time monitoring and early warning method
CN113570826A (en) * 2021-07-15 2021-10-29 长视科技股份有限公司 Method and system for realizing disaster early warning by river landslide deformation recognition
CN113781745B (en) * 2021-08-20 2023-01-31 合肥星北航测信息科技有限公司 Beidou and micromotion landslide early warning method based on K-means clustering algorithm
CN114236585B (en) * 2021-12-09 2023-04-14 国网思极位置服务有限公司 Target motion monitoring method based on Beidou navigation satellite system and storage medium
CN114333245A (en) * 2022-01-05 2022-04-12 中国地质科学院探矿工艺研究所 Multi-parameter model dynamic early warning method based on landslide three-dimensional monitoring
CN114894082A (en) * 2022-04-21 2022-08-12 北方雷科(安徽)科技有限公司 Mine-based slope deformation landslide early warning method
CN115325928B (en) * 2022-10-17 2023-01-31 西北大学 Landslide earth surface crack integrated monitoring system based on Beidou communication
CN115421172B (en) * 2022-11-04 2023-03-24 南京市计量监督检测院 Beidou deformation monitoring method based on real-time and quasi-real-time combination
CN115638767B (en) * 2022-11-07 2023-10-03 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) Ground subsidence monitoring method
CN116013035A (en) * 2023-02-02 2023-04-25 中铁第一勘察设计院集团有限公司 Early warning method and device for collapse geological disaster and electronic equipment
CN115938095B (en) * 2023-02-22 2023-05-30 湖北通达数科科技有限公司 Landslide monitoring and early warning method and system based on integrated fusion model
CN116739184B (en) * 2023-08-08 2023-11-07 四川川核地质工程有限公司 Landslide prediction method and system

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2746811A2 (en) * 2012-12-18 2014-06-25 Trimble Navigation Limited Methods for generating accuracy information on an ionosphere model for satellite navigation applications
CN104102853A (en) * 2014-08-08 2014-10-15 武汉理工大学 Slope displacement fractal forecasting method improved by grey theory
CN205066775U (en) * 2015-09-17 2016-03-02 泉州装备制造研究所 High accuracy movement track detection device
CN106441174A (en) * 2016-09-09 2017-02-22 桂林电子科技大学 High slope deformation monitoring method and system
KR20170031829A (en) * 2015-09-11 2017-03-22 금호마린테크 (주) Method and system for target acquisition and tracking using marine radar
CN106529185A (en) * 2016-11-24 2017-03-22 西安科技大学 Historic building displacement combined prediction method and system
CN106595566A (en) * 2016-11-21 2017-04-26 山东康威通信技术股份有限公司 Displacement and attitude analysis method for rubber dams
CN107292410A (en) * 2016-03-30 2017-10-24 中国石油天然气股份有限公司 Tunnel deformation Forecasting Methodology and device

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2746811A2 (en) * 2012-12-18 2014-06-25 Trimble Navigation Limited Methods for generating accuracy information on an ionosphere model for satellite navigation applications
CN104102853A (en) * 2014-08-08 2014-10-15 武汉理工大学 Slope displacement fractal forecasting method improved by grey theory
KR20170031829A (en) * 2015-09-11 2017-03-22 금호마린테크 (주) Method and system for target acquisition and tracking using marine radar
CN205066775U (en) * 2015-09-17 2016-03-02 泉州装备制造研究所 High accuracy movement track detection device
CN107292410A (en) * 2016-03-30 2017-10-24 中国石油天然气股份有限公司 Tunnel deformation Forecasting Methodology and device
CN106441174A (en) * 2016-09-09 2017-02-22 桂林电子科技大学 High slope deformation monitoring method and system
CN106595566A (en) * 2016-11-21 2017-04-26 山东康威通信技术股份有限公司 Displacement and attitude analysis method for rubber dams
CN106529185A (en) * 2016-11-24 2017-03-22 西安科技大学 Historic building displacement combined prediction method and system

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
GNSS 动态变形监测系统的设计与初步实现;陈聪等;《四川地震》;20160630(第2期);第31-33页、第47页 *
kinematic landslide monitoring with kalman filtering;M. Acar等;《Natural Hazards and Earth System Sciences》;20081231;第213-221页 *
基于Kalman滤波数据融合技术的滑坡变形分析与预测;刘超云等;《中国地质灾害与防治学报》;20151231;第26卷(第4期);第30-35页、第42页 *
基于卡尔曼滤波的灰色马尔科夫组合模型在基坑变形监测中的应用;朱军桃等;《桂林理工大学学报》;20171130;第37卷(第4期);第653-657页 *

Also Published As

Publication number Publication date
CN108332649A (en) 2018-07-27

Similar Documents

Publication Publication Date Title
CN108332649B (en) Landslide deformation comprehensive early warning method and system
CN104215249B (en) Smoothening method of driving track
CN107218923A (en) Surrounding enviroment history settles methods of risk assessment along subway based on PS InSAR technologies
CN111123300B (en) Near-real-time large-range high-precision ionosphere electron density three-dimensional monitoring method and device
CN110579787A (en) high-precision inclination monitoring method for electric power iron tower based on Beidou multi-antenna attitude measurement
CN104931022A (en) Satellite image three-dimensional area network adjustment method based on satellite-borne laser height measurement data
Davidson et al. Application of particle filters to a map-matching algorithm
WO2023134666A1 (en) Terminal positioning method and apparatus, and device and medium
CN105069295B (en) Satellite and surface precipitation measured value assimilation method based on Kalman filtering
CN109540095A (en) Roadbed settlement monitoring method based on satellite navigation and least square
CN103673995A (en) Calibration method of on-orbit optical distortion parameters of linear array push-broom camera
CN109100719B (en) Terrain map joint mapping method based on satellite-borne SAR (synthetic aperture radar) image and optical image
CN104575075A (en) City road network vehicle coordinate correcting method and device based on plough satellite
CN103217177A (en) Method, device and system for radio wave refraction correction
CN116182795B (en) Precision measurement method for vertical section of common speed railway
CN105044738A (en) Prediction method and prediction system for receiver autonomous integrity monitoring
Roberts et al. Structural dynamic and deflection monitoring using integrated GPS and triaxial accelerometers
CN109341665A (en) A kind of silt arrester fouling status investigating system and method
CN104567802A (en) Survey line land-sea elevation transfer method employing integrated shipborne gravity and GNSS
CN104280756A (en) Satellite positioning enhancing method based on receiver clock offset generalized prolongation approach method
CN117194928B (en) GNSS-based geographic deformation monitoring system
CN115638767B (en) Ground subsidence monitoring method
CN115457739B (en) Geological disaster early warning method and device, electronic equipment and storage medium
CN114088080A (en) Positioning device and method based on multi-sensor data fusion
CN102519469A (en) Planetary vehicle positioning method based on computer vision and VLBI combined adjustment

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