Deformation Analysis of the Ms 6.2 Jishishan (Gansu,China) Earthquake on the Landslide Hazard Areas
-
摘要:
地震极易诱发或加速滑坡灾害发生,2023-12-18甘肃积石山县发生Ms 6.2地震,震源深度仅为10 km,周边滑坡易发区面临突发性破坏甚至提前发生灾害的风险,亟需对易发区变形进行快速分析与评估。基于此,采集了积石山地震远场不同距离(64 km、111 km、140 km和240 km)的4处滑坡隐患区全球导航定位系统(global navigation satellite system,GNSS)和加速度计实时监测数据,采用精密单点定位、实时动态差分定位技术、GNSS加速度计自适应融合技术以及小波变换等方法综合分析了主震和余震对边坡体的破坏性影响。研究结果表明,距震中64 km的黑方台滑坡体监测到1 cm的弹性位移以及0.5 cm的永久形变,距震中240 km的舟曲滑坡体监测到明显的加速度冲击事件,捕捉到的地表振动频率在0~10 Hz之间,加速度峰值达0.035 m/s²,增加了滑坡易发区的失稳风险。此次地震对滑坡造成的位移和振动影响在水平方向上较为显著。建议相关部门尽快对震源64 km范围内的滑坡易发区进行勘察和评估,避免突发性滑坡灾害的发生。
-
关键词:
- 2023年积石山地震 /
- 远场同震位移 /
- GNSS /
- 滑坡破坏评估 /
- 滑坡易发区
Abstract:ObjectivesEarthquake events can easily induce or accelerate the landslide failure. On 18th December 2023, an Ms 6.2 earthquake occurred in Jishishan County, Gansu province, with a focal depth of only 10 km. The surrounding landslide disaster-prone areas are faced with great hidden danger of sudden or advanced failure, and it is urgent to conduct rapid analysis and assessment of deformation in prone areas.
MethodsIn this paper, the global navigation satellite system(GNSS)and accelerometer real-time observation data of four landslide prone areas at different distances (64 km, 111 km, 140 km and 240 km) from the seismic center, and the damage effects of the main shock and aftershock on the landslide were comprehensively analyzed with precise point positioning, real-time kinematic, GNSS & accelerometer adaptive coupled technology and wavelet transform.
ResultsThe results show that the elastic displacement of 1 cm and the permanent displacement of 0.5 cm aredetected in the Heifangtai landslide, 64 km away from the seismic center, while the obvious impact responses are detected in the Zhouqu landslide, 240 km away from the seismic center. The vibration frequency was between 0‐10 Hz, and the peak acceleration reached 0.035 m/s², which increased the risk of landslide failure.
ConclusionsThe displacement and vibration effects of the earthquake on the landslide are mainly concentrated in the horizontal direction, and the vertical direction is not significant. It is suggested the landslide disaster-prone areas within 64 km of Jishishan County should be investigated and evaluated as soon as possible to prevent secondary disaster damage.
-
全球导航卫星系统(global navigation satellite system,GNSS)分析中心连续向用户提供最终、快速和超快速精密轨道产品。其中,超快速轨道作为实时或近实时定位的重要参数,其精度将直接影响模糊度固定[1]和定位结果[2]。因此,分析中心对所提供的超快速轨道精度有严格的精度指标要求。例如,国际GNSS监测评估系统(international GNSS monitoring and assessment system,iGMAS)对多系统各类卫星的超快速轨道观测与预报部分精度做了详细规定。通过评估可知,GPS超快速预报轨道6 h和24 h的三维均方根误差(three-dimensional root mean square error, 3D RMS)分别达到41.7 mm与80.2 mm,与最终产品精度有明显差距,无法充分满足高精度GNSS用户的需求。为提高超快速轨道精度,学者们从预报策略[3]、最优弧段长度[4]、预报时间间隔[5]和地球自转参数误差影响[6]等方面对超快速预报轨道进行了精化。但是,这些研究主要针对的是超快速轨道的预报部分,对于观测部分的轨道精度改进还未具体展开。随着GNSS的蓬勃发展,如北斗“三步走”策略,针对各类卫星的高精度超快速轨道研究是拓展GNSS快速高效服务的必要前提。
多频多模GNSS发展、逐渐增加的卫星数目以及不断扩展的跟踪网,给分析中心超快速轨道时效性带来巨大挑战[7]。因此,时效性同样是超快速轨道生成过程中需要考虑的因素之一。为有效提高GNSS定轨参数解算效率,不同的数据处理方法,如参数消除法[8]、双差模糊度恢复法[9]、不动点理论模糊度快速固定[10]和增加参数解算间隔[11]等被详细探讨和验证,以及不同的测站优化模型[7]被提出,以解决测站分布与定轨时效性之间的关系。然而这些方法仍存在明显的缺陷,如无法顾及参数之间的相关性、不能保证产品的一致性以及无准确的函数模型等等,并不能从本质上解决分析中心超快速轨道时效性的问题。
由定轨观测方程可知,超快速轨道时效性和精度与定轨观测数据(质量、分布以及数量)直接相关[9],而分析观测数据对定轨参数解算精度的贡献主要有两个指标,即观测值精度与参数精度衰减因子(dilution of precision,DOP)[12]。其中,DOP值与定轨空间构型直接相关[13]。基于DOP值的定位最优空间构型已经得到广泛且深入的研究,由于其将地表定位最优简化为地心处(DOP值最小),因此相关理论并不适用于精密轨道确定模型。基于DOP值的地面跟踪站分布优化,有学者研究了测定地球自转参数[14]和北斗地球同步轨道(geostationary orbit satellite,GEO)卫星轨道[15]参数的最优构型。但是,这些基于特殊条件的研究不具有普适性。而实验发现,地面测站数目及空间位置对卫星轨道及其相关参数解算起到至关重要的作用[16]。因此,对于定轨弧段中某历元,其轨道空间状态参数(位置与速度)受观测数据的影响可表示为跟踪站与卫星构型的函数关系。实验表明[17],通过建立卫星轨道参数精度与其DOP值之间的函数关系,可实现对观测数据受限条件下的超快速轨道精度的间接改善,尤其是对北斗轨道精度的提升更加明显。同时,对于数据冗余区域,建立以DOP值最小为准则的跟踪站分布优化模型(以较少的测站实现精度相当的定轨过程),可顾及超快速轨道解算的精度与时效性[7]。所以,充分利用定轨DOP值是提高超快速轨道精度与时效性的有效途径之一。
本文主要针对超快速轨道观测后期由于观测数据不足而导致轨道精度降低的问题,建立超快速轨道精度与其状态参数DOP值的函数模型;通过高精度预报DOP值变化量,修正超快速轨道观测后期精度降低的问题,进一步提高预报轨道精度。
1 GNSS超快速观测轨道精化模型
1.1 基于DOP值轨道精化原理
由文献[7, 17]可知,轨道状态参数DOP值与轨道精度存在一定的相关性。因此,在超快速轨道观测后期,由于观测数据缺失而导致的轨道精度降低可通过定轨DOP值间接地改善。具体如下:设历元ti定轨观测方程(为研究轨道状态参数与DOP之间的关系,将除轨道以外的参数作为已知值,则待定参数仅为轨道状态参数)表示为:
$$\mathit{\boldsymbol{V}}\left( {{t_i}} \right) = \mathit{\boldsymbol{L}}\left( {{t_i}} \right) - \mathit{\boldsymbol{A}}\left( {{t_i}} \right)\mathit{\boldsymbol{X}}\left( {{t_i}} \right)$$ (1) 式中,L(ti)为观测值;V(ti)为改正量;X(ti)为轨道状态参数;A(ti)为系数矩阵。式(1)中,高度角大于30°采用等权模型,10°~30°选择高度角模型,剔除小于10°的观测数据,其系数矩阵可进一步表示为:
$$\mathit{\boldsymbol{A}}\left( {{t_i}} \right) = {\rm{diag}}\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{a}}_1}}&{{\mathit{\boldsymbol{a}}_2}}& \cdots &{{\mathit{\boldsymbol{a}}_{m - 1}}}&{{\mathit{\boldsymbol{a}}_m}} \end{array}} \right)$$ (2) $$ {\mathit{\boldsymbol{a}}_m} = \left[ \begin{array}{l} \begin{array}{*{20}{c}} {\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial z_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial \dot x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial \dot y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 1}^{{s_m} - {r_1}}}}{{\partial \dot z_m^{{s_m}}}}}\\ {\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial z_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial \dot x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial \dot y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, 2}^{{s_m} - {r_2}}}}{{\partial \dot z_m^{{s_m}}}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{array}\\ \begin{array}{*{20}{c}} {\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial z_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot x_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot y_m^{{s_m}}}}}&{\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot z_m^{{s_m}}}}} \end{array} \end{array} \right] $$ (3) 式中,m、n分别表示ti历元中卫星与测站数目;$\left( {\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial x_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial y_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial z_m^{{s_m}}}}} \right)$与$\left( {\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot x_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot y_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot z_m^{{s_m}}}}} \right)$分别表示sm卫星与rn测站间的几何距离关于卫星位置与速度的偏导数[17], 则:
$$\begin{array}{l} \left( {\frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial x_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial y_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial z_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot x_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot y_m^{{s_m}}}}, \frac{{\partial \rho _{m, n}^{{s_m} - {r_n}}}}{{\partial \dot z_m^{{s_m}}}}} \right) = (\frac{{x_m^{{s_m}} - x_n^{{r_n}}}}{{\rho _{m, n}^0}}, \frac{{y_m^{{s_m}} - y_n^{{r_n}}}}{{\rho _{m, n}^0}}, \frac{{z_m^{{s_m}} - z_n^{{r_n}}}}{{\rho _{m, n}^0}}, \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{\rm{\Delta }}t\left( {x_m^{{s_m}} - x_n^{{r_n}}} \right)}}{{\rho _{m, n}^0\left( {{t_r}} \right)}}, \frac{{{\rm{\Delta }}t\left( {y_m^{{s_m}} - y_n^{{r_n}}} \right)}}{{\rho _{m, n}^0\left( {{t_r}} \right)}}, \frac{{{\rm{\Delta }}t\left( {z_m^{{s_m}} - z_n^{{r_n}}} \right)}}{{\rho _{m, n}^0\left( {{t_r}} \right)}}) \end{array}$$ (4) 式中,tr为信号接收时刻;Δt表示传播时间。由式(1)可得参数协因数阵为:
$$\mathit{\boldsymbol{Q}}({t_i}) = {({\mathit{\boldsymbol{A}}^T}({t_i})\mathit{\boldsymbol{PA}}({t_i}))^{( - 1)}}$$ (5) 式中,P为权阵。第k个参数先验解算精度为:
$${\sigma _k} = {\sigma _0} \cdot \sqrt {\mathit{\boldsymbol{Q}}{{({t_i})}_{kk}}} $$ (6) 式中,σ0为中误差。由式(6)和文献[17]可知,参数解算精度与其对应的协因素之间存在一定关系:式(6)在不考虑σ0的条件下,对应参数DOP值可以衡量由测站与卫星形成的几何构型对解算参数的影响,其可表示为参数协因素矩阵对角线元素之和的平方根。因此,第k个轨道状态参数DOP值可表示为:
$$D{({t_i})_k} = \sqrt {\mathit{\boldsymbol{Q}}{{({t_i})}_{kk}}} $$ (7) 为构建以轨道状态参数DOP值为自变量的轨道精化函数模型,需要实现DOP值与轨道参数改正数一一对应,而式(7)无法确保任意历元DOP值不相等的条件成立。因此,本文定义总体DOP值,以区分各历元DOP值,即当前历元ti第k个参数的观测方程总体DOP值为前ti个DOP值平方和的平方根:
$$[D\left( {{t_i}{)_k}} \right] = \sqrt {{{[D{{({t_{i - 1}})}_k}]}^2} + {{\left( {D{{({t_i})}_k}} \right)}^2}} $$ (8) 式中,[·]表示总体DOP值。假设将卫星轨道状态改正dX表示为[D(ti)]的函数,即:
$$ {\rm{d}}\mathit{\boldsymbol{X}}\left( {{t_i}} \right) = \left[ {\begin{array}{*{20}{c}} {{\rm{d}}{X_1}}\\ {{\rm{d}}{X_2}}\\ \vdots \\ {{\rm{d}}{X_6}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {f(\left[ {D\left( {{t_i}{)_1}} \right]} \right)}\\ {f(\left[ {D\left( {{t_i}{)_2}} \right]} \right)}\\ \vdots \\ {f(\left[ {D\left( {{t_i}{)_6}} \right]} \right)} \end{array}} \right] $$ (9) 式中,f (·)表示轨道状态改正函数。基于GNSS定轨观测方程,每个历元DOP值可由式(5)精确计算(总体DOP值可相应获得)。因此,可用函数表示轨道状态改正数与总体DOP值之间的关系。由于状态改正函数受建模数据长度、采样间隔、拟合长度以及改正精度等因素共同影响(动力学模型相同的情况下),轨道修正算法中可构建以赤池信息量准则(Akaike information criterion,AIC)为基础的轨道状态改正模型优选方法。为便于讨论本文的轨道精化模型,以多项式函数为例,构建以总体DOP值为自变量的轨道状态参数精化模型,即:
$$\left[ {\begin{array}{*{20}{c}} {f(\left[ {D\left( {{t_i}{)_1}} \right]} \right)}\\ {f(\left[ {D\left( {{t_i}{)_2}} \right]} \right)}\\ \vdots \\ {f(\left[ {D\left( {{t_i}{)_6}} \right]} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {[D\left( {{t_i}{)_1}} \right]}&{{{[D{{({t_i})}_1}]}^2}}& \cdots &{{{[D{{({t_i})}_1}]}^q}}&1\\ {[D\left( {{t_i}{)_2}} \right]}&{{{[D{{({t_i})}_2}]}^2}}& \cdots &{{{[D{{({t_i})}_2}]}^q}}&1\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ {[D\left( {{t_i}{)_6}} \right]}&{{{[D{{({t_i})}_6}]}^2}}& \cdots &{{{[D{{({t_i})}_6}]}^q}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\alpha _1}}\\ {{\alpha _2}}\\ \vdots \\ {{\alpha _q}}\\ {{\alpha _0}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\xi _1}}\\ {{\xi _2}}\\ \vdots \\ {{\xi _6}} \end{array}} \right]$$ (10) 式中,α表示多项式系数;ξ为拟合残差;q为拟合阶数。通过精确求解多项式系数,可准确建立轨道精度与总体DOP值之间的函数关系式。需要注意的是,式(10)求得的轨道改正量为绝对量,通过轨道变化趋势可确定实际改正量[6]。针对超快速观测轨道后期精度降低,式(10)是以总体DOP值为自变量的函数表达式,因此通过精确获得超快速轨道后期DOP值,可间接求取轨道改正量。
1.2 DOP值精确预报方法
在测站均匀分布条件下,定轨观测方程中参数对应DOP值呈平滑趋势[17]。为描述总体DOP值的变化趋势,本文基于多项式模型拟合总体DOP变化。模型中历元为自变量,模型阶数需要根据不同参数进行实时确定。因此,设总体DOP的变化趋势为:
$$D'{({t_j})_k} = {\theta _0} + {\theta _1}{t_j} + \cdots + {\theta _b}t_j^b + e\left( {{t_j}} \right)$$ (11) 式中,θ为总体DOP曲线系数;e为模型误差;tj表示历元标记;b为多项式阶数,其由观测数据采样间隔、预报弧长和预报精度共同决定[3]。则式(11)中的矩阵形式为:
$$\mathit{\boldsymbol{D}}_k^{\rm{'}} = \mathit{\boldsymbol{G}} \cdot \mathit{\boldsymbol{\theta }} + {\rm{d}}\mathit{\boldsymbol{D}}_k^{\rm{'}}$$ (12) 式中,G为系数矩阵;dDk′表示拟合残差;θ = [θ0 θ1 θ2 … θb]T。则:
$$\mathit{\boldsymbol{\widehat \theta }} = {\left( {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}}} \right)^{ - 1}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{D}}_k^{\rm{'}}$$ (13) 令当前历元tj = 0,则预报D(tj)k为:
$$D{({t_j})_k} = {\hat \theta _0}$$ (14) 将式(14)中预报的DOP值代入式(10)中,求解当前历元轨道改正量,对轨道进行改正。
2 超快速精密轨道确定实验
2.1 超快速轨道观测部分精度分析
国际GNSS服务组织(International GNSS Service, IGS)从2000年12月正式提供12 h间隔的GPS超快速卫星轨道服务。自2004年,更新间隔提高至6 h,由48 h弧段的观测(24 h)与预报(24 h)两部分组成。目前,GNSS用户可获得IGS及其分析中心3 h延迟的GPS与GLONASS组合超快速星历;而iGMAS可提供包括BeiDou和Galileo在内的四系统超快速精密星历。为分析超快速轨道精度,取iGMAS分析中心武汉大学(Wuhan University, WHU)提供的超快速轨道以及IGS分析中心德国地学研究中心(German Research Center for Geosciences, GFZ)提供的多系统快速星历作为参考,进行精度分析。
以2017年年积日(day of year,DOY)第168天超快速轨道(whu19540_00.sp3)为例,相对应的快速多系统轨道(gbm19536.sp3)作为参考。首先,统计超快速轨道观测部分与GFZ轨道之间的残差(利用赫尔默特参数进行基准转换),其观测后期(最后3 h)对应的轨道3D RMS如图 1所示(图中G表示GPS,R表示GLONASS, E表示Galileo,C表示BeiDou)。然后,针对超快速轨道预报部分,由于预报轨道采用积分方式获取,其误差随积分时间长度的增加而逐渐增加;同时,在计算预报轨道初始状态时,观测轨道部分后期占有较大权重。因此,在分析了观测轨道后期误差对预报轨道的影响后,本文取超快速轨道观测最后1 h的数据进行轨道拟合,利用拟合后的初值预报24 h多系统轨道,并与GFZ星历进行精度对比,统计前6 h轨道(以2 h为间隔)。重复上述实验步骤,连续统计一个月(2017年DOY第168~197天)的实验结果,得到超快速观测轨道与预报轨道的精度(3D RMS)平均值分别如表 1、表 2所示。
表 1 超快速轨道观测部分残差3D RMS/cmTable 1. 3D RMS of the Observed Ultra-rapid Orbit/cm卫星系统 1~20 h 21 h 22 h 23 h 24 h GPS 3.6 3.5 4.6 6.3 7.5 GLONASS 6.1 5.6 6.1 7.5 10.2 BeiDou 12.9 12.1 12.4 15.6 21.8 Galileo 8.2 13.2 13.5 14.4 17.9 表 2 超快速轨道预报部分3D RMS/cmTable 2. 3D RMS of the Predicted Ultra-rapid Orbit/cm卫星系统 2 h 4 h 6 h 1~12 h 1~24 h GPS 7.8 9.8 10.6 10.8 16.1 GLONASS 12.5 13.4 13.6 14.3 21.9 BeiDou 23.2 39.9 61.1 70.1 139.9 Galileo 16.9 26.6 32.5 34.4 49.9 上述实验中并不是所有卫星都会出现精度明显降低的现象,图 1中仅绘制了精度出现降低的相关卫星。通过对超快速轨道观测与预报部分进行精度统计,可以得出:(1)超快速轨道观测部分的后期(最后3 h)精度出现明显降低现象。为分析超快速观测轨道后期精度降低现象,本文以分析中心可获得的一个月(2016年6月)稳定的409个测站为基础,基于自编GNSS数据预处理软件,统计了所有测站的数据质量。表 3统计了所有测站的观测数据前21 h和后3 h的数据质量的平均值,从表 3中可以看出一天内的观测数据质量无明显差别。(2)利用超快速观测部分的后期轨道(最后1 h)进行轨道预报,其会明显降低预报轨道精度。综上,为提高超快速轨道精度,有必要对超快速轨道观测部分的后期数据进行修正,改善超快速轨道精度。
表 3 连续一个月观测数据质量对比Table 3. Comparison of Data Quality During One Month指标 1~21 h 最后3 h 最小值 最大值 平均值 最小值 最大值 平均值 数据完整率/% 95.80 100.00 97.19 95.20 100.00 96.27 MP1 /m 0.03 0.29 0.21 0.05 0.28 0.21 MP2 /m 0.04 0.32 0.27 0.09 0.41 0.36 SN1 6.43 50.85 42.17 6.26 50.63 41.77 SN2 4.42 46.29 34.88 4.31 46.05 34.34 周跳比 0.02 4.22 1.49 0.02 4.81 2.48 注:MP1与MP2分别表示不同卫星系统波段1与波段2上的多路径效应平均值;SN1与SN2分别表示波段1与波段2信噪比平均值 2.2 基于定轨DOP值的轨道精化策略
本文取分析中心可稳定下载的小时文件合并得到天文件,并利用合并后的观测文件进行超快速轨道计算与预报。具体实验过程如图 2所示,主要步骤如下:(1)获取导航文件、测站列表和测站坐标,合并小时观测文件(最后3 h不合并);并对观测文件进行预处理,筛选定轨测站。(2)逐历元计算对应各参数的DOP值,并将其逐历元累加(获取各参数总DOP值)。(3)进行超快速轨道计算与预报,并将观测部分轨道与GFZ对应的多系统轨道进行对比,计算轨道残差。(4)建立各卫星状态参数与其对应历元总体DOP值的函数模型。(5)基于§2.2中模型预报超快速轨道观测部分的后期(最后3 h)DOP值。(6)将预报的DOP值代入轨道修正函数模型,利用函数模型对观测部分后期轨道进行改正,并基于改进后的单天弧段观测轨道进行单天弧段轨道预报。
2.3 基于DOP值的轨道修正模型分析
为验证基于DOP值超快速观测轨道精化模型的可行性,本文对DOP值拟合模型与轨道状态参数精化函数模型进行分析。首先,分析轨道参数DOP值与观测数据量之间的关系。选取分析中心可获得的一个月(2016年6月)稳定的409个测站,采样间隔为30 s的跟踪站观测数据。图 3(a)与图 3(b)分别表示409个全球跟踪站分布图与GFZ定轨所用测站分布图,图中用不同颜色表示可接收单系统(G: GPS)、双系统(G+R: GPS+GLONASS)、三系统(G+R+E: GPS+GLONASS+Galileo)以及四系统(G+R+E+C: GPS+GLONASS+Galileo+BeiDou)的测站分布。分别保留任意选取的200、150、100、50和0个(不保留)测站最后3 h的小时观测数据进行DOP值计算及相应定轨精度统计。由于定轨实验数据量大,图 4中仅画出具有代表的G09和C13卫星对应的一天(2016年DOY第153天)内DOP值及轨道精度变化趋势(以409个测站计算出的轨道为参考;由于200个测站方案与150个测站方案的结果十分接近,为便于区分,图中将200个测站方案轨道精度省略)。通过不同定轨方案可知:(1)随着观测部分后期测站数目减少(最后3 h的小时观测文件数目逐渐减少),DOP值出现明显增大趋势;同时,轨道精度与DOP值变化趋势保持一致。(2)轨道后期的观测值会影响整个观测弧段的定轨精度,这主要由于观测数据减少影响了轨道初始空间状态参数。(3)整个定轨弧段内,轨道状态参数DOP值呈平滑曲线。
然后,为准确描述定轨DOP值的变化趋势,本文分别基于任意选取的409、200、150、100和50个测站观测数据计算一天定轨弧段内轨道参数DOP值,并采用AIC确定最优多项式阶数。AIC由以下各项定义:
$$C = n{\rm{log}}\left( {{{\hat \sigma }^2}} \right) + 2k$$ (15) 式中,n是观测的数目;$\hat \sigma $是均方误差;k是参数的数目。在式(15)的等号右侧,第一项表示模型的优良性,第二项表示模型复杂性。因此,优选模型是具有最小AIC值的模型。以G09为例(2016年DOY第153天),其相应的模型信息量统计见表 4。从表 4中可以得出,当阶数为4时,赤池信息量最小,其拟合模型最优。在本文提出的基于定轨DOP值轨道精化模型中,应根据不同卫星选择最优的DOP预报模型阶数。
表 4 多项式模型赤池信息量Table 4. AIC of Polynomial Model阶数 409个测站 200个测站 150个测站 50个测站 1 51.11 49.70 51.43 57.09 2 36.64 35.52 34.15 48.18 3 27.48 26.41 21.22 28.32 4 18.14 16.50 15.10 18.34 5 21.51 30.40 19.42 28.22 6 28.34 30.18 25.13 38.33 同时,为讨论本文提出的轨道精化模型的可靠性,实验从观测数据长度和函数模型两个方面对状态改正模型进行验证。由§2.1实验表明,超快速轨道观测部分后期出现轨道发散现象,因此实验主要预报观测部分最后3 h的轨道改正数。具体实验方案如下。
方案1:分别基于0.5、1、2和3 d的观测数据(采样间隔30 s),预报3 h的轨道改正数(轨道间隔5 min);统计连续一个月(2016年DOY第122~151天)的预报精度;式(10)中多项式阶数采用赤池信息准则确定。相应的轨道精度统计结果(以全部测站计算轨道为参考)如图 5所示。
方案2:基于单天观测数据分别利用多项式、灰色模型、神经网络模型以及自回归模型构建轨道改正模型;同样预报3 h轨道改正数,并统计了连续一个月轨道修正精度。其相应的统计结果如图 6所示。
方案1中基于多项式函数模型,当观测数据长度为1 d时,对轨道修正效果最优;而当数据长度为0.5 d时,轨道修正效果最差。通过4组对比实验可发现:(1)本文提出的超快速轨道修正模型受观测数据长度影响,其主要是由于不同的数据长度条件下,函数修正模型拟合程度存在差异(求解预报轨道初始状态参数存在差异)。(2)轨道无法实现完全修正,由于定轨受摄动力模型以及参数估计模型共同影响,基于DOP值修正模型仅实现部分修正。(3)为了顾及超快速轨道的时效性与修正精度,本文采用1 d观测数据作为修正模型建模序列。
方案2针对不同的预报模型,从修正结果(见图 6)中可发现,神经网络模型优于其他3种预报模型,这主要是由于其他3种为线性模型。同时,4种模型修正效果无明显差异,因此在构建轨道改正模型时,可根据具体需要构建合适的函数模型。在后续研究中,将开展基于方差分量估计以及赤池信息准则的改正模型优选算法研究。
2.4 超快速观测轨道修正实验
为具体分析本文提出的超快速轨道精化模型,实验选取连续10 d(2016年DOY第141~150天)的轨道结果进行对比分析。其中轨道精化函数模型与DOP值预报模型分别采用§2.3与§1.2中的多项式模型,模型阶数均采用AIC确定。图 7中列出了DOY第141天对应的C13、E19、G09和R11的轨道修正前后精度(以GFZ轨道为参考,其计算轨道测站分布见图 3(b))和DOP值(预报值与计算值)。从实验结果可以得出,基于DOP值超快速轨道修正方法能改善观测数据缺失条件下的轨道精度(图 7(b)统计了观测部分后3 h轨道精度)。连续10 d各系统轨道修正结果如图 8所示,统计结果见表 5。
表 5 连续10天超快速轨道修正前后3D RMS对比Table 5. Comparison of 3D RMS Before and After Ultra-rapid Orbit Correction for Ten Comsecutive Days卫星类型 2016年DOY第141~150天3D RMS/mm 提升率/% 141天 142天 143天 144天 145天 146天 147天 148天 149天 150天 GPS(B) 42 43 45 48 39 48 36 42 48 42 GPS(A) 34 31 39 36 34 38 28 37 34 37 19.6 GLONASS(B) 59 56 67 61 66 66 72 61 77 78 GLONASS(A) 44 50 65 58 58 62 64 55 26 35 22.0 Galileo(B) 81 88 76 74 88 83 74 83 90 88 Galileo(A) 73 78 68 63 79 75 63 66 79 76 12.7 BDS_MEO(B) 80 75 88 67 86 68 66 79 88 88 BDS_MEO(A) 65 66 80 56 72 59 60 74 71 85 12.4 BDS_IGSO(B) 85 78 88 70 88 88 78 79 89 74 BDS_IGSO(A) 66 64 68 62 72 74 71 71 76 63 15.9 注:A表示修正后;B表示修正前;BDS(BeiDou navigation satellite system)代表北斗卫星导航系统; MEO(medium earth orbit)为中圆地球轨道卫星; IGSO(inclined geosynchronous orbit)为倾斜地球同步轨道卫星 通过本文提出的基于DOP值的超快速观测轨道后期修正模型,在一定程度上可以改善超快速轨道,其轨道观测部分后期(最后3 h)精度提升了12.4% ~22.0%;但是并不能完全对观测数据缺失情况下轨道误差进行修正,主要原因有两点:(1)精密定轨受力学模型、参数解算策略和模糊度固定率等因素综合影响,轨道参数DOP值与其精度之间相关性无法精确获得。(2)随着观测轨道后期弧段逐渐增加,修正函数模型误差逐渐增大,进一步限制了轨道修正效果。上述模拟实验中是将分析中心可下载的测站全加入定轨方程,保证了在已有测站条件下的定轨空间构型最优条件。然而,由于超快速轨道时效性限制,无法在定轨过程中将所有测站加入定轨过程中,导致定轨跟踪站与卫星之间的几何构型无法达到最优或次优,降低DOP与轨道精度之间相关性,限制了基于DOP超快速轨道修正效果;因此,超快速轨道计算过程中,测站最优或次优分布必须考虑。
3 结语
超快速轨道作为GNSS用户精密定位所需要的重要产品,其精度直接影响实时或近实时的应用。然而,通过轨道分析发现,超快速轨道观测部分后期(3 h)轨道精度逐渐降低,其将进一步影响超快速预报轨道精度。为了修正由于观测数据缺失而导致的超快速轨道观测部分发散的问题,本文从定轨空间构型角度出发,研究了观测数据不足条件下的DOP值与定轨精度之间的关系。实验结果表明,超快速观测轨道后期观测数据缺失所引起的DOP变化与轨道精度变化趋势相同。因此,本文提出了一种基于DOP值的轨道精化模型。通过对比实验,首先验证了DOP值预报模型的可行性,同时基于赤池信息准则构建了DOP预报模型最优阶数筛选方法;然后分析了不同函数模型的轨道结果,并对比了4组不同的函数模型,发现神经网络模型优于线性模型,但修复效果无明显差异;最后分析了连续10 d基于DOP值的超快速轨道观测部分的后期修正效果,以DOP为自变量的轨道修正方法可以改善由于观测数据缺失而引起的轨道精度降低现象,精度提升了12.4%~22.0%。由于轨道精度受定轨力学模型、参数解算策略以及模糊度固定模式等诸多因素影响,而空间构型不能完全对轨道精度进行表达,因此轨道精度修正效果受到了一定限制。后期将基于本文的研究探讨各因素对定轨的影响,精化超快速轨道修正模型。
http://ch.whu.edu.cn/cn/article/doi/10.13203/j.whugis20230490 -
-
[1] FENG G C, LI Z W, SHAN X J, et al. Geodetic Model of the 2015 April 25 Mw 7.8 Gorkha Nepal Earthquake and Mw 7.3 Aftershock Estimated from InSAR and GPS Data[J]. Geophysical Journal International, 2015, 203(2): 896-900.
[2] 张晓超, 裴向军, 张茂省, 等. 强震触发黄土滑坡流滑机理的试验研究: 以宁夏党家岔滑坡为例[J]. 工程地质学报, 2018, 26(5): 1219-1226. ZHANG Xiaochao, PEI Xiangjun, ZHANG Mao-sheng, et al. Experimental Study on Mechanism of Flow Slide of Loess Landslides Triggered by Strong Earthquake—A Case Study in Dangjiacha, Ningxia Province[J]. Journal of Engineering Geology, 2018, 26(5): 1219-1226.
[3] 王家鼎, 张倬元. 地震诱发高速黄土滑坡的机理研究[J]. 岩土工程学报, 1999, 21(6): 670-674. WANG Jiading, ZHANG Zhuoyuan. A Study on the Mechanism of High Speed Loess Landslide Induced by Earthquake[J]. Chinese Journal of Geotechnical Engineering, 1999, 21(6): 670-674.
[4] ZHUANG W Q, CUI D X, HAO M, et al. Geodetic Constraints on Contemporary Three-Dimensional Crustal Deformation in the Laji Shan–Jishi Shan Tectonic Belt[J]. Geodesy and Geodynamics, 2023, 14(6): 589-596.
[5] 李智敏, 李延京, 田勤俭, 等. 拉脊山断裂古地震与喇家遗址灾变事件关系研究[J]. 地震研究, 2014, 37(S1): 109-115. LI Zhimin, LI Yanjing, TIAN Qinjian, et al. Study on the Relationship Between Paleoseismic on Laji Mountain Fault and Catastrophic Event on Lajiashan Site[J]. Journal of Seismological Research, 2014, 37(S1): 109-115.
[6] 中国科学院青藏高原研究所. 甘肃临夏M 6.2级地震震源破裂过程反演初步结果[EB/OL]. (2023-12-19). https://mp.weixin.qq.com/s/SuKbTgq8-0EV2Fb NcdlYqQ. Institute of Tibetan Plateau Research, Chinese Aca-demy of Sciences. Preliminary Results of the Inversion of the Rupture Process of the Seismic Source of the Gansu Linxia M 6.2 Magnitude Earthquake[EB/OL]. (2023-12-19). https://mp.weixin.qq.com/s/SuKbTgq8-0EV2FbNcdlYqQ.
[7] 复合链生自然灾害. 同震滑坡发生概率!2023年12月18日甘肃临夏积石山县Ms 6.2地震[EB/OL]. (2023-12-20). https://mp.weixin.qq.com/s/D4nDhPw_MAoNN54YBGsQ4g. National Institute of Natural Hazards. Probability of Coseismic Landslides! December18, 2023 Ms 6.2 Earthquake in Jishishan County, Gansu, China[EB/OL]. (2023-12-20). https://mp.weixin.qq.com/s/D4nDhPw_MAoNN54YBGsQ4g.
[8] 清华大学陆新征课题组. RED-ACT |12月18日夜甘肃6.2级地震破坏力分析[EB/OL]. (2023-12-19). https://mp.weixin.qq.com/s/oZJDAaFtOPm33mF38fRIKA. Lu Xinzheng Research Group, Tsinghua University. RED-ACT |Analysis of the Destructive Power of the6.2-Magnitude Earthquake in Gansu on the Night of December 18th [EB/OL]. (2023-12-19). https://mp.weixin.qq.com/s/oZJDAaFtOPm33mF38fRIKA.
[9] 地质灾害防治与地质环境保护国家重点实验室. 甘肃积石山6.2级地震诱发地质灾害空间概率预测[EB/OL]. (2023-12-19). https://www.sklgp.cdut.edu.cn/info/1025/6599.htm. State Key Laboratory of Geohazard Prevention and Geoenvironment Protection. Spatial Probability Prediction of Geologic Hazards Induced by the 6.2 Magnitude Earthquake in Jishishan, Gansu Province [EB/OL]. (2023-12-19). https://www.sklgp.cdut.edu.cn/info/1025/6599.htm.
[10] 湖北地震局. 2023年12月18日甘肃临夏M 6.2地震GNSS快速响应[EB/OL]. (2023-12-20). https://mp.weixin.qq.com/s/xvZUVy4rHqeAql EzeulY6A. Hubei Earthquake Agency. GNSS Rapid Response to the December18, 2023 Gansu Linxia M 6.2 Earthquake [EB/OL]. (2023-12-20). https://mp.weixin.qq.com/s/xvZUVy4rHqeAqlEzeulY6A.
[11] FAN X M, SCARINGI G, KORUP O, et al. Earthquake-Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts[J]. Reviews of Geophysics, 2019, 57(2): 421-503.
[12] WU Z L, MA T F, JIANG H, et al. Multi-scale Seismic Hazard and Risk in the China Mainland with Implication for the Preparedness, Mitigation, and Management of Earthquake Disasters: An Overview[J]. International Journal of Disaster Risk Reduction, 2013, 4: 21-33.
[13] 张勤, 白正伟, 黄观文, 等. GNSS滑坡监测预警技术进展[J]. 测绘学报, 2022, 51(10): 1985-2000. ZHANG Qin, BAI Zhengwei, HUANG Guanwen, et al. Review of GNSS Landslide Monitoring and Early Warning[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(10): 1985-2000.
[14] HUANG G W, DU S, WANG D. GNSS Techniques for Real-Time Monitoring of Landslides: A Review[J]. Satellite Navigation, 2023, 4(1): 5.
[15] LIU X H, DU Y, HUANG G W, et al. Mitigating GNSS Multipath in Landslide Areas: A Novel Approach Considering Mutation Points at Different Stages[J]. Landslides, 2023, 20(11): 2497-2510.
[16] 白正伟, 张勤, 黄观文, 等. “轻终端+行业云” 的实时北斗滑坡监测技术[J]. 测绘学报, 2019, 48(11): 1424-1429. BAI Zhengwei, ZHANG Qin, HUANG Guanwen, et al. Real-Time BeiDou Landslide Monitoring Technology of “Light Terminal Plus Industry Cloud”[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(11): 1424-1429.
[17] 杜源, 王纯, 张勤, 等. 顾及黄土滑坡灾害状态特征的实时GNSS滤波算法[J]. 武汉大学学报(信息科学版), 2023, 48(7): 1216-1222. DU Yuan, WANG Chun, ZHANG Qin, et al. Real-Time GNSS Filtering Algorithm Considering State Characteristics of Loess Landslide Hazards[J]. Geomatics and Information Science of Wuhan University, 2023, 48(7): 1216-1222.
[18] 魏正发, 张俊才, 曹小岩, 等. 青海西宁南北山滑坡、崩塌成因及影响分析[J]. 中国地质灾害与防治学报, 2021, 32(4): 47-55. WEI Zhengfa, ZHANG Juncai, CAO Xiaoyan, et al. Causes and Influential Factor Analysis of Landslides and Rockfalls in North & South Mountain Areas of Xining City, Qinghai Province[J]. The Chinese Journal of Geological Hazard and Control, 2021, 32(4): 47-55.
[19] 代聪, 李为乐, 陆会燕, 等. 甘肃省舟曲县城周边活动滑坡InSAR探测[J]. 武汉大学学报(信息科学版), 2021, 46(7): 994-1002. DAI Cong, LI Weile, LU Huiyan, et al. Active Landslides Detection in Zhouqu County, Gansu Province Using InSAR Technology[J]. Geomatics and Information Science of Wuhan University, 2021, 46(7): 994-1002.
[20] DU Y, HUANG G W, ZHANG Q, et al. A New Asynchronous RTK Method to Mitigate Base Station Observation Outages[J]. Sensors, 2019, 19(15): 3376.
[21] 王铎, 黄观文, 杜源, 等. 顾及运动状态改正的GNSS滑坡监测基准站切换方法[J]. 测绘学报, 2022, 51(10): 2117-2124. WANG Duo, HUANG Guanwen, DU Yuan, et al. Switching Method of GNSS Landslide Monitoring Reference Station Considering the Correction of Motion State[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(10): 2117-2124.
[22] WANG D, HUANG G W, DU Y, et al. Stability Analysis of Reference Station and Compensation for Monitoring Stations in GNSS Landslide Monitoring[J]. Satellite Navigation, 2023, 4(1): 29.
[23] BOCK Y, MELGAR D, CROWELL B W. Real-Time Strong-Motion Broadband Displacements from Collocated GPS and Accelerometers[J]. The Bulletin of the Seismological Society of America, 2011, 101(6): 2904-2925.
[24] GENG J H, BOCK Y, MELGAR D, et al. A New Seismogeodetic Approach Applied to GPS and Accelerometer Observations of the 2012 Brawley Seismic Swarm: Implications for Earthquake Early Warning[J]. Geochemistry, Geophysics, Geosystems, 2013, 14(7): 2124-2142.
[25] 景策, 黄观文, 张勤, 等. 方差膨胀模型用于GNSS/加速度计融合变形监测[J]. 大地测量与地球动力学, 2023, 43(5): 491-497. JING Ce, HUANG Guanwen, ZHANG Qin, et al. GNSS/Accelerometer Fusion Deformation Monito⁃ring Based on Variance Inflation Model[J]. Journal of Geodesy and Geodynamics, 2023, 43(5): 491-497.
[26] JING C, HUANG G W, ZHANG Q, et al. GNSS/Accelerometer Adaptive Coupled Landslide Deformation Monitoring Technology[J]. Remote Sensing, 2022, 14(15): 3537.
[27] JING C, HUANG G W, LI X, et al. GNSS/Acce-lerometer Integrated Deformation Monitoring Algorithm Based on Sensors Adaptive Noise Modeling[J]. Measurement, 2023, 218: 113179.
[28] 匡翠林, 戴吾蛟. GPS监测高层建筑风致振动变形及小波应用[J]. 武汉大学学报(信息科学版), 2010, 35(9): 1024-1028. KUANG Cuilin, DAI Wujiao. Measurement of Wind-Induced Vibration of Tall Buildings Using GPS and Wavelet Application[J]. Geomatics and Information Science of Wuhan University, 2010, 35(9): 1024-1028.
[29] KIJEWSKI T, KAREEM A. Wavelet Transforms for System Identification in Civil Engineering[J]. Computer⁃Aided Civil and Infrastructure Enginee⁃ring, 2003, 18(5): 339-355.
[30] SHAO X Y, XU C. Earthquake-Induced Landslides Susceptibility Assessment: A Review of the State-of-the-Art[J]. Natural Hazards Research, 2022, 2(3): 172-182.
[31] NOWICKI JESSEE M A, HAMBURGER M W, FERRARA M R, et al. A Global Dataset and Model of Earthquake-Induced Landslide Fatalities[J]. Landslides, 2020, 17(6): 1363-1376.
[32] MEO M, ZUMPANO G, MENG X L, et al. Measurements of Dynamic Properties of a Medium Span Suspension Bridge by Using the Wavelet Transforms[J]. Mechanical Systems and Signal Processing, 2006, 20(5): 1112-1133.