image
Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 
Online first  /  Current issue  /  All issues  /  Top access  /  RSS

Progress in Geophysics >

2025 , Vol. 40 >Issue 4: 1349 - 1360

DOI: https://doi.org/10.6038/pg2025HH0118

Tonga volcanic in 2022 eruption monitoring based on spherical harmonic expansion precision single point positioning: the case of TONG at IGS station

  • MingHao CHANG , 1 ,
  • JinYun GUO , 1, * ,
  • HengYang GUO 1, 2 ,
  • XiaoTao CHANG 3 ,
  • YanChen CAO 1
Expand
  • 1 College of Geodesy and Geomatics, Shandong University of Science and Technology, Qingdao 266590, China
  • 2 School of Land Science and Technology, China University of Geosciences (Beijing), Beijing 100083, China
  • 3 Land Satellite Remote Sensing Application Cente, MNR, Beijing 100048, China

Received date: 2023-06-05

  Online published: 2025-08-28

Copyright

Copyright ©2025 Progress in Geophysics. All rights reserved.
Fold

Abstract

In January 2022, Tonga experienced a massive underwater volcanic eruption, causing significant impacts throughout the Pacific region. In this article, a precise positioning technique based on spherical harmonic expansion is used to decompose the errors along the propagation path of GPS signals. This technique allows for precise single-point positioning calculations of GPS data received at the TONG station within the Kingdom of Tonga. The results of these calculations are then utilized to study the response to the Tonga volcanic eruption.The TONG station's standard deviations in the East (E), North (N), and Up (U) directions are 1.32 cm, 0.63 cm, and 1.53 cm, respectively. During the volcanic eruption, the positional error time series data exhibits a noticeable response. In the computed coordinate time series data, there is a clear shift in the station's position after the volcanic eruption. By performing a fitting process, it is determined that the station moved in a southeast direction, with an eastward velocity of 0.0108 m/d, a southward velocity of 0.0023 m/d, and a downward velocity of 0.0149 m/d. By employing the Singular Spectrum Analysis (SSA) method to decompose the time series data, the principal components and residuals are obtained for each direction (E, N, U). In the E direction, as well as the U direction, there are noticeable periodic variations. The amplitude of the E direction is 0.045 cm, while the N direction has an amplitude of 0.05 cm. All three directions exhibit significant fluctuations around the time of the volcanic eruption, particularly around January 15th.The residuals in the N and U directions exhibit significant changes at the moment of the volcanic eruption. In particular, the U direction's residual decreased by 0.4 cm after the eruption compared to the pre-eruption stability. Additionally, the error correction of the propagation path, based on spherical harmonic expansion, shows a clear response to the errors caused by volcanic ash. The magnitude of these errors suggests that the volcanic ash has shifted in the southeast direction. Additionally, by analyzing the distribution of errors at different elevation angles, it is possible to observe the spread and ascent of volcanic ash. This demonstrates that this method can be employed for studying volcanic eruptions.

Key words： Tonga volcano; Precise point positioning; Spherical harmonic expansion; Volcanic ash; GPS

Cite this article

MingHao CHANG , JinYun GUO , HengYang GUO , XiaoTao CHANG , YanChen CAO . Tonga volcanic in 2022 eruption monitoring based on spherical harmonic expansion precision single point positioning: the case of TONG at IGS station[J]. Progress in Geophysics, 2025 , 40(4) : 1349 -1360 . DOI: 10.6038/pg2025HH0118

0 引言

2022年1月汤加海底火山发生了近90年最大的爆发，此次火山喷发对自然环境造成了极大的影响.火山喷发造成的海啸影响到了整个太平洋地区，环太平洋的新西兰、智利以及日本等国家监测到明显的海啸波.同时火山喷发所喷发出的大量火山灰物质，对大气影响非常严重，火山爆发之后火山灰覆盖汤加岛屿，对汤加的经济生活以及太平洋地区的航空业造成了严重的影响.火山灰主要是由火山喷发出的火山碎屑物和气体等物质构成，火山碎屑物主要成分是直径小于2 mm的岩石和矿物等；气体主要包括了水汽、CO2、SO2等(Flynn et al., 2001).火山气象台监测火山的喷发主要通过地基以及空基的设备进行观测以及预警，随着技术发展，例如激光雷达和遥感卫星等技术也逐步用于火山喷发的监测.火山灰的移动也会影响到全球定位系统(Global Positioning System，GPS)的卫星信号，因此使用GPS监测火山喷发存在很大前景.
火山喷发后，火山灰进入大气会对卫星信号的传播在大气的传播产生干扰.Houlié等(2005)对2000年8月爆发的日本三宅岛火山进行研究，使用了地震层析成像算法，根据无电离层延迟测量的异常，将异常与火山灰温度都映射在时空图中，确认了热火山灰对于GPS信号造成了延迟.GPS掩星技术利用无线电信号穿过大气层时的弯曲信息对大气信息进行反演(Gulbis et al., 2006)，用来探测地球大气以及电离层(Jin et al., 2014).随着GPS卫星技术逐渐成熟，使用GPS进行大气探测得到发展(Yunck et al., 2000).1995年，美国UCAR(University Corporation for Atmospheric Research)主导的地球大气掩星探测项目GPS/MET搭载卫星升空并成功进行了大气探测，论证了使用GPS掩星技术探测地球大气层的可行性(Ware et al., 1996).可以根据反演得到的大气电离层电子密度以及对流层的大气参数，如温度、密度等研究火山喷发后进入大气的火山灰的各种信息.
GPS信噪比(Signal-to-Noise Ratio，SNR)在火山喷发的监测中也有很大的前景(Ohta and Iguchi, 2015).卫星信号通过火山灰时会造成衰减，而SNR数据是接收信号功率与噪声功率的比值，对水汽并不敏感但是会受到火山灰的影响.Larson (2013)对2008年的阿拉斯加州奥克莫克火山以及2009年的力道特火山进行研究，使用SNR数据进行火山喷发的研究，结果表明在火山喷发初期的时候信噪比会有一定程度的下降.在2012年的新西兰汤加里罗火山喷发中，Fournier和Jolly(2014)使用了信噪比探测方法以及相位残差法进行火山探测，证实了GPS数据在探测火山喷发的可行性.Grapenthin等(2013)使用GPS信号的相位残差，对卫星信号穿过火山灰与监测站之间的信号延迟联系，提出利用相位残差来检测火山的方法.但上述方法依然存在着不少问题，GPS信噪比方法虽然不受水汽的影响，但是受到电离层的影响，并不是任何的SNR监测都能够检测到火山灰(Aranzulla et al., 2015).非差拟合相位残差的方法则是容易受到大气中水汽以及温度的影响，并且由于火山灰的高度不一定，难以适用于所有火山爆发情况(Ohta and Iguchi, 2015).
精密单点定位方法(Precise Point Positioning，PPP)能够精密定位(张小红等，2017)，并且对于卫星信号在传播中的误差进行了精密的修正，因此可以使用大气传播中修正参数对传播路径上的大气进行研究，也可以得到火山喷发后进入大气的物质信息.在卫星信号传播的误差中，大体分为电离层误差、对流层误差(Liu et al., 2016；)以及多路径效应(Ning and Yuan, 2017；蒋光伟等, 2018).
在对流层延迟上，PPP定位技术无法使用双差的方法来消除影响，所以一般使用参数估计法，经验模型法以及外部修正等方法来进行改正(闫建巧等，2016).较为广泛的有Saastamoinen模型(Saastamoinen，1972)，它考虑到了对流层模型的高度问题，将温度这一参数改变成为常数并且按照高度进行了划分，这样分层进行计算，模型将在12 km以下的对流层使用梯度函数进行计算，按照每千米6.5 ℃进行计算，在更高层次的对流层则是设置为一个常数.相比较于单一温度设置的Hopfield的模型(Hopfield，1969)，能够更加符合现实，计算精准.UNB3模型，则是利用气象数据资料得到一个大气的参数表.利用之前的测量数据以及资料，不需要当时的实测数据，利用计算好的格网数据可以内插得出需要的气象数据.
在电离层延迟方面，单频的精密单点定位技术无法组成消电离层组合，并且电离层很受到地域性的影响，难以直接消除(李征航和张小红，2009)，多使用模型进行改正，常用的有克罗布歇模型(Klobuchar，1987)，格网模型还有历元差分模型等等, 也有人利用电离层测高仪来估计(赵聪等，2022).克罗布歇模型利用振幅、周期和所处的纬度以及根据日期以及辐射常量来设定的常数来计算得到接收信号的天顶方向的电离层延迟.GIM模型是很常见的电离层模型(Feltens，2003)，它利用IGS网站提供的全球电离层格网模型文件来得到信号传入点的具体位置，同时利用网格数据内插得到电子含量以此来求得电离层延迟.SEID模型则是根据观测值的历元的变化信息，对每颗卫星都构建电离层的延迟模型(Deng et al., 2009).双频的精密单点定位可以使用双频消电离层组合的方法来消除信号传播中的电离层延迟一阶项的误差，非常简单方便快速，是非常简单以及普遍的方法(殷志祥和孙明博，2018).
传统的误差延迟改正方法在信号传播过程中都存在各种局限性，不能普遍适用.本文使用一种基于球谐展开的精密定位方法，利用球谐展开的方法来表示与测站和卫星之间高度角以及方位角的误差项，修正信号传播过程中带来的误差，以便能够快速准确地进行精密单点定位，对汤加火山喷发进行研究.

1 研究区与数据

1.1 研究区域

洪阿哈阿帕伊岛海底火山(175.38°W，20.57°S) 位置如图 1所示.火山位于南太平洋岛国汤加王国境内，绝大部分位于海平面以下，这座岛屿本身是一座火山，最高海拔149 m.该火山位于汤加的福努阿福欧岛的东南方向大概30 km处，位于汤加王国最大的岛屿汤加塔布岛的西北方向65 km处.洪阿哈阿帕伊岛海底火山是是汤加-克尔马德克群岛火山弧的一部分，处于太平洋板块与印度洋板块的交界处.由于这里是太平洋板块沿着汤加海沟俯冲进地幔深处之后形成的俯冲带，造成这里浅源地震、火山活动频发.北京时间的2022年1月14日、15日，洪阿哈阿帕伊岛海底火山发生了猛烈喷发.
图1 汤加站点与火山位置图

Fig 1 Location map of Tonga station and volcano

1.2 数据

根据汤加火山喷发时间，本文使用2021年12月到2022年1月两个月的汤加站点TONG的GPS数据.本文实验使用的数据均来自IGS综合分析中心(Analysis Center Coordinator，ACC)，通过Crustal Dynamics Data Information System(CDDIS)下载，包括采样间隔为30 s的GPS观测文件(https://cddis.nasa.gov/archive/gnss/data/daily/)；导航电文(https://cddis.nasa.gov/archive/gnss/data/daily/)；混合精密星历(https://cddis.nasa.gov/archive/gnss/products/mgex/)以及钟差数据文件(https://cddis.nasa.gov/archive/gnss/products/).观测文件使用RINEX2版本，混合观测数据文件.精密星历文件使用IGS组织提供的最终精密星历，星历使用了IGb14的参考框架(刘路等，2021).
IGS汤加站点Tonga(175.179°W，21.145°S)，海拔为56.3 m，位于汤加王国的努库阿洛法(Nuku'alofa)，地处于汤加塔布岛的北海岸.该站点位于火山的东南方向，相距约67.26 km，此站点可以接收GPS卫星信号数据.

2 原理方法

本文采用了基于球谐展开的精密单点定位方法，利用球谐展开式表示了卫星与测站之间路径的相关误差，进行定位得到坐标.基于球谐展开的精密单点定位方法，既能够在定位精度上保证在厘米级别，同时又可以利用球谐展开式对火山喷发情况做出响应.
传统的精密单点定位利用载波相位观测值(Zumberge et al., 1997)，加入传播过程中的误差改正项以及其余误差项(Carlin et al., 2021)，建立观测方程为：
$\begin{align*}\varphi_{i}^{j}= & R_{i}^{j}+c \cdot\left\lfloor\left(t_{\mathrm{r}}\right)_{i}-\left(t_{\mathrm{s}}\right)_{i}^{j}\right\rfloor+\left(\rho_{\text {trop }}\right)_{i}^{j}+ \\& \left(\rho_{\text {ion }}\right)_{i}^{j}+\lambda N_{i}^{j}+\varepsilon_{\mathtt{φ}}, \end{align*}$
式中，$\varphi_{i}^{j}$表示卫星j在第i观测历元的载波相位观测值；$R_{i}^{j}$表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；$t_{\mathrm{r}}$表示接收机钟差，$t_{\mathrm{s}}$ 表示卫星钟差；($\rho_{\text {trop }}$)表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；($\rho_{\text {ion }}$)表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；N表示了载波相位观测值的整周模糊度参数，使用了Lambda算法作为模糊度固定策略；c表示光速；λ表示载波的波长；$\varepsilon_{\mathtt{φ}}$表示为观测数据残差.
信号传播路径误差同卫星与测站两者之间的方位角和高度角有关，因此将公式(1)中的关于传播路径的误差项使用球谐展开式表达，可得到基于球谐展开的精密单点定位观测方程：
$\begin{align*}\varphi_{i}^{j}= & R_{i}^{j}+c \cdot\left[\left(t_{\mathrm{r}}\right)_{i}-\left(t_{\mathrm{s}}\right)_{i}^{j}\right]+ \\& \sum\limits_{n=0}^{N_{\max }} \sum\limits_{m=0}^{n} P_{n m}\left(\cos \left(e_{i}^{j}\right)\right)\left[C_{n m} \cos \left(m\left(\alpha_{i}^{j}\right)\right)+\right. \\& \left.S_{n m} \sin \left(m\left(\alpha_{i}^{j}\right)\right)\right]+\lambda N_{i}^{j}+\varepsilon_{\mathtt{φ}}, \end{align*}$
式中，n为球谐展开的阶数，m为球谐展开的次数，为球谐展开的最大阶数；$P_{n m}\left(\cos \left(e_{i}^{j}\right)\right)$ 表示n阶m次的缔合勒让德多项式；Cnm和Snm表示n阶m次的球谐函数模型的系数，Cnm和Snm为球谐展开的待求参数；$e_{i}^{j}$和$\alpha_{i}^{j}$ 分别表示测站与j卫星在第i观测历元之间的高度角和方位角(山东科技大学，2022).由公式(2)的球谐展开精密单点定位观测方程可以得到此方法的误差方程：
$\begin{align*}\left(\nu_{\varphi}\right)_{i}^{j}= & R_{i}^{j}+c \cdot\left[\left(t_{\mathrm{r}}\right)_{i}-\left(t_{\mathrm{s}}\right)_{i}^{j}\right]+ \\& \sum\limits_{n=0}^{N_{\max }} \sum\limits_{m=0}^{n} P_{n m}\left(\cos \left(e_{i}^{j}\right)\right)\left[C_{n m} \cos \left(m\left(\alpha_{i}^{j}\right)\right)+\right. \\& \left.S_{n m} \sin \left(m\left(\alpha_{i}^{j}\right)\right)\right]+\lambda N_{i}^{j}-\varphi_{i}^{j}, \end{align*}$
式中$\left(\nu_{\varphi}\right)_{i}^{j}$ 表示卫星j在第i观测历元的载波相位观测数据的改正值.对公式(3)进行线性化处理，在测站近似坐标(X0, Y0, Z0)处进行泰勒级数展开并保留一阶项，可得到此方法的精密单点定位方程的线性化表达式：
$\begin{align*}\left(v_{\mathtt{φ}}\right)_{i}^{j}= & I_{i}^{j} \mathrm{~d} X+J_{i}^{j} \mathrm{~d} Y+K_{i}^{j} \mathrm{~d} Z+d\left(t_{r}\right)_{i}+\mathrm{d} N_{i}^{j}+ \\& \left(A_{n m}\right)_{i}^{j} \mathrm{~d} C_{n m}+\left(B_{n m}\right)_{i}^{j} \mathrm{~d} S_{n m}-\left(\omega_{\varphi}\right)_{i}^{j}, \end{align*}$
式中，$I_{i}^{j}$ 是测站近似坐标与精密星历中卫星j在第i观测历元中的坐标计算得到的测站坐标X0的改正数dX的系数；$J_{i}^{j}$ 是测站近似坐标与精密星历中卫星j在第i观测历元中的坐标计算得到的测站坐标Y0的改正数dY的系数；$K_{i}^{j}$是测站近似坐标与精密星历中卫星j在第i观测历元中的坐标计算得到的测站坐标Z0的改正数dZ的系数.$\mathrm{d} N_{i}^{j}$表示整周模糊度参数的改正数；$d\left(t_{r}\right)_{i}$表示在第i观测历元测站接收机钟差的待求参数，因为接收机钟差没有办法使用多项式拟合的很好(东旭华等，2021)，因此将其也作为待求参数，并且系数也设置为1.Anm表示代表了球谐展开式中参数Cnm的系数，Bnm表示代表了球谐展开式中参数Snm的系数.i表示观测历元数，i=1, 2, …, epoch，epoch为最大的观测历元数.
由上述误差公式可以看出，观测方程有epoch×j个观测方程，观测方程数目根据观测历元以及每个观测历元能够观测到的卫星数目来决定.其中的待求参数有3个位置参数：dX, dY, dZ；epoch(观测历元数目)个接收机钟差$\mathrm{d} t_{\mathrm{r}}$ 以及(Nmax+1)2个球谐展开式的参数Cnm、Snm.当选择的观测历元数目足够多，改正的误差项可以线性化描述，就可以联立方程组通过最小二乘来解出待求参数.通过设定滑动窗口来选择的观测历元数目i，再根据每个历元观测的j颗卫星数目，可以组成i×j个方程.
根据上述得到的球谐展开的精密单点定位误差方程，得到此误差方程的系数矩阵，改正数矩阵，未知参数矩阵以及常数项矩阵分别为矩阵B、V、X、L.根据平差原理，将误差方程的线性化表达利用矩阵形式表达(Dai，2002)，表达式为：
$\boldsymbol{V}=\boldsymbol{B} \boldsymbol{X}-\boldsymbol{L}, $
则未知参数X的最小二乘估计为：
$\boldsymbol{X}=\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{L}, $
因为使用了球谐展开式来代替信号传播的路径误差，所以给解算方程增加了球谐的相关未知数，因此单个历元中的观测卫星数量不足以解算求得方程中所有的未知数.在解算时，需要多个历元的数据来进行联立方程组解算坐标；这样就需要一个历元窗口来滑动选择合适的历元长度来解算坐标位置.
使用序贯的方法，将所有观测历元的观测值构成一个数据集，通过滑动历元窗口来选择合适的历元数目进行方程组的构建以及解算.在解算完成一组方程后，将历元窗口向后进行滑动，将原来观测数据矩阵中的第一个观测历元的数据推走，按顺序推进下一观测历元的观测数据.同时在平差解算过程中的B、V、L矩阵中也做同样的处理，将矩阵中的前面第一行数据剔除掉，将整体向下推进，最后一行读取新的观测历元的矩阵数值，这样得到了新的观测历元矩阵，可以继续按照原来方法进行方程组的解算得到坐标结果值.
在最后的结果输出上，将得到的地心地固直角坐标系坐标下的结果，转换成为站心坐标系下的结果.因为站心坐标系是一个更方便判断误差的局部坐标系，只要站心点的位置选择合理，表达的地理坐标都会是很小的值，非常便于计算.将选择的站点真值作为站点坐标系的站点位置坐标，解算得到的大地坐标依此转换为站心坐标系，得到的结果就是E、N、U方向的坐标，可以相较于大地坐标更加直接的看出各个方向的误差以及偏移.

3 结果和分析

3.1 解算结果

球谐展开的阶数不同也会导致方程未知数的不同，球谐展开式的阶数n同球谐展开式内的未知数m的关系为m=(n+1)2－1.高阶次的球谐展开会带来更多的未知数，因此需要更多的观测历元来参与解算，并且导致方程解算时间也会增长.
为了选择更好的解算阶数，首先对基于球谐展开式的精密单点定位方法的各个阶数进行实验.选择汤加站点的2021年4月10日的GPS数据进行实验性解算，以2021年4月10日IGS公布的天解坐标为真值，得到在相同的15个观测历元下不同球谐阶数的解算结果，表 1为解算对比结果.
表1 5~8阶球谐展开TONG站解算结果

Table 1 The results of TONG solved by PPP when spherical harmonic expanding to degree 5~8

阶数 Max/m Min/m Mean/m STD/m RMS/m
5阶 E -0.0469 -0.0035 -0.0234 0.0070 0.0244
N -0.0142 0.0152 -0.0002 0.0060 0.0059
U 0.0411 0.0800 0.0776 0.0061 0.0579
6阶 E -0.0424 -0.0077 -0.0242 0.0057 0.0249
N -0.0090 0.0119 0.0004 0.0044 0.0045
U 0.0503 0.0720 0.0621 0.0043 0.0623
7阶 E -0.0400 -0.0137 -0.0248 0.0048 0.0253
N -0.0074 0.0099 0.0008 0.0036 0.0037
U 0.0565 0.0708 0.0646 0.0033 0.0647
8阶 E -0.0354 -0.0221 -0.0277 0.0032 0.0278
N -0.0047 0.0077 0.0008 0.0030 0.0031
U 0.0612 0.0678 0.0643 0.0018 0.0644
根据表 1，基于5阶球谐展开的定位结果精度在E方向上为7 cm，N和U方向都为6 cm左右；基于6阶球谐展开的定位精度在N和U方向能达到4 cm；基于7阶和8阶球谐展开的定位精度在N和U方向能达到3 cm.可以看出随着球谐展开式阶数的增加，定位精度也在逐渐提高，标准差也在逐渐减小，说明数据更为精确；均方根误差并未随着阶数的增加而减少，说明更高阶的算法可能会造成特大或特小的误差使精度降低.同时，更高阶数的解算也会花费更长时间.其中得到5阶球谐展开的定位结果需要149 s，6阶需要213 s，7阶需要914 s，8阶需要一天的时间.因此在综合考虑各个阶次解算的精度以及所花费的时间之后，选择使用6阶球谐展开式15观测历元的解算策略，既能够得到较高的定位精度也能花费较少的时间.
转换为站心坐标系时使用IGS公布的周解坐标产品作为真值，利用2021年2186到2190的5个GPS周的周解坐标产品作为坐标的精确坐标.利用内华达大地测量实验室的汤加站点数据(http://geodesy.unr.edu/NGLStationPages/GlobalStationList)作为对比数据(Blewitt et al., 2018)，内华达实验室时序数据框架为IGS14，频率为5 min，公布的站点点位坐标为纬度：21.145°S，经度：175.179°W，高度：56.310 m.分别统计解算的站点坐标各个方向的最大值(MAX)、最小值(MIN)、平均值(MEAN)、标准差(STD)以及均方根(RMS).
表 2为利用球谐展开式计算精密单点定位方法解算的汤加站点误差统计结果.由表 1可知，最大的误差在U方向，达到分米级别，各个方向的平均值能够达到厘米级别，标准差精度也达到厘米级别.E方向的定位精度在1.32 cm左右，N方向的定位精度为0.6 cm左右，U方向的定位精度也有1.5 cm.内华达大地测量实验室的站点数据，各个方向上的定位精度都在厘米左右.对比可以看到使用球谐展开方法的精密定位技术和实验室数据的精度相等，使用球谐展开方法的精密定位技术的各个方向上误差最大最小值略好于实验室数据.
表2 TONG站点定位结果对比

Table 2 The results of spherical harmonic expansion single point positioning of TONG

MAX MIN MEAN STD RMS
球谐展开单点定位结果/m E 0.0358 -0.0549 -0.0105 0.0132 0.0169
N 0.0262 -0.0490 0.0028 0.0063 0.0068
U 0.1286 -0.0134 0.0367 0.0153 0.0396
内华达大地测量实验室点位结果/m E 0.1717 -0.1845 -0.0003 0.0183 0.0183
N 0.1188 -0.1233 0.0006 0.0170 0.0170
U 0.6814 -0.3894 0.0022 0.0512 0.0513
利用基于球谐展开的精密单点定位方法求得的汤加站点在2021年12月到2022年1月的站点点位数据结果如图 2所示，图 3为对照数据，为内华达大地测量实验室数据.两图中蓝色虚线表示1月15日火山最大规模喷发时刻，为UTC时的凌晨4点10分左右.解算原始数据在2022年1月26日以及27日有缺失；对照数据在2021年12月6日、2022年1月1日、12月26、27日有缺失.
图2 球谐展开单点定位坐标时序图

Fig 2 Spherical harmonic expansion single point positioning error timing diagram

图3 内华达大地测量实验室点位坐标时序图

Fig 3 Time sequence diagram of point error in Nevada geodetic laboratory

对照来看，相比于内华达物理实验室的误差数据，使用球谐展开的方法对TONG站点进行静态精密单点定位解算，在误差结果上精度达到厘米级别；在E方向以及N方向上，使用新方法得到的结果能够大致保持在2 cm左右，在天顶方向，此方法的精度保持在3.6 cm左右.整体来看，使用球谐展开式方法能够体现出向东南方向偏移的趋势，对比的内华达州大地测量实验室则并无趋势.并且从时序图中可以看出，在火山最大规模喷发时刻，解算结果在E方向以及U方向有明显跳动，E方向在火山做大规模喷发之前，误差保持在0左右，在最大规模喷发时，E方向有很明显的误差波动，在喷发过后不久，E方向上误差迅速有一个明显的突变，并且在这一个月内保持间断的误差并且存在很多跳变；U方向根据惠灵顿VAAC(Volcanic Ash Advisory Centers，火山灰咨询中心)的说法，在火山喷发后的一天后，火山灰的羽流就已经达到了30 km处左右，并且根据全球雷电数据集GLD360的监测，在火山喷发后，GLD360在羽流之中检测到了将近400000次的闪电事件，在一月喷发后的误差突增以及跳变可能是因为火山喷发的羽流以及频发的闪电影响了卫星信号在大气层中的传播；N方向上变化幅度较小，一直保持在0左右.
取喷发后三天的点位坐标时序数据进行拟合，来判断火山喷发对IGS站点位置的影响.对各个方向的时序数据进行简单的多项式线性一次拟合，拟合后的一次函数的斜率即为移动的速率，并且斜率的正负代表了移动的方向.通过拟合后的函数得出各个方向上偏移方向以及偏移速度，如图 4所示：E方向向东偏移，速度为0.0108 m/d；N方向向南偏移，速度为0.0023 m/d；U方向向下移动，速度为0.0149 m/d.
图4 站点移动趋势图

Fig 4 Site movement trend

3.2 时序数据SSA分析

使用奇异谱分析(Singular Spectral Analysis, SSA)来对解算得到的时序数据结果进行分析.奇异谱分析是一种处理非线性时间序列数据的一种方法，通过对所要研究的时间序列的轨迹矩阵进行分解、重构等操作，提取出时间序列中的不同成分序列，从而进行对时间序列进行分析和去噪(Schoellhamer，2001).利用W相关法判断主成分的阶数，然后原始数据减去主成分得到数据残差.
根据经验，原始数据的长度为M，则窗口长度L一般需要满足$\frac{M}{3} < L < \frac{M}{2}$, 经过多次尝试，剔除掉不合理数据后，将2.5天的原始数据长度作为周期，选择窗口长度为6948.从图 5中可以看出，N、E、U三个方向的数据主成分分解各不相同，E方向有2个频率分解较好的信号；N方向分解较好，在前20阶里包含了2个频率分辨率较好的信号；U方向有两个分辨率一般的信号.
图5 W相关分析结果

Fig 5 Results of W-correlation analysis

按照W相关分析图的结果，E、N方向取前三阶为主项，U方向取前四阶为主项；将主项重构，重构后的坐标时序图如图 6所示.E方向在1月之前周期平稳，振幅在0.045 cm，在15号左右开始变小；N方向上周期不稳定，但是在1月之后波动明显增大，振幅不稳定并且没有明显周期；U方向周期规律明显，在1月之前定位结果的振幅在0.05 cm左右，坐标结果平稳无异常；在1月8号左右周期振动幅度明显增大，在1月15号左右周期振幅发生明显变化，振幅达到0.15 cm，随后振幅变小.
图6 主成分重构图

Fig 6 Principal component reconstruction diagram

将原本的坐标时序数据减去主成分，得到坐标的残差数据，如图 7所示.在E方向上的残差没有明显响应；N方向的残差在12月20日左右有一个明显跳变，并且在1月15日有一个小的跳变；U方向1月15日有一个跳变，并且在之后有明显的波动.
图7 残差图

Fig 7 Residual plot

3.3 球谐展开误差项同火山灰分析

根据球谐展开式，来改正信号传播中有关于方位角以及高度角的误差项.将球谐展开变成关于卫星与接收机之间方位角以及高度角度的函数，利用求解得到的球谐系数来绘制改正量与方位角和高度角的变化图，此误差项的函数表达式为：
$\begin{align*}A(\varphi, \lambda)= & \sum\limits_{n=0}^{N_{\max }} \sum\limits_{m=0}^{n} P_{n m}(\cos (\theta))\left[C_{n m} \cos (m(\lambda))+\right. \\& \left.S_{n m} \sin (m(\lambda))\right], \end{align*}$
式中函数的自变量为φ(高度角)和λ(方位角), 其中θ=90°-φ.因为火山喷发后火山灰会进入大气成为信号传播路径上的误差，因此利用此函数可以探究火山喷发后火山灰对信号传播的影响.
UTC时间的13号15：20左右，汤加火山发生喷发，产生的蒸汽和气体羽流上升到海拔20 km.14日的傍晚开始，火山断断续续喷发，根据惠灵顿VAAC的报道，羽流达到了14 km.在UTC时间的15日凌晨4：00，发生了最大规模的喷发，火山喷发羽流上升的非常快，在5：19左右，羽流已经上升到15.2 km左右；到下午的14：43时，羽流已经到达30 km左右的高度.
根据VAAC的报道，选择火山最大爆发的1月15日的爆发节点进行绘图，得到误差随时间的变化图.UTC凌晨4点至凌晨5点45分的变化图如图 8所示.
图8 UTC时4点—5点45分误差图

Fig 8 Error chart from 4:15 to 5:45 at UTC

根据球谐展开进行绘图，将其表示为极坐标形式，其中角坐标表示为火山相对于TONG站点的方位角，径坐标表示为火山相对于TONG站点的高度角，球谐展开代表的误差大小利用颜色的不同表现出来，其中颜色越深表示误差越大，红色虚线表示了火山相对于IGS站点的方位角度.
根据报道，4时20分左右开始火山爆发，位置位于汤加站点的西北方向，在喷发初期，火山灰还不够高，在信号传播路径上的影响多是因为之前的小规模喷发所造成的，图中误差较大的方位逐渐向东偏移，符合VAAC的报道，喷发的火山灰向东漂移.
在UTC时5点左右，根据VAAC的报道，火山灰已经上升到了将近20 km，并且根据卫星图火山灰羽流的顶部直径已经达到了600 km，火山灰的扩散非常快.从误差图上可以看出在5点左右，火山相对于TONG站点的方位附近，高度角30°至40°之间，存在着非常大的误差，和报道相吻合.并且在低高度角处，误差相较于4点开始变大，说明火山灰相较于上一时间节点已经上升了，符合报道，说明火山灰的蔓延对卫星信号在传播过程中造成了不小的误差.
UTC时14点变化图如图 9所示.在喷发之后的14点40分左右，在火山相对于TONG站点的方位角度上，低高度角附近明显出现了很高的误差，说明在火山位置，火山灰已经上升到了非常高的高度，符合对于此时间段的火山灰报道.并且在更高的高度角处有着很明显的误差影响，说明火山灰蔓延的范围非常广.
图9 UTC时14点40分—14点50分误差图

Fig 9 Error chart from14:40 to 14:50 at UTC

4 讨论

因为GPS信号在卫星信号传播过程中，存在各种不同类型的误差延迟，并且存在类似对流层湿延迟的误差项难以利用模型模拟改正.利用卫星与地面接收机之间的方位角以及高度角来构建的球谐展开来对卫星信号传播的路径进行整体模拟，将球谐展开式线性化表达后作为误差项来构建精密单点定位的误差方程组来进行解算，可以避免这些问题，对难以模拟的误差延迟进行整体化处理.
通过利用球谐展开式构建误差方程对汤加站点进行解算，得到的站点坐标结果和内达华大地测量实验室的站点时序结果相比，定位精度大致相同，都可以达到厘米级别.球谐展开解算得到的站点坐标，E方向的平均误差为1.05 cm，N方向平均误差为0.28 cm，U方向的平均误差在3.67 cm；并且在标准差中各个方向的误差都达到厘米级别，相对于参考数据，此方法在定位上的精度可以得到保证.使用基于球谐展开的精密单点定位方法，利用站点点位观测数据以及卫星数据可以快速获取站点位置信息，对比使用气象数据构建经验模型的改正方法更加的方便快捷简单.并且球谐展开得到的坐标时序数据从一定程度上能够响应火山喷发，根据时序数据的波动以及SSA分解之后周期和残差的变化，能够判断火山的喷发以及对于站点位置的影响.
利用精密单点定位在传播路径上的误差改正分析，可以判断火山灰的扩散蔓延上升的情况.在对几个关键时间节点上的精密单点定位传播路径误差进行绘图，发现在火山喷发后，相关误差项的数值以及方位上的变化，较为符合现实中的火山灰上升速度以及方位.并且GPS卫星可以确保持续的时空覆盖率，可以做到长时间测量.

5 总结

本文使用基于球谐展开的精密单点定位方法对汤加站点进行解算，研究此方法对于火山喷发的响应以及监测的可行性.
使用此方法得到的E、N、U各方向的站点平均误差为1.05 cm、0.28 cm、3.67 cm；E、N、U各个方向上的标准差分别为1.32 cm、0.63 cm、1.53 cm.在精度上达到了精密单点定位的要求.
在点位误差的时序数据中，E、U方向上对于火山喷发的响应表现明显，在汤加火山的最大喷发时刻都表现出了较大波动，通过球谐展开单点定位坐标时序图可以看出，在15日爆发之后三天内，E、U方向上的坐标值相比于爆发时刻变化了大约4 cm.经过一次拟合得到站点向东南移动，向东速度为0.0108 m/d，向南移动速度为0.0023 m/d，向下移动为0.0149 m/d.经过SSA分解之后，各个方向的周期以及残差也都在火山喷发时出现了明显的变化.
对于精密单点定位中使用的球谐展开式所改正的大气传播误差，同火山灰的上升扩散情况也表现出来非常大的相关性.在喷发后对于火山方位出现了1.05 m的延迟误差，并且在火山喷发扩散一段时间后，低高度角同样出现了高误差，反算可得火山灰喷发高度已经达到30 km.
本文使用基于球谐函数的精密单点定位的方法对汤加火山进行喷发监测，此方法得到的时序点位数据中同火山喷发存在响应，同时利用球谐展开修正的传播路径误差也与火山喷发后火山灰的上升扩散趋势相关联，验证了使用GPS定位技术来研究火山喷发的可行性.同时此方法也存在一定局限性，在现有的改正方法下还难以对响应得到的火山喷发现象进行更为精细的探索以及定量的研究，后续工作还需要更加完善的处理方法来对于响应的误差项进行研究.

感谢武汉IGS数据中心、IGS提供的观测数据、精密星历、精密钟差等文件，感谢内华达州大地测量实验室的时序数据，感谢RTKLIB提供的开源代码，感谢审稿专家提供的修改意见和编辑部的大力支持.

References

Publishing order | Descend order by publishing year | Descend order by cited within
Aranzulla M , Cannavò F , Scollo S . Detection of Volcanic Plumes by GPS: the 23 November 2013 Episode on Mt. Etna. Annals of Geophysics, 2015, 57: 115- 132.

DOI

Blewitt G , Hammond W , Kreemer C . Harnessing the GPS data explosion for interdisciplinary science. Eos, 2018, 99

DOI

Carlin L , Hauschild A , Montenbruck O . Precise point positioning with GPS and Galileo broadcast ephemerides. GPS Solutions, 2021, 25 (2): 77

DOI

Dai L. 2002. Augmentation of GPS with GLONASS and pseudolite signals for carrier phase-based kinematic positioning[Ph. D. thesis]. Sydney: University of New South Wales, 126-129.

Deng Z , Bender M , Dick G , et al. Retrieving tropospheric delays from GPS networks densified with single frequency receivers. Geophysical Research Letters, 2009, 36 (19): L19802

DOI

Dong X H , Guo J Y , Kong Q L , et al. High-precision positioning method of miniature GNSS array based on side length constraint network adjustment. Journal of Shandong University of Science and Technology (Natural Science), 2021, 40 (2): 21- 30.

DOI

Feltens J . The international GPS service (IGS) ionosphere working group. Advances in Space Research, 2003, 31 (3): 635- 644.

DOI

Flynn L P , Harris A J L , Wright R . Improved identification of volcanic features using Landsat 7 ETM+. Remote Sensing of Environment, 2001, 78 (1-2): 180- 193.

DOI

Fournier N , Jolly A D . Detecting complex eruption sequence and directionality from high-rate geodetic observations: The August 6, 2012 Te Maari eruption, Tongariro, New Zealand. Journal of Volcanology and Geothermal Research, 2014, 286: 387- 396.

DOI

Grapenthin R , Freymueller J T , Kaufman A M . Geodetic observations during the 2009 eruption of Redoubt Volcano, Alaska. Journal of Volcanology and Geothermal Research, 2013, 259: 115- 132.

DOI

Gulbis A A S , Elliot J L , Person M J , et al. Charon's radius and atmospheric constraints from observations of a stellar occultation. Nature, 2006, 439 (7072): 48- 51.

DOI

Hopfield H S . Two-quartic tropospheric refractivity profile for correcting satellite data. Journal of Geophysical Research, 1969, 74 (18): 4487- 4499.

DOI

Houlié N , Briole P , Nercessian A , et al. Sounding the plume of the 18 August 2000 eruption of Miyakejima volcano (Japan) using GPS. Geophysical Research Letters, 2005, 32 (5): L05302

DOI

Jiang G W , Cheng C L , Wang B , et al. Effection of data quality relativity analysis of multipath. Progress in Geophysics, 2018, 33 (2): 461- 466.

DOI

Jin S G , Cardellach E , Xie F Q . GNSS Remote Sensing: Theory, Methods and Applications. Dordrecht: Springer, 2014 85- 99.

Klobuchar J . Ionospheric time-delay algorithm for single-frequency GPS users. IEEE Transactions on Aerospace and Electronic Systems, 1987, AES-23 (3): 325- 331.

DOI

Larson K M . A new way to detect volcanic plumes. Geophysical Research Letters, 2013, 40 (11): 2657- 2660.

DOI

Li Z H , Zhang X H . New Techniques and Precise Data Processing Methods of Satellite Navigation and Positioning. Wuhan: Wuhan University Press, 2009 85- 90.

Liu L , Guo J Y , Zhou M S , et al. A real-time correction method of GNSS broadcast ephemeris orbit. Journal of Shandong University of Science and Technology (Natural Science), 2021, 40 (3): 9- 17.

DOI

Liu Z M , Li Y Y , Guo J Y , et al. Influence of higher-order ionospheric delay correction on GPS precise orbit determination and precise positioning. Geodesy and Geodynamics, 2016, 7 (5): 369- 376.

DOI

Ning Y F , Yuan Y B . A modified geometry- and ionospheric-free combination for static three-carrier ambiguity resolution. GPS Solutions, 2017, 21 (4): 1633- 1645.

DOI

Ohta Y , Iguchi M . Advective diffusion of volcanic plume captured by dense GNSS network around Sakurajima volcano: A case study of the vulcanian eruption on July 24, 2012. Earth, Planets and Space, 2015, 67 (1): 157

DOI

Saastamoinen J . Contributions to the theory of atmospheric refraction. Bulletin Géodésique (1946-1975), 1972, 105 (1): 279- 298.

DOI

Schoellhamer D H . Singular spectrum analysis for time series with missing data. Geophysical Research Letters, 2001, 28 (16): 3187- 3190.

DOI

Shandong University of Science and Technology. 2022-03-11. GNSS single-point positioning method based on spherical harmonic expansion (in Chinese). China, CN113093242B.

Ware R , Exner M , Feng D , et al. GPS Sounding of the Atmosphere from Low Earth Orbit: Preliminary Results. Bulletin of the American Meteorological Society, 1996, 77 (1): 19- 40.

DOI

Yan J Q , Chen M J , Wang W , et al. Analysis of different troposphere correction comparison on PPP positioning effect. Engineering of Surveying and Mapping, 2016, 25 (4): 28- 32.

DOI

Yin Z X , Sun M B . Impact of high-order ionospheric delay on station coordinates in different areas. Science of Surveying and Mapping, 2018, 43 (4): 38- 42. 38-42, 65

DOI

Yunck T P , Liu C H , Ware R . A history of GPS sounding. Terrestrial, Atmospheric and Oceanic Sciences, 2000, 11 (1): 1- 20.

DOI

Zhang X H , Li X X , Li P . Review of GNSS PPP and its application. Acta Geodaetica et Cartographica Sinica, 2017, 46 (10): 1399- 1407.

Zhao C , Jiang C H , Liu J , et al. Comparative study of ionosonde TEC and GNSS-TEC. Progress in Geophysics, 2022, 37 (4): 1528- 1534.

DOI

Zumberge J F , Heflin M B , Jefferson D C , et al. Precise point positioning for the efficient and robust analysis of GPS data from large networks. Journal of Geophysical Research: Solid Earth, 1997, 102 (B3): 5005- 5017.

DOI

东 旭华 , 郭 金运 , 孔 巧丽 , 等. 基于边长约束网平差的微型GNSS阵列高精度定位方法. 山东科技大学学报(自然科学版), 2021, 40 (2): 21- 30.

DOI

蒋 光伟 , 程 传录 , 王 斌 , 等. 影响GNSS数据质量的多路径效应相关性分析. 地球物理学进展, 2018, 33 (2): 461- 466.

DOI

李 征航 , 张 小红 . 卫星导航定位新技术及高精度数据处理方法. 武汉: 武汉大学出版社, 2009 85- 90.

刘 路 , 郭 金运 , 周 茂盛 , 等. 一种GNSS广播星历轨道实时改正方法. 山东科技大学学报(自然科学版), 2021, 40 (3): 9- 17.

DOI

山东科技大学. 2022-03-11. 一种基于球谐展开的GNSS单点定位方法: 中国, CN113093242B.

闫 建巧 , 陈 明剑 , 汪 威 , 等. 不同对流层改正方法对PPP定位效果影响的对比分析. 测绘工程, 2016, 25 (4): 28- 32.

DOI

殷 志祥 , 孙 明博 . 高阶电离层延迟对不同区域测站坐标的影响. 测绘科学, 2018, 43 (4): 38- 42. 38-42, 65

DOI

张 小红 , 李 星星 , 李 盼 . GNSS精密单点定位技术及应用进展. 测绘学报, 2017, 46 (10): 1399- 1407.

赵 聪 , 姜 春华 , 刘 静 , 等. 基于测高仪的TEC与GNSS-TEC的比较研究. 地球物理学进展, 2022, 37 (4): 1528- 1534.

DOI

Options
Outlines

/

〈 〉