CN106371092A

CN106371092A - Deformation monitoring method based on GPS and strong-motion seismograph observation adaptive combination

Info

Publication number: CN106371092A
Application number: CN201610726219.9A
Authority: CN
Inventors: 涂锐; 卢晓春; 成芳; 张鹏飞; 张兴刚
Original assignee: National Time Service Center of CAS
Current assignee: National Time Service Center of CAS
Priority date: 2016-08-25
Filing date: 2016-08-25
Publication date: 2017-02-01
Anticipated expiration: 2036-08-25
Also published as: CN106371092B

Abstract

The invention provides a deformation monitoring method based on the self-adaptive combination of GPS and strong motion instrument observation, which integrates the high-resolution strong motion instrument acceleration observation into the GPS-PPP positioning model, and takes the baseline drift error of the strong motion instrument as The unknown parameters are estimated together with other positioning parameters in real time; at the same time, the optimization of parameter estimation is realized by adaptively adjusting the dynamic noise and the solution time window, which not only efficiently solves the correction problem of baseline drift, but also enhances the solution strength of GPS. The high-frequency noise is reduced, the complementary advantages of the two observation techniques are realized, and high-precision wide-band deformation information can be obtained in real time, which has broad application prospects. The invention has the characteristics of advanced algorithm, high precision, strong stability, high degree of automation, etc., and has the overall packaging of multi-satellite systems such as BDS, GPS, GLONASS, etc., is simple to implement, and can be applied to various application occasions.

Description

A Deformation Monitoring Method Based on Adaptive Combination of GPS and Strong Motion Instrument Observation

技术领域technical field

本发明涉及一种地表形变监测方法。The invention relates to a method for monitoring ground deformation.

背景技术Background technique

GPS和强震仪观测是获取高精度地表形变(位移、速度、加速度)的两种有效手段，它们已广泛应用于自然灾害监测且各具特色。GPS易于获取高精度位移，但存在采样频率低、高频信噪比低、信号稳定性差的缺陷；同时，强震仪易于获取高分辨率加速度，但因基线漂移误差的存在，其积分后的速度和位移常存在偏差。当前的数据处理方式大多是单传感器模式，导致多传感器的观测资源没有充分利用。如何将两类观测数据进行有机组合，实现优势互补，实时提供高精度宽频带的形变信息，更好的服务于灾害等形变监测具有重要价值。GPS and strong motion instrument observations are two effective means to obtain high-precision surface deformation (displacement, velocity, acceleration). They have been widely used in natural disaster monitoring and have their own characteristics. GPS is easy to obtain high-precision displacement, but it has the defects of low sampling frequency, low high-frequency signal-to-noise ratio, and poor signal stability. Displacement is often biased. Most of the current data processing methods are single-sensor mode, which leads to the underutilization of multi-sensor observation resources. How to organically combine the two types of observation data to achieve complementary advantages, provide high-precision broadband deformation information in real time, and better serve deformation monitoring such as disasters is of great value.

发明内容Contents of the invention

为了克服现有技术的不足，本发明提供一种基于GPS与强震仪观测自适应组合的形变监测方法，将高分辨率的强震仪加速度观测融入到GPS-PPP定位模型中，并将强震仪的基线漂移误差当作未知参数与其它定位参数一起进行实时估计，同时通过自适应调整动态噪声和求解时间窗口来实现参数估计的最优化，既高效解决了基线漂移的校正问题，同时又增强了GPS的求解强度，降低了高频噪声，实现了两种观测技术的优势互补，可以实时获取高精度宽频带的形变信息。In order to overcome the deficiencies of the prior art, the present invention provides a deformation monitoring method based on the adaptive combination of GPS and strong motion instrument observations, which integrates the high-resolution strong motion instrument acceleration observations into the GPS-PPP positioning model, and combines the strong motion The baseline drift error of the seismograph is used as an unknown parameter for real-time estimation together with other positioning parameters. At the same time, the optimization of parameter estimation is realized by adaptively adjusting the dynamic noise and solving the time window, which not only efficiently solves the correction problem of baseline drift, but also The solution strength of GPS is enhanced, high-frequency noise is reduced, and the complementary advantages of the two observation technologies are realized, and high-precision wide-band deformation information can be obtained in real time.

本发明解决其技术问题所采用的技术方案是：通过两种观测数据的融合解算得到的基本的结果信息后，自适应的调整求解策略，实现参数估计最优化和两种观测手段的优势互补，具体包括以下步骤：The technical solution adopted by the present invention to solve the technical problem is: After the basic result information obtained through the fusion of two observation data, the solution strategy is adaptively adjusted to realize the optimization of parameter estimation and the complementary advantages of the two observation methods , including the following steps:

第一步，获得GPS的相位/伪距观测值、强震仪的加速度观测值、以及GPS数据处理需要的辅助产品，包括GPS卫星轨道、卫星钟差和地球自转参数；The first step is to obtain the phase/pseudorange observation value of GPS, the acceleration observation value of strong motion instrument, and the auxiliary products required for GPS data processing, including GPS satellite orbit, satellite clock error and earth rotation parameters;

第二步，对第一步获得的GPS的相位/伪距观测值、强震仪的加速度观测值、以及GPS数据处理需要的辅助产品进行数据检查、粗差剔除和周跳探测，得到干净的数据；对所述干净的数据进行相对论、潮汐、天线相位中心、对流层和地球自转误差的修正；In the second step, the phase/pseudo-range observations of GPS obtained in the first step, the acceleration observations of the strong motion instrument, and the auxiliary products required for GPS data processing are carried out for data inspection, gross error elimination and cycle slip detection to obtain a clean data; corrections for relativistic, tidal, antenna phase center, tropospheric and earth rotation errors are applied to said clean data;

第三步，建立GPS与强震仪观测的紧组合模型，所述GPS和强震仪观测的紧组合模型如下：The third step is to establish a tight combination model of GPS and strong motion instrument observations. The tight combination model of GPS and strong motion instrument observations is as follows:

${[\begin{matrix} {L L}_{c c} \\ {P P}_{c c} \end{matrix}]}_{k k} = = [\begin{matrix} {e e}_{k k} & {m m}_{z z,, k k} & 11 & 11 \\ {e e}_{k k} & {m m}_{z z,, k k} & 11 & 00 \end{matrix}] [\begin{matrix} {x x}_{k k} \\ {z z}_{k k} \\ {δt δt}_{k k} \\ b b \end{matrix}] + + [\begin{matrix} {ϵ ϵ}_{{L L}_{c c}} \\ {ϵ ϵ}_{{P P}_{c c}} \end{matrix}],, \begin{matrix} {ϵ ϵ}_{{L L}_{c c}} ~ ~ N N ((00,, {σ σ}_{{L L}_{c c}}^{22})) \\ {ϵ ϵ}_{{P P}_{c c}} ~ ~ N N ((00,, {σ σ}_{{P P}_{c c}}^{22})) \end{matrix}$

$\begin{matrix} {a a}_{k k} = = {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} - - {u u}_{k k} + + {ϵ ϵ}_{{a a}_{k k}} & {ϵ ϵ}_{{a a}_{k k}} ~ ~ N N ((00,, {σ σ}_{{a a}_{k k}}^{22})) \end{matrix}$

其中，L_c和P_c是扣除卫星和接收机间几何距离后的载波相位残差和伪距码观测值残差；e是卫星与接收机之间的单位矢量；角标k表示历元数；m_z,k是对流层湿延迟天顶方向投影函数；x、分别是接收机的位移和加速度矢量；z、δt和b分别为对流层天顶延迟、接收机钟差和相位模糊度；ε是测量噪声，ε_Lc和ε_Pc分别表示载波和伪距观测的噪声，其对应的方差为δ²；a和u分别为强震仪加速度观测值和基线漂移误差；基线漂移误差参数在设定宽度的每个估计窗口内作为常数估计，所述的天顶对流层延迟作为常数进行估计或者表现为随机游走过程，所述的接收机钟差作为高斯白噪声逐历元估计，所述的相位模糊度在连续无周跳情况下作为常数估计；Among them, L _c and P _c are the carrier phase residuals and pseudorange code observation residuals after deducting the geometric distance between the satellite and the receiver; e is the unit vector between the satellite and the receiver; the subscript k represents the number of epochs ; m _z,k is the projection function of the tropospheric wet delay zenith direction; x, are the displacement and acceleration vectors of the receiver, respectively; z, δt, and b are the tropospheric zenith delay, receiver clock error, and phase ambiguity, respectively; ε is the measurement noise, and ε _Lc and ε _Pc represent the noise of carrier and pseudorange observations, respectively , and its corresponding variance is δ ² ; a and u are the acceleration observation value of the strong motion instrument and the baseline drift error respectively; the baseline drift error parameter is estimated as a constant in each estimation window of the set width, and the zenith tropospheric delay Estimated as a constant or manifested as a random walk process, the receiver clock error is estimated as Gaussian white noise epoch by epoch, and the phase ambiguity is estimated as a constant under continuous cycle-slip-free conditions;

接收机的运动和基线漂移误差的状态方程如下：The state equations of the receiver motion and baseline drift error are as follows:

$[\begin{matrix} {x x}_{k k} \\ {\overset{\cdot &Center Dot;}{x x}}_{k k} \\ {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} \\ {u u}_{k k} \end{matrix}] = = [\begin{matrix} I I & 0.5 0.5 I I τ τ & 0.5 0.5 {Iτ Iτ}^{22} & 00 \\ 00 & I I & I I τ τ & 00 \\ 00 & 00 & I I & 00 \\ 00 & 00 & 00 & I I \end{matrix}] [\begin{matrix} {x x}_{k k - - 11} \\ {\overset{\cdot \cdot}{x x}}_{k k - - 11} \\ {a a}_{k k} - - {\overset{^^}{u u}}_{k k - - 11} \\ {u u}_{k k - - 11} \end{matrix}] + + {w w}_{k k}$

其中，为接收机的速度矢量，I是3x3的单位矩阵，τ是GPS和强震仪的最小采样率，w_k是动态噪声，其数学期望E(w_k)＝0，方差q_u为基线漂移的功率密度，k-1表示上一历元；采用基线漂移改正后加速度来代替加速度预测值a_k-1；给定观测值权比和动态噪声；采用卡尔曼滤波进行参数估计，同时进行残差编辑和精度统计，得到接收机的位置、速度、加速度、基线漂移、对流层延迟和接收机钟差信息及精度指标；in, is the velocity vector of the receiver, I is the 3x3 identity matrix, τ is the minimum sampling rate of GPS and strong motion instrument, w _k is the dynamic noise, its mathematical expectation E(w _k )=0, variance q _u is the power density of baseline drift, k-1 means the previous epoch; the acceleration after baseline drift correction to replace the acceleration prediction value a _k-1 ; given the weight ratio of observations and dynamic noise; use Kalman filter to estimate parameters, and at the same time carry out residual editing and precision statistics to obtain the receiver's position, velocity, acceleration, baseline drift, Tropospheric delay and receiver clock information and accuracy indicators;

第四步，根据第三步的求解信息计算基线漂移的时间序列，并构建四个学习统计量，确定统计量的标准；四个统计量分别为：The fourth step is to calculate the time series of baseline drift based on the solution information in the third step, and construct four learning statistics to determine the standard of the statistics; the four statistics are:

基线漂移的位移影响 Displacement Effects of Baseline Drift

原始加速度的位移影响 Displacement Effect of Raw Acceleration

改正后的加速度的位移影响 Corrected acceleration for displacement effects

信号强度 signal strength

其中，t₁、t₂是定义的估计窗口的起始时间和终止时间，分别为起始时间和终止时间对应的历元数，是a的均值，每个学习统计量的标准A₀，B₀，C₀，D₀是静态观测期间每个学习统计量在设定时间序列中的三倍标准差；Among them, t ₁ and t ₂ are the start time and end time of the defined estimation window, are the epoch numbers corresponding to the start time and end time, respectively, is the mean value of a, the standard of each learning statistic A ₀ , B ₀ , C ₀ , D ₀ is three times the standard deviation of each learning statistic in the set time series during static observation;

第五步，根据第四步确定的标准将基线漂移误差参数分成四个阶段，并计算新的动态噪声和时间估计窗口；四个阶段的确定如下：In the fifth step, the baseline drift error parameters are divided into four stages according to the standard determined in the fourth step, and a new dynamic noise and time estimation window are calculated; the determination of the four stages is as follows:

初始化阶段：满足(A₂&B₁)或(B₁&C₂)，Initialization phase: satisfy (A ₂ &B ₁ ) or (B ₁ &C ₂ ),

静态阶段：满足(A₁)，Static phase: satisfy (A ₁ ),

瞬时基线漂移阶段：满足(A₂&D₂)，Instantaneous baseline drift phase: satisfy (A ₂ &D ₂ ),

永久基线漂移阶段：满足(A₂&D₁)&(B₂&C₁)；Permanent baseline drift phase: satisfy (A ₂ &D ₁ )&(B ₂ &C ₁ );

在瞬时基线漂移阶段，动态噪声调整如下：During the transient baseline drift phase, the dynamic noise is adjusted as follows:

${q q}_{u u}^{' '} = = \sqrt{\frac{| | {u u}_{{t t}_{22}} - - {u u}_{{t t}_{11}} | |}{(({t t}_{22} - - {t t}_{11}))}}$

其中，q′_u为基线漂移的新的功率密度；基线漂移估计窗口的长度改变为每秒一次；对于其它三个阶段，动态噪声q_u＝0.0001-0.001m^1/2/s^3/2，窗口长度保持为5秒；Among them, q′ _u is the new power density of the baseline drift; the length of the baseline drift estimation window is changed once per second; for the other three stages, the dynamic noise q _u =0.0001-0.001m ^1/2 /s ^3/2 , The window length is kept at 5 seconds;

第六步，根据第五步中确定的新的动态噪声和新的时间估计窗口，采用卡尔曼滤波进行自适应组合解算，获得最优组合结果。In the sixth step, according to the new dynamic noise and the new time estimation window determined in the fifth step, Kalman filter is used for adaptive combination solution to obtain the optimal combination result.

本发明的有益效果是：The beneficial effects of the present invention are:

第一，实现两种观测手段的优势互补，提供高精度宽频带的形变信息。First, realize the complementary advantages of the two observation methods, and provide high-precision and wide-band deformation information.

本发明将高精度GPS和高分辨率强震仪观测进行紧组合处理，弥补了GPS技术信噪比差和强震仪技术基线漂移的局限性，实现了他们的优势互补，其组合结果可以实时输出高精度宽频带的位移、速度、加速度信息，可直接应用于灾害等形变监测。The invention combines the observations of high-precision GPS and high-resolution strong motion instruments tightly, which makes up for the limitations of the poor signal-to-noise ratio of GPS technology and the baseline drift of strong motion instrument technology, realizes their complementary advantages, and the combined results can be real-time Output high-precision broadband displacement, velocity, and acceleration information, which can be directly applied to deformation monitoring such as disasters.

第二，有效识别并校正了强震仪的基线漂移，为基线漂移的校正提供了新途径。Second, the baseline drift of the strong motion instrument is effectively identified and corrected, which provides a new way for the correction of the baseline drift.

强震仪的基线漂移误差的识别是校正是地震领域棘手的问题，通常采用经验的方法进行完成，但经验的方法不能准确描述基线漂移误差的特性，校正后的结果存在不同程度的系统偏差，且不能实时进行。本发明将强震仪的基线漂移误差进行模型化参数估计，既可以有效的识别并校正强震仪的基线漂移，同时为基线漂移的校正提供了新途径，对旋转地震学的研究具有重要意义。The identification and correction of the baseline drift error of the strong motion instrument is a thorny problem in the field of seismicity, which is usually completed by empirical methods, but the empirical method cannot accurately describe the characteristics of the baseline drift error, and there are different degrees of systematic deviation in the corrected results. and cannot be performed in real time. The present invention carries out modeling parameter estimation on the baseline drift error of the strong motion instrument, which can effectively identify and correct the baseline drift of the strong motion instrument, and at the same time provides a new way for the correction of the baseline drift, which is of great significance to the research of rotational seismology .

第三，有效的降低了GPS的观测噪声，加快了GPS定位的收敛速度。Third, it effectively reduces the observation noise of GPS and accelerates the convergence speed of GPS positioning.

高分辨率的加速度观测融合到GPS定位方程中，可以约束GPS求解强度，加快GPS定位的收敛速度，同时可以有效的降低GPS的噪声，提高GPS定位精度和稳定性。The high-resolution acceleration observation is integrated into the GPS positioning equation, which can constrain the strength of the GPS solution and speed up the convergence speed of the GPS positioning. At the same time, it can effectively reduce the noise of the GPS and improve the accuracy and stability of the GPS positioning.

附图说明Description of drawings

图1是紧组合处理中的四个阶段示意图；Figure 1 is a schematic diagram of the four stages in the compact combination process;

图2是GPS与强震仪观测组合框图；Figure 2 is a combination block diagram of GPS and strong motion instrument observation;

图3是GPS与强震仪自适应组合框图。Figure 3 is a block diagram of the adaptive combination of GPS and strong motion instrument.

具体实施方式detailed description

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。The present invention will be further described below in conjunction with the accompanying drawings and embodiments, and the present invention includes but not limited to the following embodiments.

本发明将高分辨率的强震仪加速度观测融入到GPS-PPP定位模型中，并将强震仪的基线漂移误差当作未知参数与其它定位参数一起进行实时估计，同时通过自适应调整动态噪声和求解时间窗口来实现参数估计的最优化，既高效解决了基线漂移的校正问题，同时又增强了GPS的求解强度，降低了高频噪声，实现了两种观测技术的优势互补，可以实时获取高精度宽频带的形变信息。The invention integrates the high-resolution strong motion instrument acceleration observation into the GPS-PPP positioning model, and uses the baseline drift error of the strong motion instrument as an unknown parameter for real-time estimation together with other positioning parameters, and at the same time adjusts the dynamic noise through self-adaption and solve the time window to achieve the optimization of parameter estimation, which not only efficiently solves the problem of baseline drift correction, but also enhances the solution strength of GPS, reduces high-frequency noise, realizes the complementary advantages of the two observation techniques, and can obtain real-time High-precision wide-band deformation information.

本发明的技术方案主要包括两大核心算法：The technical solution of the present invention mainly includes two core algorithms:

(1)GPS与强震仪观测的紧组合处理(1) Tight combination processing of GPS and strong motion instrument observations

GPS和强震仪观测的紧组合采用下列观测方程：The tight combination of GPS and strong motion meter observations uses the following observation equations:

${[\begin{matrix} {L L}_{c c} \\ {P P}_{c c} \end{matrix}]}_{k k} = = [\begin{matrix} {e e}_{k k} & {m m}_{z z,, k k} & 11 & 11 \\ {e e}_{k k} & {m m}_{z z,, k k} & 11 & 00 \end{matrix}] [\begin{matrix} {x x}_{k k} \\ {z z}_{k k} \\ {δt δt}_{k k} \\ b b \end{matrix}] + + [\begin{matrix} {ϵ ϵ}_{{L L}_{c c}} \\ {ϵ ϵ}_{{P P}_{c c}} \end{matrix}],, \begin{matrix} {ϵ ϵ}_{{L L}_{c c}} ~ ~ N N ((00,, {σ σ}_{{L L}_{c c}}^{22})) \\ {ϵ ϵ}_{{P P}_{c c}} ~ ~ N N ((00,, {σ σ}_{{P P}_{c c}}^{22})) \end{matrix} - - - - - - ((11))$

$\begin{matrix} {a a}_{k k} = = {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} - - {u u}_{k k} + + {ϵ ϵ}_{{a a}_{k k}} & {ϵ ϵ}_{{a a}_{k k}} ~ ~ N N ((00,, {σ σ}_{{a a}_{k k}}^{22})) \end{matrix} - - - - - - ((22))$

其中，L_c和P_c是扣除卫星和接收机间几何距离后的载波相位和伪距码观测值残差，e是卫星与接收机之间的单位矢量，k表示历元。m是对流程湿延迟投影函数，x，分别是接收机位移、加速度矢量，z，δt和b分别为对流层天顶延迟、接收机钟差和相位模糊度，ε是测量噪声，其对应的方差为δ²。a和u分别为强震仪加速度观测值和基线漂移误差，k表示历元数。基线漂移参数在每个估计窗口中作为常数估计，其它的参数，如天顶对流层延迟，在短时间内(一个小时或两个小时)可以作为常数进行估计，或者表现为随机游走过程，而接收机钟差当作高斯白噪声逐历元估计，模糊度参数在连续无周跳情况下作为常数估计。Among them, L _c and P _c are the residuals of carrier phase and pseudorange code observations after deducting the geometric distance between the satellite and the receiver, e is the unit vector between the satellite and the receiver, and k represents the epoch. m is the wet delay projection function for the process, x, are receiver displacement and acceleration vector, z, δt and b are tropospheric zenith delay, receiver clock error and phase ambiguity respectively, ε is measurement noise, and its corresponding variance is δ ² . a and u are the acceleration observation value and baseline drift error of the strong motion instrument, respectively, and k represents the number of epochs. The baseline drift parameter is estimated as a constant in each estimation window, other parameters, such as the zenith tropospheric delay, can be estimated as a constant in a short period of time (one hour or two hours), or behave as a random walk process, while The receiver clock error is estimated as Gaussian white noise epoch by epoch, and the ambiguity parameters are estimated as constants in the continuous cycle-slip condition.

测站运动和基线漂移的状态方程表示如下：The state equations of station motion and baseline drift are expressed as follows:

$[\begin{matrix} {x x}_{k k} \\ {\overset{\cdot \cdot}{x x}}_{k k} \\ {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} \\ {u u}_{k k} \end{matrix}] = = [\begin{matrix} I I & 0.5 0.5 I I τ τ & 0.5 0.5 {Iτ Iτ}^{22} & 00 \\ 00 & I I & I I τ τ & 00 \\ 00 & 00 & I I & 00 \\ 00 & 00 & 00 & I I \end{matrix}] [\begin{matrix} {x x}_{k k - - 11} \\ {\overset{\cdot \cdot}{x x}}_{k k - - 11} \\ {a a}_{k k} - - {\overset{^^}{u u}}_{k k - - 11} \\ {u u}_{k k - - 11} \end{matrix}] + + {w w}_{k k},, - - - - - - ((33))$

其中为速度矢量，I是3x3的单位矩阵，τ是GPS和强震仪的最小采样率，w_k是动态噪声，其数学期望E(w_k)＝0，方差q_u为基线漂移的功率密度。k-1表示上一历元，这里采用基线漂移改正后加速度来代替加速度预测值a_k-1，避免加速度和基线漂移具有大的动态噪声。应用公式(1)，(2)及(3)，给定观测值权比和动态噪声(这里经验权重由测量噪声来确定，通常的取值范围分别为：动态噪声给定为常数，q_u＝0.003～0.01m^1/2/s^3/2，窗口时间长度为1-10秒)，可用卡尔曼滤波进行参数估计既可以得到位置，速度，加速度、基线漂移、对流层延迟和接收机钟差。另外，上述随机模型可自适应地调整使得参数估计最优化。in is the velocity vector, I is the identity matrix of 3x3, τ is the minimum sampling rate of GPS and strong motion instrument, w _k is the dynamic noise, its mathematical expectation E(w _k )=0, variance q _u is the power density of baseline drift. k-1 represents the previous epoch, where the acceleration after baseline drift correction is used to replace the acceleration prediction value a _k-1 , to avoid acceleration and baseline drift with large dynamic noise. Apply formulas (1), (2) and (3), given the weight ratio of observations and dynamic noise (here the empirical weight is determined by the measurement noise, and the usual value ranges are: The dynamic noise is given as a constant, q _u = 0.003～0.01m ^1/2 /s ^3/2 , the window time length is 1-10 seconds), and Kalman filter can be used for parameter estimation to obtain position, velocity, acceleration, baseline Drift, tropospheric delay and receiver clock bias. In addition, the above stochastic model can be adaptively adjusted to optimize parameter estimation.

(2)GPS与强震仪观测自适应组合处理(2) Adaptive combined processing of GPS and strong motion instrument observations

对于GPS和强震仪观测的紧组合的稳健估计处理，需要解决两个关键问题。第一个是基线漂移是否发生？如果发生，则需要估算基线漂移的持续时间。第二个是基线漂移的变化情况，于是本发明可以确定状态向量的最佳动态噪声。For robust estimation processing of tight combinations of GPS and strong motion meter observations, two key issues need to be addressed. The first is does baseline drift occur? If it occurs, the duration of the baseline drift needs to be estimated. The second is the variation of the baseline drift, so the present invention can determine the optimal dynamic noise of the state vector.

为了解决这些问题，常采用自适应滤波处理。首先定义一些学习的统计量来检测基线漂移是否发生，并确定它是否显著，然后调整动态噪声来实现参数求解的自适应估计。在卫星导航领域中，这些统计量通常是基于状态差异、方差分量和预测残差来定义。但它们不能很好的应用于GPS和强震紧组合处理。一方面，对于GPS和强震仪紧组合处理，引入了基线漂移和相位模糊参数，单历元的有限观测量估计的几何状态参数偏差较大甚至不能估计。在这种情况下，状态差异或方差分量的统计特性会很差。另一方面，强震记录的基线漂移对积分后速度位移非常敏感。虽然基线漂移的值相比峰值加速度的值通常是非常小的，但它们将显著影响积分后的速度和位移。此外，基线漂移误差主要是由地面倾斜或旋转造成，并有一个重要特征，在没有运动时为零或恒定值，因此，选择基线漂移作为学习统计量更为合适。再者，高分辨率的加速度信号也可以用于检测运动状态。这里先假设静态状态下基线漂移和加速度为零，然后定义了四个学习统计量，来识别是否存在基线漂移以及基线漂移变化情况。由于强震仪采样率比GPS高得多，选择GPS采样间隔持续的时间作为基线漂移估计窗口的长度。同时，为避免随机误差的影响，在整个基线漂移估计窗口将基线漂移和加速度积分成位移进行分析。In order to solve these problems, adaptive filtering is often used. We first define some learned statistics to detect whether baseline drift occurs and determine whether it is significant, and then adjust the dynamic noise to achieve adaptive estimation of the parameter solution. In the field of satellite navigation, these statistics are usually defined based on state differences, variance components and prediction residuals. But they can't be well applied to the combined processing of GPS and strong earthquakes. On the one hand, for the tight combination processing of GPS and strong motion instruments, baseline drift and phase ambiguity parameters are introduced, and the geometric state parameters estimated by the limited observations of a single epoch have large deviations or even cannot be estimated. In this case, the statistical properties of the state differences or variance components would be poor. On the other hand, the baseline drift of strong motion records is very sensitive to the integrated velocity displacement. Although the value of baseline drift is usually very small compared to the value of peak acceleration, they will significantly affect the integrated velocity and displacement. In addition, the baseline drift error is mainly caused by ground tilt or rotation, and has an important characteristic, which is zero or a constant value when there is no motion, so it is more appropriate to choose baseline drift as a learning statistic. Furthermore, high-resolution acceleration signals can also be used to detect motion states. Here, it is assumed that the baseline drift and acceleration are zero in the static state, and then four learning statistics are defined to identify whether there is baseline drift and the change of baseline drift. Since the sampling rate of the strong motion instrument is much higher than that of GPS, the duration of the GPS sampling interval is selected as the length of the baseline drift estimation window. At the same time, in order to avoid the influence of random errors, the baseline drift and acceleration are integrated into displacement for analysis in the entire baseline drift estimation window.

统计量1: Statistics 1:

统计量2: Statistics 2:

统计量3: Statistics 3:

统计量4: Statistics 4:

其中t₁，t₂是定义的估计时间窗口(一般定义为采样间隔的倍数)，为历元数，是a的均值，A代表基线漂移的位移影响，B代表原始加速度的位移影响，C代表改正后的加速度的位移影响，D代表信号强度。每个学习统计量根据定义的标准A₀，B₀，C₀，D₀都有两种情况。这些标准是在静态观测期间，通过每个学习统计量(即A，B，C，D)的时间序列(3-5分钟)的三倍标准差计算得到。Where t ₁ and t ₂ are the defined estimated time windows (generally defined as multiples of the sampling interval), is the number of epochs, is the mean value of a, A represents the displacement effect of the baseline drift, B represents the displacement effect of the original acceleration, C represents the displacement effect of the corrected acceleration, and D represents the signal strength. Each learning statistic has two cases according to the defined criteria A ₀ , B ₀ , C ₀ , D ₀ . These criteria are calculated by three standard deviations of the time series (3-5 min) of each learned statistic (i.e., A, B, C, D) during static observations.

基于所定义的统计量和标准，可以分析基线漂移估值并自适应识别。它们可以分成四个阶段，如图1所示。Baseline drift estimates can be analyzed and adaptively identified based on defined statistics and criteria. They can be divided into four stages, as shown in Figure 1.

初始化阶段(1)：满足(A₂&B₁)或(B₁&C₂)，图1.b中的0-47秒。数据解算的开始。静止状态，原始加速度几乎为零(如图1.a)，但解算的基线漂移不是零(如图1.b)。解算基线漂移需要时间来收敛，完成初始化。基线漂移的动态噪声应该是非常小的值。Initialization stage (1): satisfy (A ₂ &B ₁ ) or (B ₁ &C ₂ ), 0-47 seconds in Fig. 1.b. Start of data evaluation. At rest, the original acceleration is almost zero (as shown in Figure 1.a), but the calculated baseline drift is not zero (as shown in Figure 1.b). Solving baseline drift takes time to converge and complete initialization. The dynamic noise of the baseline drift should be a very small value.

静态阶段(2)：满足(A₁)，图1.b中的47-95秒。在此期间，解算的基线漂移几乎为零，基线漂移的动态噪声也非常小。Static phase (2): Satisfy (A ₁ ), 47-95 seconds in Fig. 1.b. During this period, the baseline drift of the solution is almost zero, and the dynamic noise of the baseline drift is also very small.

瞬时基线漂移阶段(3)：满足(A₂&D₂)，图1.b中的95-150秒。在此期间，检测到基线漂移，并且信号显著；基线漂移的动态噪声应调整为较大的值。Instantaneous Baseline Drift Phase (3): Satisfy (A ₂ & D ₂ ), 95-150 seconds in Fig. 1.b. During this period, baseline drift was detected and the signal was significant; the dynamic noise of the baseline drift should be adjusted to a larger value.

永久基线漂移阶段(4)：满足(A₂&D₁)&(B₂&C₁)，图1.b中的150秒至最后。在此期间，也检测到基线漂移，但没有显著信号，基线漂移的动态噪声应该是较小的值。Permanent baseline drift phase (4): satisfy (A ₂ &D ₁ ) & (B ₂ &C ₁ ), 150 seconds to the end in Fig. 1.b. During this period, baseline drift was also detected, but without a significant signal, the dynamic noise of the baseline drift should be a small value.

在检测到瞬时基线漂移的时间间隔内，动态噪声调整如下：During the time interval in which the transient baseline drift is detected, the dynamic noise is adjusted as follows:

${q q}_{u u}^{' '} = = \sqrt{\frac{| | {u u}_{{t t}_{22}} - - {u u}_{{t t}_{11}} | |}{(({t t}_{22} - - {t t}_{11}))}} - - - - - - ((88))$

其中q′_u为基线漂移的新的功率密度，反映了加速度基线漂移的实际变化率。另外，基线漂移估计窗口的长度可以根据信号存在情况自适应的调整：如果不存在信号，如情况D₁时，基线漂移主要是由环境变化引起，并且变化缓慢，所以本发明可以在几秒或更长的时间间隔来估计基线漂移；否则，如情况D₂，基线漂移主要是由地面运动的倾斜和/或旋转引起，且变化非常迅速，因此估计的频率应该为一秒或更少的时间间隔，来发现基线漂移的变化并调整动态噪声，本发明采用的是一秒钟估计一次。where q′ _u is the new power density of the baseline drift, which reflects the actual rate of change of the acceleration baseline drift. In addition, the length of the baseline drift estimation window can be adaptively adjusted according to the existence of the signal _: if there is no signal, such as in case D1, the baseline drift is mainly caused by environmental changes, and the change is slow, so the present invention can be within a few seconds or longer time intervals to estimate the baseline drift _; otherwise, as in case D2, the baseline drift is mainly caused by the tilt and/or rotation of the ground motion and changes very rapidly, so the estimated frequency should be one second or less The interval is used to find the change of the baseline drift and adjust the dynamic noise. The present invention uses an estimation once a second.

本发明的实施例包括以下步骤：Embodiments of the present invention include the following steps:

第一步，GPS与强震仪观测的紧组合处理，如图2所示。分为数据输入，数据处理和结果输出三个核心模块。在数据输入模块，需要提供GPS的相位/伪距观测值，强震仪的加速度观测值，以及GPS数据处理需要的辅助产品如轨道、钟差、地球自转参数产品。在数据解算模块，首先是数据的预处理，进行数据检查，粗差剔除和周跳探测得到干净的数据；随后，进行相对论、潮汐、天线相位中心、对流层、地球自转误差的修正工作；接着进行紧组合模型的建立，将强震仪的基线漂移误差当作未知参数和其它定位参数一起估计；最后进行残差编辑和精度统计得到组合解(即基本组合解)。The first step is the tight combination processing of GPS and strong motion instrument observations, as shown in Figure 2. It is divided into three core modules: data input, data processing and result output. In the data input module, it is necessary to provide the phase/pseudorange observation value of GPS, the acceleration observation value of the strong motion instrument, and the auxiliary products required for GPS data processing such as orbit, clock difference, and earth rotation parameter products. In the data calculation module, firstly, data preprocessing, data inspection, gross error elimination and cycle slip detection are performed to obtain clean data; then, relativity, tides, antenna phase center, troposphere, and earth rotation errors are corrected; then The establishment of the tight combination model is carried out, and the baseline drift error of the strong motion instrument is estimated as an unknown parameter together with other positioning parameters; finally, the combination solution (ie, the basic combination solution) is obtained by residual editing and precision statistics.

第二步，GPS与强震仪自适应组合处理，如图3所示。自适应组合处理采用两个卡尔曼滤波器来实现。第一个滤波器用于基本组合解算，即是第一步的工作，用来判断运动状态和确定最佳的动态噪声以及基线漂移估计窗口的长度。第二个滤波器用于自适应组合解算，以获得稳健的组合结果。由经验权重和动态噪声(这里经验权重由测量噪声来确定，通常的取值范围分别为：比如动态噪声给定为常数，比如q_u＝0.003-0.01m^1/2/s^3/2，窗口时间长度为1-10秒)，首先采取第一个滤波器来估计基线漂移，然后计算和确定其时间序列以及四种学习统计数据标准。通过分析，若判定为瞬时基线漂移阶段，动态噪声根据瞬时基线漂移的变化进行改变，基线漂移估计窗口的长度改变为每秒一次。对于其它状态期间，动态噪声限制为一个较小的值(取值范围为q_u＝0.0001-0.001m^1/2/s^3/2)，窗口长度保持为5秒。然后，根据确定的最佳动态噪声和时间估计窗口，使用第二个滤波器进行自适应组合解算，以获得最优组合结果。The second step is adaptive combination processing of GPS and strong motion instrument, as shown in Figure 3. The adaptive combination process is implemented using two Kalman filters. The first filter is used for the basic combined solution, that is, the first step, which is used to judge the motion state and determine the optimal dynamic noise and the length of the baseline drift estimation window. The second filter is used for adaptive combination solving to obtain robust combination results. By empirical weight and dynamic noise (here the empirical weight is determined by the measurement noise, the usual value ranges are: for example The dynamic noise is given as a constant, such as q _u =0.003-0.01m ^1/2 /s ^3/2 , and the window time length is 1-10 seconds), firstly adopt the first filter to estimate the baseline drift, and then calculate and determine Its time series as well as four learning statistics standards. Through analysis, if it is determined to be the phase of instantaneous baseline drift, the dynamic noise changes according to the change of instantaneous baseline drift, and the length of the baseline drift estimation window is changed to once per second. For other state periods, the dynamic noise is limited to a small value (the value range is q _u =0.0001-0.001m ^1/2 /s ^3/2 ), and the window length is kept at 5 seconds. Then, according to the determined optimal dynamic noise and time estimation window, the second filter is used for adaptive combination solution to obtain the optimal combination result.

Claims

1. A deformation monitoring method based on GPS and strong motion instrument observation adaptive combination, is characterized in that comprising the following steps:

The first step is to obtain the phase/pseudorange observation value of GPS, the acceleration observation value of strong motion instrument, and the auxiliary products required for GPS data processing, including GPS satellite orbit, satellite clock error and earth rotation parameters;

In the second step, the phase/pseudo-range observations of GPS obtained in the first step, the acceleration observations of the strong motion instrument, and the auxiliary products required for GPS data processing are carried out for data inspection, gross error elimination and cycle slip detection to obtain a clean data; corrections for relativistic, tidal, antenna phase center, tropospheric and earth rotation errors are applied to said clean data;

The third step is to establish a tight combination model of GPS and strong motion instrument observations. The tight combination model of GPS and strong motion instrument observations is as follows:

{[\begin{matrix} {L L}_{c c} \\ {P P}_{c c} \end{matrix}]}_{k k} = = [\begin{matrix} {e e}_{k k} & {m m}_{z z,, k k} & 11 & 11 \\ {e e}_{k k} & {m m}_{z z,, k k} & 11 & 00 \end{matrix}] [\begin{matrix} {x x}_{k k} \\ {z z}_{k k} \\ {δt δt}_{k k} \\ b b \end{matrix}] + + [\begin{matrix} {ϵ ϵ}_{{L L}_{c c}} \\ {ϵ ϵ}_{{P P}_{c c}} \end{matrix}],, \begin{matrix} {ϵ ϵ}_{{L L}_{c c}} ~ ~ N N ((00,, {σ σ}_{{L L}_{c c}}^{22})) \\ {ϵ ϵ}_{{P P}_{c c}} ~ ~ N N ((00,, {σ σ}_{{P P}_{c c}}^{22})) \end{matrix}

\begin{matrix} {a a}_{k k} = = {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} - - {u u}_{k k} + + {ϵ ϵ}_{{a a}_{k k}} & {ϵ ϵ}_{{a a}_{k k}} ~ ~ N N ((00,, {σ σ}_{{a a}_{k k}}^{22})) \end{matrix}

Among them, L _c and P _c are the carrier phase residuals and pseudorange code observation residuals after deducting the geometric distance between the satellite and the receiver; e is the unit vector between the satellite and the receiver; the subscript k represents the number of epochs ; m _z,k is the projection function of the tropospheric wet delay zenith direction; x, are the displacement and acceleration vectors of the receiver, respectively; z, δt, and b are the tropospheric zenith delay, receiver clock error, and phase ambiguity, respectively; ε is the measurement noise, and ε _Lc and ε _Pc represent the noise of carrier and pseudorange observations, respectively , and its corresponding variance is δ ² ; a and u are the acceleration observation value of the strong motion instrument and the baseline drift error respectively; the baseline drift error parameter is estimated as a constant in each estimation window of the set width, and the zenith tropospheric delay Estimated as a constant or manifested as a random walk process, the receiver clock error is estimated as Gaussian white noise epoch by epoch, and the phase ambiguity is estimated as a constant under continuous cycle-slip-free conditions;

The state equations of the receiver motion and baseline drift error are as follows:

[\begin{matrix} {x x}_{k k} \\ {\overset{\cdot \cdot}{x x}}_{k k} \\ {\overset{\cdot\cdot \cdot\cdot}{x x}}_{k k} \\ {u u}_{k k} \end{matrix}] = = [\begin{matrix} I I & 0.5 0.5 I I τ τ & 0.5 0.5 {Iτ Iτ}^{22} & 00 \\ 00 & I I & I I τ τ & 00 \\ 00 & 00 & I I & 00 \\ 00 & 00 & 00 & I I \end{matrix}] [\begin{matrix} {x x}_{k k - - 11} \\ {\overset{\cdot &Center Dot;}{x x}}_{k k - - 11} \\ {a a}_{k k} - - {\overset{^^}{u u}}_{k k - - 11} \\ {u u}_{k k - - 11} \end{matrix}] + + {w w}_{k k}

in, is the velocity vector of the receiver, I is the 3x3 identity matrix, τ is the minimum sampling rate of GPS and strong motion instrument, w _k is the dynamic noise, its mathematical expectation E(w _k )=0, variance q _u is the power density of baseline drift, k-1 means the previous epoch; the acceleration after baseline drift correction to replace the acceleration prediction value a _k-1 ; given the weight ratio of observations and dynamic noise; use Kalman filter to estimate parameters, and at the same time carry out residual editing and precision statistics to obtain the receiver's position, velocity, acceleration, baseline drift, Tropospheric delay and receiver clock information and accuracy indicators;

The fourth step is to calculate the time series of baseline drift based on the solution information in the third step, and construct four learning statistics to determine the standard of the statistics; the four statistics are:

Displacement Effects of Baseline Drift

Displacement Effect of Raw Acceleration

Corrected acceleration for displacement effects

signal strength

Among them, t ₁ and t ₂ are the start time and end time of the defined estimation window, are the epoch numbers corresponding to the start time and end time, respectively, is the mean value of a, the standard of each learning statistic A ₀ , B ₀ , C ₀ , D ₀ is three times the standard deviation of each learning statistic in the set time series during static observation;

In the fifth step, the baseline drift error parameters are divided into four stages according to the standard determined in the fourth step, and a new dynamic noise and time estimation window are calculated; the determination of the four stages is as follows:

Initialization phase: satisfy (A ₂ &B ₁ ) or (B ₁ &C ₂ ),

Static phase: satisfy (A ₁ ),

Instantaneous baseline drift phase: satisfy (A ₂ &D ₂ ),

Permanent baseline drift phase: satisfy (A ₂ &D ₁ )&(B ₂ &C ₁ );

During the transient baseline drift phase, the dynamic noise is adjusted as follows:

{q q}_{u u}^{' '} = = \sqrt{\frac{| | {u u}_{{t t}_{22}} - - {u u}_{{t t}_{11}} | |}{(({t t}_{22} - - {t t}_{t t}))}}

Among them, q′ _u is the new power density of the baseline drift; the length of the baseline drift estimation window is changed once per second; for the other three stages, the dynamic noise q _u =0.0001-0.001m ^1/2 /s ^3/2 , The window length is kept at 5 seconds;

In the sixth step, according to the new dynamic noise and the new time estimation window determined in the fifth step, Kalman filter is used for adaptive combination solution to obtain the optimal combination result.