An Algorithm of GPS Single-Epoch Kinematic Positioning Based on Doppler Velocimetry
-
摘要: 在传统的GPS单历元动态定位参数估计中,由于先验约束信息较少,模糊度求解和定位性能往往较低。据此,研究了多普勒测速在GPS单历元动态定位中的应用,提出了一种附有GPS多普勒测速信息约束的坐标更新方法,即利用移动载体先验坐标和多普勒测速信息预测当前历元坐标。考虑到该方法进行坐标更新的精度较高和传统单点定位坐标更新的稳健性较强等特点,在此基础上给出了相应的单历元整周模糊度参数求解策略,进一步提高了该方法的定位性能。试验结果表明,所提出的算法相对于传统的无速度信息坐标约束的GPS单历元动态定位算法,其模糊度浮点解精度、模糊度固定率和平均定位精度均有了进一步提高,尤其是在观测卫星数较少、卫星几何结构较差的情况下更为明显。Abstract: Due to the less prior constraint information, the ambiguity resolution and positioning performance are not very good in the GPS single-epoch kinematic positioning parameter estimation. Thus, a new method of GPS single-epoch kinematic positioning based on Doppler velocimetry is proposed in this paper. The coordinate initialization using the velocity information is firstly discussed, the current epoch's coordinate of the moving vehicle is predicted with the priori coordinate and velocity information. Compared with the conventional method of single point positioning (SPP) coordinate initialization, it is expected to obtain higher precision, but poorer robustness as well. Thus, the corresponding strategy for single-epoch ambiguity resolution is further introduced so as to improve the positioning performance of the method in this paper. Contrasting with the conventional algorithm of GPS single-epoch kinematic positioning, the experimental results show that, the method proposed can improve the precision of the float ambiguity, the final ambiguity fixed rate and the average positioning accuracy, especially in the case that, the number of GPS satellites is not sufficient or the geometric structure of the satellites is not very good.
-
Keywords:
- Doppler velocimetry /
- single-epoch /
- kinematic positioning /
- ambiguity resolution
-
据中国地震台网发布,2022-01-08T01:45(北京时间),中国青海省海北州门源地区发生了Mw 6.7地震,震源深度约为10 km;地震震中(101.26° E,37.77° N)距离门源县城54 km,距离西宁市136 km,距离2016-01-26门源Mw 5.9地震的震中不足40 km。图 1为两次门源地震的区域构造背景,其中断层数据参见文献[1-2],白色实线框为不同轨道的合成孔径雷达(synthetic aperture radar,SAR)影像的覆盖范围。2022年门源地震的孕震背景为青藏高原东北缘的昌马堡-古浪-海原构造带,历史地震资料显示该区域地质构造背景复杂且强烈,曾发生过两次8级以上和多次7级以上大地震,如1920年8.5级海原地震、1927年8.0级古浪地震(见图 1)。该构造带内部断裂纵横,其中冷龙岭断裂带由一组几乎平行的北西西向分布的左旋逆走滑性质的断裂组成,长约120 km,本次地震位于冷龙岭断裂带西北端与托莱山断裂的阶区[3-4]。自21世纪以来,本次地震震中100 km范围内发生了13次5级及以上的中强震(https://earthquake.usgs.gov/)。
距离2022年地震最近的一次地震是震源机制与其完全不同的2016年Mw 5.9门源逆冲型地震,震中处于冷龙岭断裂中部以北。2016年地震北侧是民乐-大马营断裂,以挤压逆冲为主、略带左旋走滑性质(见图 1)。文献[5]利用高分遥感影像与震前高分影像推测2016年地震的发震断层为冷龙岭北断裂;文献[6]通过对省级固定地震台站及部分流动地震台站记录到的波形资料进行联合精定位,推测该地震为北西向冷龙岭断裂与北西向民乐-大马营断裂之间的一条盲断层;文献[7]利用合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)技术精确反演了发震断层的几何参数及滑移分布,并指出冷龙岭断裂复杂的几何形态导致了该地震断层滑移集中在较小区域,因此造成了较高水平的应力降,加剧了局部应力积累,增加了区域地震灾害风险。
2016年和2022年两次相继发生的门源地震为研究青藏高原东北缘的断层系统提供了机会。亟需解决的关键问题一是2022年地震震源机制的快速解算,二是判断2022年地震是否由2016年地震触发。门源地区位于河西走廊,是“一带一路”的重要纽带,也是诸多重大国家工程的必经之路,如西气东输等工程。因此对2022年门源地震的发震机理进行研究,厘清该区域的断层构造活动特征,亦可为国家重大工程的建设与运维作出贡献。本着快速响应与精确解算的原则,且考虑到发震区域地处高原山区,常规的大地测量手段难以获取连续、全覆盖的地表形变场,本文首先采用了InSAR技术分别获取两次门源地震的同震地表形变场,该技术具有全天时、覆盖范围广等显著优点,经过30余年的探索和发展,已经成为研究同震形变和反演震源参数必不可少的大地测量技术[8-13];然后利用弹性半空间的位错模型确定了两次事件的断层几何参数,并基于分布式滑动模型反演确定了两次地震断层面上的滑动分布;最后解算了区域库仑应力场,由此探讨了2016年门源地震对2022年门源地震的触发关系。表 1为不同研究机构或学者获取的两次门源地震的震源参数,其中USGS(United States Geological Survey)表示美国地质勘探局,GCMT(global centroid moment tensor)表示全球矩心矩张量项目。
表 1 2016年和2022年门源地震震源参数Table 1. Source Parameters for the 2016 and 2022 Menyuan Earthquakes地震 研究来源 Mw 走向/(°) 倾角/(°) 滑动角/(°) 长度/km 宽度/km 最大滑动量/m 最大滑动深度/km 2016年门源地震 GCMTa 5.9 146 43 83 ― ― ― ― 335 47 96 USGSb 5.9 141 50 79 ― ― ― ― 337 41 103 文献[14] 5.9 134 43 68 24 20 0.45 9.5 文献[7] 5.9 127 45 73 ― ― 0.75 ― 本文研究c 5.9 137 46 79 4.4 5.1 0.9 ― 本文研究d 5.9 137 44 79 16 16 0.53 9.0 2022年门源地震 GCMT 6.7 104 82 1 ― ― ― ― 14 89 172 USGS 6.6 104 88 15 ― ― ― ― 13 75 178 AB段c ― 104 80 0 8 13.7 1.07 ― BC段c 6.6e 109 80 5 8 13.7 2.84 ― AB段d ― 104 80 0 10 16 2.5 5 BC段d 6.7e 109 80 5 20 16 3 5 注:a表示数据来源于https://www.globalcmt.org/; b表示数据来源于https://earthquake.usgs.gov/;c表示均匀滑动分布;d表示分布式滑动分布;e表示AB、BC两段断层矩张量之和对应的矩震级 1 InSAR数据处理与同震形变场
1.1 InSAR数据处理
本文获取了覆盖2016年和2022年两次门源地震区域的Sentinel-1A升降轨影像数据,用于InSAR分析的SAR影像参数见表 2。Sentinel-1A的SAR影像处理采用GAMMA软件,利用二轨法的InSAR技术获取地震同震地表形变场[15-17]。干涉处理过程中采用美国宇航局发布的30 m空间分辨率的SRTM(shuttle radar topography mission)消除地形相位和进行地理编码。为了提升信噪比(signal-to-noise ratio,SNR),首先采用了距离向:方位向为10∶2的多视以及滤波处理;然后采用最小费用流(minimum cost flow,MCF)算法进行相位解缠[18],进而获取视线向的地表形变量,利用通用型卫星雷达大气改正系统(generic atmospheric correction online service for InSAR,GACOS)来降低大气噪声对形变场精度的影响[19-22];最后分别获取2016年、2022年两次门源地震的同震地表形变场,如图 2所示。
表 2 2016年和2022年门源地震用于InSAR分析的SAR影像参数Table 2. Parameters of SAR Images for InSAR Analysis on the 2016 and 2022 Menyuan Earthquakes地震发生年份/年 轨道 飞行方向 获取时间 空间基线/m 时间基线/d 主影像 辅影像 2016 128 升轨 2016-01-13 2016-02-06 14.2 24 033 降轨 2016-01-18 2016-02-11 6.0 24 2022 026 升轨 2021-12-29 2022-01-10 -107.6 12 033 降轨 2021-12-29 2022-01-10 55.1 12 1.2 同震形变场
许多学者利用InSAR技术对2016年门源地震进行了研究[2,7,14]。文献[7]利用InSAR技术获取的升降轨最大视线向的形变为5.8 cm和6.8 cm;并通过反演认为本次地震可能发生在海原断裂的一个次级断裂上,且该断层在地下10 km处和海原断裂相连。文献[14]利用升降轨Sentinel-1数据获取了同震地表形变场,最大视线向形变均为7.0 cm;同时综合利用余震信息,确定了发震断层参数,认为发震断层呈铲状。本文得到的形变场与文献[14]类似,图 2(a)和2(b)分别为2016年门源地震同震升、降轨视线向形变场。升降轨同震形变场的形变信息连续,基本完整,表明2016年门源地震的活动断层并未破裂到地表,升、降轨影像同震形变场视线向最大隆升值为6.7 cm和7.0 cm。此外,两个轨道的同震形变场结构较为单一均匀,方向较为一致,均显示地表靠近SAR卫星,符合逆断层型地震的运动特征。图 2(c)和2(d)分别为2022年门源地震同震升、降轨视线向形变场。由于026升轨SAR影像的覆盖限制,图 2(c)只覆盖了本次地震的西部区域,其最大形变为67 cm。完整的降轨视线向同震形变场显示,本次地震造成的地表形变范围达30 km×20 km,其最大形变为78 cm。形变场中间分割区大致呈NWW-SEE方向,这与USGS、GCMT等机构发布的震源参数的走向基本一致,且该方向与冷龙岭断裂基本平行。
2022年地震升降轨同震形变场的上下两盘表现出了相反的形变态势,同时同一轨道的形变场上下两盘也表现出了相反的运动状态,这种现象表明本次地震引起的地表形变以水平移动为主,符合走滑断层型地震的运动特征。同震视线向形变场的断层附近出现了较为严重失相干现象,主要原因为发震断层破裂到了地表,形变梯度较大,且引起了形变相位的不连续,这一现象与光学地表解译结果相吻合。图 3为2022-01-08高分七号影像探测到的2022年门源地震地表破裂带。从图 3可以看出,高分七号地表破裂解译结果显示本次地震造成了大于20 km的地表破裂,其中发震区域的西部地形平缓,地表破裂较连续,出现了较明显的拉张裂隙、挤压鼓包等地表破裂现象(见图 3(b)中的(1)~(3)区域);而发震区域的东部因穿过高山、且受山上冰雪覆盖等因素影响,地表破裂呈现零散分布(见图 3(b)中的(4)~(8)区域)。另外,通过光学影像解译获得的地表破裂为后续2022年地震的反演,特别是其走向参数提供了重要参考。
2 断层几何参数及滑动分布反演
本文利用Okada提出的弹性半空间位错理论[23],联合升降轨同震形变场数据,分别反演了2016年和2022年两次门源地震的断层信息。为了抑制噪声和提高反演的效率,反演前对同震形变场进行了四叉树降采样[24]。反演借助PSOKINV软件包进行了两步法反演策略:(1)假设矩形断层面上为均匀滑动,采用非线性搜索算法获取断层的几何参数,如发震断层的位置、断层走向、断层倾角等;(2)将发震断层面离散为一定尺寸的子断层执行分布式滑动反演,获取子断层各自的运动参数[25-27]。2016年和2022年门源地震同震滑动分布结果如图 4所示。
2.1 2016年门源地震反演
2016年门源地震在进行上述两步法反演的第一步均匀滑动反演时,根据整理的资料信息将断层的长度变化和宽度变化均设置为3~10 km,倾角变化设置为30°~60°。均匀反演得到的断层参数见表 1。
为了获取断层面上的精细滑动分布,首先将发震断层沿着走向和倾向分别扩展为16 km,然后将延长后的发震断层面离散为1 km×1 km的子断层进行分布式滑动反演。由图 4(a)可以看出,2016年门源地震发震断层的最大滑动量为0.53 m,位于地下9 km深度处,滑动分布主要集中在地下4~12 km区域,表明地震活动并未破裂到地表。子断层的滑移方向表明发震断层是一个具有少量左旋走滑分量的逆断层。设置本地区的剪切模量为30 GPa,则反演得出的矩震级为5.9。本文反演结果与文献[2,14]公布的数据吻合较好。
图 5为2016年门源地震同震地表形变场、模型拟合形变场及残差,其中左列为地表形变场观测值,中间列为均匀反演或分布式滑动反演结果模拟的同震形变场,右列为拟合残差,黑色矩形为投影到地面的断层面,而蓝色实线为断层面延长与地面的交线。由于2016年地震的地表形变较小,故成图时采用了解缠后的结果。由图 5可以看出,InSAR观测的升降轨形变场与相应的拟合结果基本一致,验证了本次地震断层反演的可靠性。
2.2 2022年门源地震反演
虽然026升轨影像只提供了2022年门源地震的部分地表形变场,但是升降轨SAR影像对走滑断层破裂的地表形变成像具有互补的特点,因此,本文联合了升、降轨地表形变场进行反演。在反演时,由于InSAR解译的形变场已经破裂到地表,故将断层的上边界设置为0 m。光学解译的地表破裂结果描绘出一条NWW-SEE向的断层迹线(见图 3(a)),结合此迹线走向和InSAR同震形变场的分布将反演模型设置为AB和BC两个断层,其中AB段的走向固定为104°,BC段的走向固定为109°。长度变化设置为8~40 km,宽度变化设置为6~20 km,断层倾角变化设置为78°~90°。均匀反演得到的断层参数见表 1。
在获得发震断层的几何参数后,将AB段发震断层沿着走向延长至10 km,沿着倾向延伸至16 km,将BC段发震断层沿着走向延长至20 km,沿着倾向延伸至16 km,并将延长后的发震断层离散成1 km×1 km的子断层进行分布式滑动反演。图 4(b)为分布式反演得到的断层滑动分布,子断层的运动方向表明该断层主要受左旋滑动控制,属于高倾角走滑型地震。子断层最大滑动量为3.5 m,出现在地下约4 km处,滑动主要集中在地下2~7 km处。反演获得矩震级为6.7,与USGS、GCMT等机构公布的发震断层数据吻合较好。
图 6为2022年门源地震同震地表形变场、模型拟合形变场及残差。由图 6可以看出,InSAR观测的升降轨形变场与相应的拟合结果基本一致。值得指出的是,最右列的残差在断层近场的两侧较为明显,主要原因包括:(1)近场大梯度形变引起的低相干乃至失相干;(2)低相干引起的解缠误差;(3)活动断层几何结构的简化。
综合研究结果,可以初步判定2022年门源地震的发震断层主体为冷龙岭西段西北端,且本次地震活动极有可能向西破裂到了其西侧的托莱山断裂。理由如下:(1)基于InSAR同震地表形变场反演确定的震源机制与冷龙岭断裂的运动性质吻合;(2)光学解译结果(图 1)显示,本次地震引起的地表破裂在冷龙岭断裂西端稍向南偏转,延伸至托莱山断裂(AB段),而且通过滑移分布建模显示该断层表现为左旋走滑类型,与托莱山断裂的运动性质吻合。
3 区域库仑应力变化对2022地震诱发关系探讨
强震发生后,发震区域的应力状态将会发生改变,进而延缓或者促进邻近区域活动断层上地震的发生,大量的地震触发关系研究表明地震之间存在相互作用[28-30]。本文运用Coulomb 3.3软件包进行库仑应力分析,在处理过程中泊松比设置为0.25,摩擦系数设置为0.4[31-33]。以2016年门源地震断层滑动分布作为源断层,本文反演得到的2022年门源地震的断层几何作为接收断层,分别采用5 km和10 km的深度计算2016年门源地震引起的库仑应力变化对2022年门源地震地震的影响。
2022年门源地震发生在2016年门源地震震后的第6年,且两次地震震中距离不足40 km。图 7为2016年门源地震引起的邻近区域静态库仑应力变化。由图 7可以看出,2016年地震断层滑移在5 km和10 km的深度对2022年发震断层面上均产生了正库仑应力变化,且2022年地震发震断层整体处于库仑应力变化的增强区,这说明2016年地震对2022年地震的发生具有一定的促进作用。类似的库仑应力触发地震事件在巴颜喀拉块体东端的地震序列、2017年Iran地震、2020年Mexico地震等震例中均有研究[34-36]。另外,文献[37]发现,发生在冷龙岭断裂北部的2016年门源地震的逆冲运动为冷龙岭断裂上地震活动腾出了空间。
因为2022年门源地震所处的构造背景中分布着较多左旋走滑性质的断裂,因此本文计算的库仑应力结果对分析发震区域整体的应力分布也具有参考意义。从图 7分析得出,2016年门源地震引起的库仑应力的增强区主要分布在其西北部的Ⅰ、Ⅱ区域和东南部的Ⅲ区域,其中,冷龙岭北断裂、肃南-祁连断裂、皇城-双塔断裂西端、民乐-大马营断裂东南端、托莱山断裂、冷龙岭断裂东端以及门源断裂东端均存在较为明显的库仑应力增强(大于0.1 bar),地震发生风险增高,灾害风险加大,因此有必要对这些地区开展持续监测和详细的风险评估。
4 结语
本文对2022年门源地震做了快速响应,利用Sentinel-1A卫星升、降轨影像数据分别获取了2016年和2022年两次门源地震高精度同震地表形变场,经反演确定了两次门源地震发震断层的几何参数和滑动分布模型,探讨了2016年地震事件对2022年门源地震的影响,得出以下结论:
1)2016年门源地震升、降轨最大视线向隆升值分别为6.7 cm和7.0 cm,且同震地表形变场方向较为一致,符合逆断层型地震的运动特征。2022年门源地震同震地表形变场显示,地表破裂沿NWW-SEE方向延伸,地震造成的地表形变范围达30×20 km,其最大形变为78 cm。
2)2016年门源地震的滑动分布显示,发震断层的最大滑动量为0.53 m,断层滑动主要集中在地下4~12 km的区间。根据2022年门源地震的滑动分布,子断层最大滑动量为3.5 m,出现在地下约4 km处。初步判定2022年地震发震断层主要为冷龙岭断裂的西段,且极有可能破裂到了其西北端西侧的托莱山断裂。
3)以2022年门源地震的发震断层为接收断层,计算了2016年门源地震造成的库仑应力变化。结果显示2016年门源地震在2022年地震断层面上产生了明显的正库仑应力变化,在一定程度上触发了2022年地震的发生。同时,2016年门源地震引起了冷龙岭北断裂、肃南-祁连断裂、皇城-双塔断裂西端、民乐-大马营断裂东南端、托莱山断裂、冷龙岭断裂东端以及门源断裂东端的库仑应力增加,这些断裂的地震危险性在今后需要得到持续关注。
本文利用近实时数据对2022年地震进行快速响应,对发震断层进行了初步建模,为救灾救援,灾害评估,应急响应提供第一手的数据支撑。后期将收集更全面的数据,利用多种技术进行处理,尤其是利用像素偏移量跟踪法获取近场大形变数据,对断层几何模型,滑移分布及震源机制作出更精细的建模和分析,评估研究区域的中长期地震灾害风险。
-
表 1 3种方法模糊度解算结果对比
Table 1 Comparison of Ambiguity Resolution Results Computed by the Three Methods
方法 平均
Ratio值平均
ADOP值模糊度
固定率/%方法1 3.6 0.24 76.5 方法2 6.1 0.14 87.7 方法3 13.1 0.09 91.3 表 2 3种方法平均定位精度对比/m
Table 2 Comparison of Average Positioning Accuracies Computed by the Three Methods/m
方法 dX dY dZ 方法1 0.042 3 0.052 1 0.047 6 方法2 0.025 2 0.032 0 0.023 4 方法3 0.014 9 0.020 6 0.013 4 -
[1] 陈永奇.单历元GPS变形监测数据处理方法的研究[J].武汉测绘科技大学学报, 1998, 4(4):324-328 http://www.cqvip.com/Main/Detail.aspx?id=3290759 Chen Yongqi. Development of the Methodology for Single Epoch GPS Deformation Monitoring[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1998, 4(4):324-328 http://www.cqvip.com/Main/Detail.aspx?id=3290759
[2] Teunissen P J G, Odolinski R, Odijk D. Instanta-neous BeiDou+GPS RTK Positioning with High Cut-off Elevation Angles[J]. Journal of Geodesy, 2014, 88(4):335-350 doi: 10.1007/s00190-013-0686-4
[3] Parkins A. Increasing GNSS RTK Availability with a New Single-Epoch Batch Partial Ambiguity Resolution Algorithm[J]. GPS Solutions, 2011, 15(4):391-402 doi: 10.1007/s10291-010-0198-0
[4] Wu Z, Bian S, Ji B, et al. Short Baseline GPS Multi-frequency Single-Epoch Precise Positioning:Utilizing a New Carrier Phase Combination Method[J]. GPS Solutions, 2016, 20(3):373-384 doi: 10.1007/s10291-015-0447-3
[5] Wang L, Feng Y, Guo J. Reliability Control of Single-Epoch RTK Ambiguity Resolution[J]. GPS Solutions, 2016, 21(2):1-14 doi: 10.1007/s10291-016-0550-0
[6] 陈远, 于兴旺, 叶聪云, 等. GPS多普勒观测值测速的精度分析[J].全球定位系统, 2008, 33(1):31-34 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chtb201303003 Chen Yuan, Yu Xingwang, Ye Congyun, et al. Analysis of Single-Point GPS Velocity Determination with Doppler[J]. GNSS World of China, 2008, 33(1):31-34 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chtb201303003
[7] 肖云, 孙中苗, 程广义.利用GPS多普勒观测值精确确定运动载体的速度[J].武汉测绘科技大学学报, 2000, 25(2):113-118 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=whchkjdxxb200002005 Xiao Yun, Sun Zhongmiao, Cheng Guangyi. Precise Determination of Velocity for Airborne Gravimetry Using the GPS Doppler Observations[J]. Journal of Wuhan Technical University of Surveying and Mapping, 2000, 25(2):113-118 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=whchkjdxxb200002005
[8] Kubo N, Pullen S. Instantaneous RTK Positioning Based on User Velocity Measurements[C]. The 21st International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2008), Savannah, GA, 2008
[9] Kubo N. Advantage of Velocity Measurements on Instantaneous RTK Positioning[J]. GPS Solutions, 2009, 13(4):271-280 doi: 10.1007/s10291-009-0120-9
[10] 杨元喜.自适应动态导航定位[M].北京:测绘出版社, 2006 Yang Yuanxi. Adaptive Navigation and Kinematic Positioning[M].Beijing:Surveying and Mapping Press, 2006
[11] Zhao S, Cui X, Guan F, et al. A Kalman Filter-Based Short Baseline RTK Algorithm for Single-Frequency Combination of GPS and BDS[J]. Sensors, 2014, 14(8):15415-15433 doi: 10.3390/s140815415
[12] Dai F, Mao X. BDS/GPS Dual Systems Positioning Based on Kalman Filter in Urban Canyon Environments[C]. The 17th International IEEE Conference on Intelligent Transportation Systems, Qingdao, China, 2014
[13] Teunissen P J G. The Least-Squares Ambiguity Decorrelation Adjustment:A Method for Fast GPS Integer Ambiguity Estimation[J]. Journal of Geodesy, 1995, 70(1-2):65-82 doi: 10.1007/BF00863419
[14] Teunissen P J G. An Optimality Property of the Integer Least-Squares Estimator[J]. Journal of Geodesy, 1999, 73(11):587-593 doi: 10.1007/s001900050269
[15] Teunissen P J G, Dennis O. Ambiguity Dilution of Precision: Definition, Properties and Application[C]. ION GPS 1997, Kansas City, MO, 1997
-
期刊类型引用(2)
1. 赵彬彬,汤鑫,唐忠立,王安. 基于交集成分细分的面/面目标拓扑关系扩展模型. 长沙理工大学学报(自然科学版). 2021(02): 66-75 . 百度学术
2. 杜晓初,蔡贤. 栅格空间中不确定空间拓扑关系的判别. 测绘科学. 2020(10): 119-126 . 百度学术
其他类型引用(2)