A quality-aware fixed-lag sliding window fusion method for tailings dam displacement monitoring

Abstract

This invention discloses a robust fixed-hysteresis sliding window fusion method for online displacement monitoring of tailings dams, combining RTK, dual millimeter-wave radar, and IMU quality sensing. The method constructs the sliding window state and multi-source factors through time synchronization, coordinate transformation, and asynchronous correlation. The covariance is dynamically adjusted based on sensor quality, and innovative gating, robust kernel, and progressive overriding are employed to suppress abnormal observations. Joint optimization is performed within a fixed window, and new priors are generated through edge-mapping, outputting three-dimensional displacement, quality scores, and anomaly indicators. This invention improves the continuity, stability, and interpretability of displacement monitoring in complex environments while controlling the computational load at the edge.

Description

Technical Field

[0001] This invention belongs to the technical fields of tailings dam safety monitoring, multi-source navigation and positioning and robust state estimation. Specifically, it relates to an online three-dimensional displacement monitoring method that unifies the modeling of RTK absolute position, dual millimeter-wave radar distance and IMU short-time motion constraints, and performs quality-aware robust optimization within a fixed hysteresis sliding window. Background Technology

[0002] Tailings dam displacement is characterized by small amplitude, slow change, long monitoring period, and high safety risk. RTK can provide absolute three-dimensional position, but due to obstruction caused by dam slope, mining equipment, and complex terrain, as well as multipath and communication fluctuations, its solution state may switch between fixed solution, floating-point solution, single-point solution, or invalid state. Directly using RTK coordinates can easily lead to abrupt jumps or interruptions. Millimeter-wave radar can provide non-contact range observation, but it is still affected by obstruction from cooperative targets, reflection changes, and geometric degradation. IMU has a high output frequency and can maintain short-term motion continuity, but it will accumulate drift during integration due to noise and zero bias.

[0003] Existing single-step filtering or simple weighting methods typically employ fixed observation covariance, making it difficult to simultaneously utilize solution state, number of satellites, accuracy factor, millimeter-wave signal quality, and residual statistics. Furthermore, when RTK recovers a fixed solution, immediately applying stronger constraints may cause a quadratic jump in the estimation results. On the other hand, the computational cost of global batch optimization continuously increases with data volume, making it unsuitable for long-term online operation at edge devices. Therefore, a fusion method with bounded computational scale, capable of utilizing recent historical states, identifying and continuously suppressing anomalous observations, and outputting quality indicators is needed. Summary of the Invention

[0004] (a) Purpose of the invention This invention aims to solve the problems of discontinuity, instability, and unexplainable anomalies in three-dimensional displacement sequences caused by time-varying RTK quality, millimeter-wave ranging obstruction or jumps, IMU integral drift, multi-source asynchronous operation, and limited computing resources at the edge, all of which contribute to complex tailings dam monitoring environments.

[0005] (II) Technical Solution To solve the above problems, the present invention adopts... Figure 1 The overall algorithm flow is shown below. RTK position and quality data, dual millimeter-wave distance data, and IMU inertial data first enter the time synchronization, coordinate transformation, and asynchronous association stages; then, the sliding window state establishment, IMU propagation or pre-integration, quality-aware covariance calculation, NIS gating, robust kernel weighting, and RTK asymptotic takeover are completed sequentially; multi-factor nonlinear joint optimization is performed within a fixed hysteresis window, and the oldest state is marginalized by Schur, finally outputting three-dimensional displacement, quality score, and anomaly flags.

[0006] (III) Beneficial Effects Compared with fixed covariance fusion, this invention enables RTK and millimeter-wave observations of different quality states to enter the optimization with different amounts of information, and continuously reduces the weight of moderate anomalies and rejects or weakens extreme anomalies; gradual takeover avoids the sudden dominance of estimation after RTK recovers a fixed solution; dual millimeter-wave distance forms planar redundancy with RTK and IMU, and the vertical gradual variation factor compensates for the structural deficiency of millimeter-wave not directly observing the Z-axis; fixed hysteresis sliding window reuses recent historical constraints, and Schur complement retains information of removed states, thus taking into account robustness, continuity and edge-end computational constraints. Attached Figure Description

[0007] Figure 1 This is a flowchart of the overall algorithm of the present invention; Figure 2 For fixed hysteresis sliding window factor plot; Figure 3 This is a flowchart for quality perception and anomaly suppression. Detailed Implementation

[0008] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings. Figure 1 Give the overall processing sequence from data input to displacement output; Figure 2 The connections between the prior factor, IMU factor, RTK factor, dual millimeter-wave distance factor and vertically slowly varying factor within a fixed hysteresis sliding window are given, as well as the process of marginalization of the oldest state and addition of new states; Figure 3 This paper presents a three-level NIS decision process, dynamic covariance, robust kernel, and gradual takeover procedure driven by RTK quality metrics, millimeter-wave quality metrics, and prediction residuals. The embodiments described are for illustrative purposes only and are not intended to limit the scope of protection of this invention.

[0009] S101: Multi-source data acquisition, time synchronization, and coordinate unification An RTK antenna, IMU, and two millimeter-wave radars are fixed to the node to be monitored, with the two millimeter-wave radars pointing towards the cooperative target whose location is known. A local monitoring coordinate system L is established based on the initial stable position or field calibration point. X The Y-axis is located on the primary monitoring plane, and the Z-axis is vertical. RTK data must contain at least three-dimensional position. z R,k Solution State s k Satellite number n sat,k HDOP or PDOP and optional reporting standard deviation; the j-th millimeter wave data includes distance. z M,j,k And signal strength, valid indicators or equivalent quality metrics; IMU data includes comparisons. f m angular velocity oh m And timestamps.

[0010] The original timestamps of each sensor are first corrected according to their calibration deviation relative to a unified clock, and then the observations are associated with the nearest state node, or short-time propagation is used between adjacent nodes; observations that exceed the allowable delay retain their quality marks and are rejected from entering the current window.

[0011] The RTK global coordinates are first subtracted from the reference position, and then converted to local coordinates through a rotation matrix from the global coordinate system to the local monitoring coordinate system.

[0012] S102: State Variables and Fixed-Hysteresis Sliding Window

[0013] x k =[ p k T, v k T, q k T, b a,k T, b g,k T ] T (1) In the formula, p k ∈ R 3 , v k ∈ R 3 Position and velocity, respectively. q k For unit attitude quaternions, b a,k ∈ R 3 , b g,k ∈ R 3 The accelerometer and gyroscope have zero bias, respectively. For low-dynamic implementations where the attitude remains essentially unchanged, this can be fixed. q k And omitted b g,k However, the residuals and dimension of the entire text should be simplified simultaneously.

[0014] X k ={ x k-N+1 , x k-N+2,…, x k (2) In the formula, N is the number of nodes in the sliding window state; the window duration is determined by N and the node interval, and is selected based on the computing power of the edge processor, the sensor frequency, and the expected rate of change of the tailings dam. Figure 2 Display the relationships between factors and state shifting in and out within the window.

[0015] S103: IMU State Propagation and Inertia Factor In the formula, n g,k and n a,k These are white noise from the gyroscope and accelerometer, respectively. Let the time interval between adjacent nodes be... Δt k , R ( q k ) is the rotation matrix corresponding to the attitude quaternion. g L Let be the gravity vector in the local coordinate system. After zero bias and white noise correction, the position, velocity, and attitude are propagated according to the following formula:

[0016] p k+1 - = p k + v k Δt k +1 / 2[ R ( q k ) f k + g L ] Δt k 2 ; v k+1 - = v k +[ R ( q k ) f k + g L ] Δt k ; q k+1 - = q k ⊗Exp(1 / 2 oh k Δt k (3) In the formula, the superscript "-" indicates propagation prediction, ⊗ indicates quaternion multiplication, and Exp(·) indicates rotation increment mapping. w ba,k and w bg,k It is a zero-biased random walk. The IMU residuals of adjacent states are constructed from the difference between the propagation prediction and the next state, and the information matrix is ​​obtained from the inverse of the process covariance.

[0017] The process covariance is obtained from the IMU noise density, zero-bias random walk, and node time interval propagation. When the sampling frequency is high, pre-integration of multiple IMU frames can also be performed between two state nodes, and the progressively propagated residual can be replaced with the pre-integrated residual.

[0018] S104: RTK Observation Model and Quality-Perceived Covariance RTK nominal covariance S R,0 According to quality standards α R,k Adjustments. To avoid confusion between the concepts of weight and covariance, we define a larger covariance scale as indicating less reliable observations, and the information matrix is ​​the inverse of the covariance:

[0019] S R,k = α R,k S R,0, Oh R,k = S R,k -1 ; α R,k = α s ( s k ) α n ( n sat,k ) α d (DOP k ) α σ ( s R,k (4) In the formula, α s To solve the state scale, satisfying α s (FIX)< α s (FLOAT)< αs (SINGLE) α s (INVALID); α n The number of satellites does not decrease. α d It does not decrease as HDOP or PDOP increases. α σ It does not decrease as the standard deviation of the RTK report increases. Invalid observations can cause... α R,k It tends to a predetermined large value or does not construct a factor directly.

[0020] The precision factor scale can be a squared ratio relative to the reference precision factor, and its range can be limited by upper and lower limit truncation functions; alternatively, a piecewise function or a lookup table with the same monotonicity can be used.

[0021] S105: Dual Millimeter Wave Range Model and Geometric Quality Let the position of the j-th cooperative target in the local coordinate system be... a j The lever arm of the radar ranging reference point relative to the monitoring node reference point is... l j If the change in posture is negligible, then R ( q k ) l j It can be replaced by a fixed lever arm.

[0022] h M,j ( x k )=‖ p k + R ( q k ) l j - a j ‖+ β j ; r M,j,k = z M,j,k - h M,j ( x k (5) In the formula, β j This represents the zero offset in ranging obtained through static calibration. Let... u j,kLet be the unit vector pointing from the cooperative target to the radar ranging reference point. The first-order Jacobian of the range residual with respect to position can be obtained from this unit direction vector.

[0023] Angle between two ranging directions i k Used to describe planar geometric conditions. When |sin i k When the value is small, the two factors are nearly collinear, the planar positioning information is weakened, and therefore the geometric covariance scale is increased accordingly.

[0024] S M,j,k = α sig,j,k· α jump,j,k· α geo,k· s M,j,0 2 (6) In the formula, ε is a positive number to prevent division by zero. α sig As a measure of signal quality. α jump This is the distance jump scale. s M.j,0 This represents the nominal ranging standard deviation. For distances that are out of lock, invalid, or exceed the physical change limit, the corresponding factor can be directly rejected.

[0025] S106: Innovation and NIS Level 3 Gating like Figure 3 As shown, the RTK solution state, number of satellites, and HDOP, as well as the millimeter-wave signal quality and ranging geometry, are read respectively, and the observation residuals are calculated in combination with the propagation state; an innovation vector is formed for each RTK factor or millimeter-wave factor. r i Jacobian matrix J i and innovation covariance S i Then, NIS statistics are used to perform a three-level judgment of normal, suspicious and abnormal.

[0026] c i = r i T S i -1 r i (7) In the formula, P k - To propagate state covariance, S i The observed covariance after quality perception. c i To normalize the innovation square. Based on the residual dimension. d i Given a selected false alarm probability, the first threshold is obtained from the chi-square distribution. t 1 and second threshold t 2, and t 2> t 1.

[0027] w NIS ( c )=1, c ≤ t 1; ( t 2- c ) / ( t 2- t 1), t 1< c ≤ t 2;0, c > t 2 (8) Normal observations retain their weights; suspicious observations are weighted using a continuous function to avoid output jitter caused by hard gating; abnormal observations are not factored, or are set as extremely weak constraints if the graph structure must be maintained. When both RTK and dual millimeter wave are effective at the same time, the corresponding planar displacement increments can be compared, and the covariance scale of suspicious observations can be increased accordingly.

[0028] S107: Robust Cores and Gradual RTK Takeover NIS provides statistical gating before factors enter the optimization process; the robust kernel continues to suppress outliers based on the updated residuals during iteration. Let the normalized residual norm be... s i =√( r i T Oh i r i The Huber kernel maintains a unit weight within the threshold and continuously reduces the weight after exceeding the threshold; the Cauchy kernel, on the other hand, decays smoothly as the residual increases.

[0029] In the formula, d H and d C These are the Huber threshold and the Cauchy scale, respectively. Preferably, the millimeter-wave distance uses the Huber kernel to preserve the continuous effect of moderate deviations, while the large jumps in the RTK plane use the Cauchy kernel to enhance the suppression of large residuals.

[0030] w i = w state,i · w quality,i · w NIS,i · w robust,i ; Oh i,eff = w i Oh i,0 (9) In the formula, w state To solve for the state weights, w quality The quality weights are formed by factors such as the number of satellites, DOP, signal quality, and geometric conditions. w NIS For gating weights, w robust For robust weights, Oh i,0 This is the nominal information matrix.

[0031] When RTK recovers from a FLOAT, SINGLE, or invalid state to a FIX state, let the number of consecutive stable fixed epochs after recovery be . c k The pipe length is L r :

[0032] or k =min(1, c k / L r ); Oh R,k,eff = or k w R,k Oh R,k (10) or k Gradually increase the value from near zero to 1 to achieve a smooth recovery of the RTK absolute constraint. Further adjustments can be made... or k Replace with a cubic smoothing function 3 u ²-2 u ³, of which u =clip(( c k ) / ( L r),0,1), so that the first-order rate of change of the take-up endpoint is zero.

[0033] S108: Fixed-hysteresis sliding window joint optimization like Figure 2 As shown, adjacent state nodes within the window are connected via IMU factors. Each valid state node can be connected to an RTK position factor and two millimeter-wave distance factors. The latest state is also connected to a vertically gradually varying factor. The prior factor formed by the marginalization result of the previous window constrains the first state of the window. When a new state enters, the oldest state is marginalized. Vertically gradually varying prediction can use a constant value model or a local constant velocity model.

[0034] When the difference between the RTK vertical observations and the prediction exceeds the threshold set by combining the window duration and the upper limit of short-term physical changes in the tailings dam, the amount of RTK vertical information is reduced, rather than forcibly assuming long-term linear settlement. The joint objective function is:

[0035] X k *=argmin X {‖ r p || Ωp 2 + S || r I || ΩI 2 + Wed R (‖ r R || ΩR,eff 2 )+ S jΣρM (‖ r M,j || ΩM,j,eff 2 )+ S || r Z || ΩZ 2 (11) In the formula, r p The prior residuals generated by the marginalization of the previous window. r R and r M For the corresponding robust loss, Oh Z This is a vertically varying information matrix. All summations only cover factors relevant to valid observations within the current window.

[0036] Nonlinear least squares solution At the t-th iteration, each residual is linearized to the first order.r i (X⊞δx)≈ r i (X)+ J iδx By accumulating these, we obtain the Hessian approximation matrix H and the gradient g, and then solve the following equation:

[0037] ( H + λI ) δx =- g (12) In the formula, l is the non-negative damping coefficient, I is the identity matrix, and ⊞ indicates that ordinary addition is used for position, velocity, and zero bias, and rotating manifold increment is used for attitude. The process stops when the increment norm is less than the convergence threshold, the change in the objective function is less than the threshold, or the maximum number of iterations is reached.

[0038] S109: Schur marginalization When a new state enters and the number of window nodes exceeds N, the oldest state and only the factors connected to it are linearized. The normal equations are divided into blocks according to the variable to be marginalized δxm and the retained variable δxr, and the variable to be marginalized is eliminated by Schur elimination.

[0039] H p = H rr - H rm H mm -1 H mr ; b p = b r - H rm H mm -1 b m (13) In the formula, H mm The sub-block corresponding to the variable to be marginalized H rr To preserve the sub-blocks corresponding to the variables, H mr and H rm These are cross-sub-blocks. To improve numerical stability, Cholesky decomposition or solving linear equations can be used instead of explicit inversion. H p and b pThis forms the prior knowledge of the retained variables in the next sliding window, allowing the window length to be fixed while retaining historical information.

[0040] S110: Displacement, Mass, and Abnormal Output

[0041] d k = p k *- p 0; d X = k = e X T d k, d Y = k = e Y T d k, d Z = k = e Z T d k ; d 3D,k =√( d X,k 2 + d Y,k 2 + d Z,k 2 (14) In the formula, p k * indicates the latest optimized location. p 0 represents the initial stable position or the field reference position. e X , e Y and e Z This is a unit vector in a local coordinate system. The output also includes the receiving status of each sensor, NIS value, effective weight, RTK takeover coefficient, solution status, and overall quality score. The overall quality score can be obtained by a normalized combination of the effective factor ratio, Hessian condition index, and the quality weights of each sensor; its specific mapping does not limit the scope of protection of this invention.

Claims

1. A quality-sensing fixed-hysteresis sliding window fusion method for tailings dam displacement monitoring, characterized in that, include: S101: Acquire RTK 3D position and quality information, distance and quality information from two millimeter-wave radars, and IMU inertial data; perform unified time stamping, coordinate transformation, asynchronous association, and caching on the multi-source data. S102: Establish a state including the monitoring point's position, velocity, attitude, accelerometer zero bias, and gyroscope zero bias, and form a fixed-hysteresis sliding window from the N most recent states. S103: Establish the inertial propagation factor between adjacent states based on IMU data. S104: Establish the quality-sensing RTK position factor based on the RTK solution state, number of satellites, accuracy factor, reporting standard deviation, and consecutive fixed epochs. S105: Based on the two millimeter-wave radars and their corresponding coordinates... S106: Establish a dual-distance factor based on the geometric relationship between the targets; S107: Perform three-level gating (normal, suspicious, and abnormal) on RTK observations and millimeter-wave observations based on innovation vector, innovation covariance, and normalized innovation square; S108: Combine the quality weight with the robust kernel weight, and apply an asymptotic takeover coefficient when the RTK recovers from the degenerate solution to the fixed solution; S109: Perform nonlinear least squares joint optimization on the prior factor, inertial propagation factor, RTK position factor, dual-distance factor, and vertically slowly varying factor in the fixed hysteresis sliding window; S100: Perform Schur-complement marginalization on the oldest state in the sliding window and use the marginalization result as the prior for the next sliding window. S110 outputs the three-dimensional displacement, quality score, and anomaly markers of the tailings dam monitoring points based on the difference between the latest optimized position and the benchmark position.

2. The method according to claim 1, characterized in that... In step S102 k The states are x k =[ p k T , v k T , q k T , b a,k T , b g,k T ] T ,in p k , v k and q k These are the position, velocity, and attitude quaternions in the local monitoring coordinate system, respectively. b a,k and b g,k These are zero biases for the accelerometer and gyroscope, respectively; the fixed hysteresis sliding window is... X k = x (k-N+1) ,…, x k N is a positive integer. Step S103 includes: removing the corresponding zero bias from the IMU angular velocity and specific force; propagating the position, velocity and attitude based on the attitude rotation matrix and gravity vector; modeling the accelerometer zero bias and gyroscope zero bias as random walks; constructing the IMU residual based on the difference between the propagation state and the adjacent states within the window, and determining the information matrix of the IMU residual from the process noise.

3. The method according to claim 1, characterized in that, Step S104 multiplies the nominal RTK covariance by a quality scale factor, which consists of at least two of the following: solution state scale, satellite number scale, accuracy factor scale, reporting standard deviation scale, and continuous fixed epoch scale. The scale corresponding to a fixed solution is smaller than the scales corresponding to floating-point solutions and single-point solutions. The RTK position factor corresponding to an invalid solution is rejected or set as a weak constraint. In step S105, the predicted value of the j-th millimeter-wave range is the Euclidean distance between the j-th radar ranging reference point and the j-th cooperative target on the mobile monitoring node. The range residual is constructed based on the difference between the predicted value and the measured distance, and the covariance or information matrix of the dual range factors is adjusted according to the millimeter-wave signal quality, range jump variable, and the angle between the two ranging directions.

4. The method according to claim 1, characterized in that, Step S106: Calculate the normalized innovation square. γ = r T S (-1) r When γ is not greater than the first threshold, the observation weight is maintained; when γ is greater than the first threshold but not greater than the second threshold, the observation weight is continuously reduced; when γ is greater than the second threshold, the observation is rejected or set as a predetermined weak constraint. The first and second thresholds are determined by the residual dimension and the confidence level of the chi-square distribution.

5. The method according to claim 1, characterized in that, Step S107 applies a Huber kernel or Cauchy kernel to the normalized residuals to obtain robust weights, and multiplies the robust weights, solution state weights, quality index weights, NIS gating weights, and geometric quality weights to scale the nominal information matrix corresponding to the observation factors. When the RTK recovers from a floating-point solution, a single-point solution, or an invalid state to a fixed solution, a continuous stable fixed epoch count is set. c k and pipe length L r and make the gradual takeover coefficient η k =min1, c k / L r Alternatively, a smoothing function with the same monotonic and continuous characteristics can be used to gradually increase the effective information content of the RTK position factor as the number of continuous, stable, and fixed epochs increases.

6. The method according to claim 1, characterized in that, In step S108, the vertical gradual variation factor is predicted as a constant or locally linear value based on the accepted vertical state within the window. When the jump variable of the RTK vertical observation relative to the predicted value exceeds the physical change threshold corresponding to the window duration, the RTK vertical information content is reduced, and short-term continuity is maintained by the vertical gradual variation factor. Step S108 is solved using the Gauss-Newton method, the damped Gauss-Newton method, or the Levenberg-Marquardt method; the linearized normal equation is ( H + λI ) δx =- g And update the increment δx to the window state using manifold addition. Step S109 divides the linearization equation into blocks according to the variables to be marginalized and the variables to be retained, and uses... H p = H rr - H rm H mm (-1) H m r and b p = b r - H rm H mm (-1) b m Generate priors about the retained variables, where H mm This refers to the Hessian sub-block corresponding to the variable to be marginalized.

7. An online displacement monitoring device for tailings dams, comprising a processor, a memory, an RTK receiving unit, two millimeter-wave radars and an IMU, wherein the memory stores a program, and when the program is executed by the processor, it implements the method described in any one of claims 1 to 6.

8. A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method of any one of claims 1 to 6.

Images