摘要:利用信噪比(SNR)数据进行逆建模可以反演海面高度及其变化,但反演精度、稳定性依赖初值的精度和SNR数据的时间连续性,基于多模多频SNR数据逆建模反演海平面变化的性能及其在潮汐分析中的应用有待进一步研究。本文为LSP(Lomb-Scargle periodog
本文内容来源于《测绘学报》2024年第11期(审图号GS京(2024)2421号)
基于多模多频SNR数据逆建模反演海面高度变化
陈灵秋1,2, 柴洪洲,1, 暴景阳3, 王敏1,2, 郑乃铨41.信息工程大学地理空间信息学院,河南 郑州 450001
2.
3.
4.
基金项目
作者简介
第一作者:陈灵秋(1988—),女,博士生,讲师,研究方向为GNSS遥感。 E-mail:
通讯作者: 柴洪洲 E-mail:chaihz1969@163.com
摘要
利用信噪比(SNR)数据进行逆建模可以反演海面高度及其变化,但反演精度、稳定性依赖初值的精度和SNR数据的时间连续性,基于多模多频SNR数据逆建模反演海平面变化的性能及其在潮汐分析中的应用有待进一步研究。本文为LSP(Lomb-Scargle periodogram)反演结果引入海面动态改正,用于逆建模过程参数初始化,获得稳定、均匀的高精度海面高度反演值,并由反演海面高度开展潮汐调和分析。选取MAYG、BRST和SC02这3个大潮差测站,对其1 a的多模多频SNR数据开展逆建模反演试验,与验潮站实测海面高度对比分析进行算法验证。结果表明,MAYG站逆建模反演海面高度的均方根误差(RMSE)为5.97 cm,BRST站为8.78 cm,SC02站为2.38 cm,海平面变化反演精度达到厘米级;与实测海面高度的潮汐调和分析结果对比,年度、逐月拟合残差中误差高度一致,提取的分潮振幅平均绝对误差(MAE)优于1 cm,提取的迟角MAE在3°以内,潮汐分析提取的潮汐成分和非潮汐水位均具有高度一致性。多模多频SNR数据逆建模反演海面高度能够替代验潮站实测海面高度用于潮汐调和分析。
关键词
GNSS-IR海面高度反演逆建模LSPSNR数据
本文引用格式
陈灵秋, 柴洪洲, 暴景阳, 王敏, 郑乃铨.
CHEN Lingqiu, CHAI Hongzhou, BAO Jingyang, WANG Min, ZHENG Naiquan.
海面高度是建立海洋垂直基准[1]、平均海平面模型和海洋潮汐模型[2]等的基础数据,对海洋科学研究[3]、全球气候变化研究[4]、海洋防灾减灾等具有重要意义,高时空分辨率的海面高度数据是建立高精度模型的前提[5]。沿海验潮站和卫星高度计是当今海面高度变化观测的主要手段。沿海验潮站可提供高精度和高采样的单点海面高度变化[6];卫星高度计可以按较理想的空间分布轨迹测量全球范围海面高度变化,在开阔海域能够实现2~3 cm的测高精度,但趋岸海域精度快速下降,在沿岸1~10 km海域测高精度为8~30 cm[7-8],且受卫星重访周期影响其时空分辨率受限。此外,配备GNSS的浮标也可用于海面高度测量,提供2 cm左右的测量精度[9],但是受通信、能源条件限制难以实现大规模布设和长期监测,主要应用于星载高度计的验证[10]。近年来,得益于GNSS技术的快速发展和不断完善,GNSS干涉遥感(GNSS interferometric reflectometry,GNSS-IR)技术被认为具有扩展沿海海面高度测量网络的潜力。
GNSS接收机接收卫星直射信号用于精确三维位置解算,反射信号与直射信号的干涉效应被量化于信噪比(signal-to-noise ratio,SNR)数据中。GNSS-IR技术对SNR数据进行处理,提取干涉信号的相位、振幅和频率等特征参数以反演测站周围地表环境变化,如海平面高度[11]、风速[12]、雪深[13]、土壤湿度[14]等。对干涉信号的频率参数进行估计可获得镜面反射点与天线相位中心之间的垂直距离,以实现海面高度及其变化测量。GNSS-IR反演海面高度变化的经典方法是利用LSP(Lomb-Scargle periodogram)频谱分析法进行频谱分析提取干涉信号频率,并通过对反演结果加海面动态改正[15]、对流层改正[16],以及利用多模多频反演结果在滑动窗口内进行最小二乘平差以提高精度和时间分辨率[17]。然而,在大潮差测站海面变化速率大,基于LSP频谱分析法的反演结果很难达到优于10 cm的精度,提高海面高度反演精度以达到验潮站的观测精度水平,仍是充分发挥GNSS-IR技术优势的重要研究目标。
GNSS-IR反演海面高度变化的另一类方法是对SNR数据进行逆建模,通过函数拟合来估计海面高度参数。SNR数据逆建模法基于多路径建模理论提出,并被应用于雪深估计[18-19]。考虑到海平面的动态变化,文献[20]在逆建模方法中引入B样条函数,采用GPS、GLONASS的L1、L2信号的SNR数据,在潮差1 m左右海域实现了优于5 cm的较高精度海面高度反演。文献[21]采用类似方法处理了5个站点的GPS L1信号的SNR数据,结果表明,与LSP频谱分析法相比,逆建模法精度能提高60%。文献[22]采用多个低成本天线的策略减小噪声影响,在潮差不足1 m的3个测站实现了高精度水位测量。在此基础上,文献[23]联合无迹卡尔曼滤波和逆建模方程,提出一种实时海面高度估计算法。
在GNSS-IR技术反演海面高度变化的理论发展过程中,国内外学者不断尝试将反演海面高度应用于潮汐调和分析。从GPS的SNR数据反演海面高度开展潮汐分析的结果来看,能够实现主要潮汐成分的提取[24],反演结果与验潮站实测数据的振幅和相位之间存在良好的一致性[16,25]。文献[26]表明,除太阴太阳赤纬全日分潮(K1)、太阳主要半日分潮(S2)和太阴太阳赤纬半日分潮(K2)外,其他主要分潮与验潮站结果之间存在毫米级差异。文献[27]在4个站点开展了GNSS-IR海面高度反演,也发现K1、K2潮汐成分存在误差,推断或与GPS轨道周期有关。文献[28]结合GPS、GLON ASS和Galileo的数据进行海面高度反演,与仅采用GPS相比,K1、K2分潮振幅的提取精度得以提高。
逆建模法为高度非线性函数的参数估计,依赖初值的精度和SNR数据的时间连续性以确保参数估计的稳定性。然而,大潮差海域的海面变化速率较大,难以保证初值的精度;多模多频SNR数据的运用能够改善数据的时间连续性,但是缺少该方面的研究。目前国内外对逆建模法反演海面高度的研究偏少,仅在小潮差测站采用了GPS L1或者GPS、GLON ASS的L1、L2信号进行研究。本文将逆建模法扩展为采用GPS、GLONASS、Galileo和Beidou的多频SNR数据,针对大潮差测站提出了可靠的初值设置方案,并对3个测站1 a的海面高度反演结果进行评估;此外,目前国内外针对LSP频谱分析法反演海面高度开展了潮汐调和分析相关研究,本文采用多模多频SNR数据的逆建模反演结果进行潮汐调和分析,评估其提取非潮汐水位和潮汐成分的性能。本文的研究工作展现了GNSS-IR技术拓展海面高度测量网络的潜力,有利于推动GNSS-IR反演海面高度的应用。
1 原理与方法1.1 SNR数据的逆建模原理GNSS接收机记录的SNR数据为信号功率与噪声功率的比值,对SNR数据进行二阶多项式拟合后便可去除直射信号部分的影响,余下部分可视作直射信号与反射信号的干涉结果,在低高度角范围相干性强,呈现信号震荡,如图1所示。该信号的频率与反射面高度有关,可建模为[15]
式中,DNSR为去趋势项的SNR序列;A为信号振幅,在进行高度反演时将其视作常数;H为天线相位中心到反射面的垂直距离,即相对于天线相位中心的海面高度;e为GNSS卫星高度角;λ为GNSS信号波长;ϕ0为信号相位残差。
图1图1 GNSS-IR反演海面高度变化原理
Fig.1 The inversion principle of sea surface height based on GNSS-IR
将DNSR序列的时间变量视作sin e,并考虑海面高度与时间有关,则信号的频率为
式中,和分别为海面高度和卫星高度角的时间导数。由于序列的时间变量不是均匀变化的,对于此不均匀时间序列采用LSP频谱分析法获得显著频率f。因此,顾及动态改正的海面高度为式中,可通过先验的海面高度变化、拟合潮波函数及最小二乘估计[29]
来确定,可得到顾及动态改正的LSP反演结果。SNR数据的逆建模法把式(1)中海面高度H的估计视作非线性函数的参数估计问题。以衰减因子表示DNSR振荡幅度随高度角增加而减小的过程,表示为[23]
(4)式中,κ是波数;s是反射面高度的标准偏差。由此,式(1)中DSNR建模为
式中,ts=sine;同相、异相项C1和C2代替式(1)中的振幅和相位,以保证建模过程中的数值稳定性;令Λ=s2也是基于相同的考虑。
为确保反射信号来自海面,DSNR序列在一定的高度角和方位角范围内截取,导致其在时间上是不连续的。因此,要获得连续的海面高度估计需要用一个函数来表示海面高度变化,然后估计该函数的参数来获得连续的海面高度变化。
任何考虑时变特征的分析函数都能表示海面高度随时间变化。B样条函数具有较好柔韧性,在考虑主要的半日潮和长期海面高度变化时能够提供足够的可变性[20]。p阶B样条函数是由n个控制点决定的p阶分段多项式函数[30],定义为
式中,H(u)是用于表示海平面随时间u变化的B样条曲线;PI是控制点,控制点数n=(tn+1-tp)/space,space为节点间距;Ni,p为节点向量T=[t1t2…tn+p]定义的p阶B样条基函数。为确保起止时刻的稳定性,采用准均匀节点向量,其前p个节点和后p个节点相同,即t1=t2=…=tp,且tn+1=tn+2=…tn+p,节点数为n+p。B样条基函数Ni,p递归计算得到
由一定的时间窗口[tp,tn+1]内的多模多频DSNR时间序列估计式(5)和式(6)中未知参数Λ、C1、C2和Pi。其中,未知参数Λ、C1、C2与卫星信号强度、接收机特性和反射面电磁属性等因素有关,较短的时间窗口内可视作常数。另外,C1、C2对于不同卫星系统频率其数值不同,如使用4个GNSS系统的相应频率数为fG、fR、fE和fC,则C1、C2的个数均为Nf=fG+fR+fE+fC,未知参数的总数为2Nf+n+1。设结合式(5)和式(6)的非线性函数Y(x)=[y(1 x)y(2 x)…ym(x)]T,m为时间窗口[tp,tn+1]中DSNR数据数量,未知参数x为
(8)
而
式中,j={1,2,…,Nf},为DSNRi的系统频率序号;、Hi是与时间有关的量,为DSNRi对应的时间和海面高度。对于特定的时间序列,反射面高度Hi(P1,P2,…,Pn)为控制点的函数,通过估计控制点参数获得海面高度变化。根据非线性最小二乘估计理论迭代计算未知参数,即(10)
高次的非线性方程组的迭代算法中,高斯-牛顿法是经典的方法之一,即最小二乘平差中所采用的方法。然而,当雅可比矩阵奇异或充分接近奇异时,高斯-牛顿方法难以求解,LM(Levenberg-Marquardt)方法是高斯-牛顿方法的一种改进[31-32]。LM算法的迭代公式为
(11)
式中,Yk=Y(xk);Jk=Y'(xk)是Y(x)在xk的雅可比矩阵;εk≥0是LM方法的迭代参数。当Jk奇异或者接近奇异时,引入的εk克服了迭代无法求解的困难。
逆建模法需设置未知参数的初值和建模的时间窗口大小。为了保证计算效率,将C1、C2初始值设置为0,Λ的初始值设置为1。控制点Pi决定曲线的走向,其初始值的准确性对拟合结果影响较大,本文采用顾及动态改正反演结果拟合B样条曲线得到Pi的初始值。时间窗口[tp,tn+1]和节点间距需根据反射面的变化特性确定,详细讨论见文献[22]。考虑潮汐半日周期的主要特征,本文窗口长度均设置为6 h,节点间距设置为2 h。如图2所示,为确保初值的稳定性,选择该时间窗口中间时段的LSP反演结果来进行B样条插值。经迭代拟合后输出该时间窗口内的参数估值,时间窗口后移其
步长,逐次计算每一窗口内未知参数。全部计算完成后,选取每个时间窗口内中间时段的控制点依据式(6)确定海面高度。
图2图2 逆建模时间窗口
Fig.2 The time windows of inverse modeling process
1.2 逆建模反演海面高度的潮汐调和分析为进一步评估逆建模法反演海面高度变化的能力,本文对反演结果开展潮汐调和分析。海水受到的引潮力是日地距离、日月距离、太阴赤纬、太阳赤纬等多个变量的非线性函数,在潮汐分析中海面高度被认为是一系列余弦函数线性叠加的结果[33]
式中,S0为平均海面高度;J为分潮数量;σj为第j个分潮的角速率,由地球、月亮和太阳的天体运行规律决定;Vj为第j个分潮的格林威治初始相角;Gj为第j个分潮的迟角;hj为第j个分潮的振幅;uj为第j个分潮的交点订正角;ζj为第j个分潮的交点因子。根据海面高度变化时间序列,通过最小二乘法估计各分潮振幅和迟角来实现海面高度变化预测。
2 测站与数据GNSS-IR技术反演海面高度的精度与测站潮差大小存在一定关系,对于经典LSP频谱分析法,潮差大的测站容易受到动态改正精度的影响导致反演精度受限。此外,逆建模法要对时间窗口内的DSNR进行非线性最小二乘拟合,时间窗口内的DSNR数据分布由该时间段内测站接收到的卫星信号决定,采用多系统多频的SNR数据能改善时间窗口内的DSNR分布。因此,为评估逆建模法海平面变化的反演性能,本文选取了3个潮差较大的沿海多GNSS测站(multi-GNSS experiment,MGEX)-MAYG站、BRST站和SC02站。
表1总结了3个测站的GNSS设备、潮差、距离海面的高度等测站信息。MAYG站位于法国马约特岛,该海域潮差约4.0 m,测站附近约10 m处有一验潮站Dzaoudzi。BRST测站位于法国的布雷斯特港,海域潮差约7.0 m,距离测站293 m处有一验潮站BREST。Dzaoudzi和BREST验潮站海面高度数据可在法国潮汐观测网络(réseaux de référence des observations marégraphiques,REFMAR)官方网站(
表1测站概况
Tab.1
测站名坐标潮差/m距海面高度/m方位角取值高度角取值设备数据间隔/s接收机天线MAYG12.78°S,45.26°E4.060°~180°,330°~360°5°~15°TRIMBLE ALLOYTRM59800.0030BRST48.38°N,4.50°W7.017130°~330°8°~20°TRIMBLE ALLOYTRM57971.0030SC0248.55°N,123.01°W3.5550°~240°5°~13°TRIMBLE NETR9TRM59800.8015新窗口打开| 下载CSV
为确保反射信号来自海面,需确定测站SNR数据的高度角和方位角范围,使用文献[34]提供的软件工具在谷歌地球上绘制第一菲涅尔反射区,如图3所示。根据测站周围环境,考虑逆建模过程中计算效率和拟合精度,高度角范围应在能保证LSP频谱分析的前提下尽量缩小以减小SNR数据的噪声影响,而方位角范围应尽量扩大以保证反演的数量。基于此考虑,BRST站高度角范围取值为8°~20°,方位角范围为130°~330°;MAYG站高度角范围取值为5°~15°,方位角范围为0°~180°和330°~360°;SC02站高度角范围取值为5°~13°,方位角范围为50°~240°。其中,SC02站可观测方位角开阔,海面无明显干扰物,具有较好的观测条件。
图3图3 第一菲涅尔反射区
Fig.3 Screenshot of the first Fresnel zones for MAYG and BRST
本文使用2023年MAYG、BRST站和2022年SC02站的GNSS观测数据,测站接收的GNSS信号和可用的SNR数据情况见表2。SC02站观测文件格式为RINEX 2.11,SNR数据类型标识与另两个测站不同,且不记录BeiDou系统数据。此外,由于GPS的S2W数据的频谱分析存在双峰问题[35],本文不再使用该数据。BRST、MAYG和SC02站可用的SNR类型分别为17、16和9个。图4为MAYG测站1 d的多模多频DSNR数据分布情况,不同GNSS系统的DSNR数据之间加入了80 dB-Hz的偏差。采用多模多频SNR数据使得滑动时间窗口中数据量增加,数据在时间轴上分布更密集,有利于逆建模反演的精度和稳定性。
表2测站GNSS信号频点和SNR数据情况
Tab.2
L1:1 575.42 MHZS1CS1CS1GPSL2:1 227.60 MHZS2XS2XS2L5:1 176.45 MHZS5XS5XS5GLON-ASSG1:1602+K×9/16 MHZS1C,S1PS1C,S1PS1G2:1246+K×7/16 MHZS2C,S2PS2C,S2PS2E1:1 575.42 MHZS1XS1XS1E5a:1 176.45 MHZS5XS5XS5GalileoE5b:1 207.140 MHZS7XS7XS7E5(E5a+E5b):1 191.795 MHZS8XS8XS8E6:1 278.75 MHZS6X——B1:1 561.098 MHZS2IS2I—B3:1 268.52 MHZS6IS6I—BeiDouB2:1 207.140 MHZS7IS7I—B1C:1 575.42 MHZ(BDS-3)S1XS1X—B2a:1 176.45 MHZ(BDS-3)S5XS5X—注:K的取值为0~23。
新窗口打开| 下载CSV
图4图4 MAYG测站2023年DOY191的多模多频DSNR数据分布
Fig.4 The temporal distribution of the DSNRs of multi-GNSS and multi-frequency at MAYG on DOY 191 of 2023
3 逆建模法反演海面高度逆建模法反演海面高度变化的数据处理流程如图5所示。首先由经典LSP频谱分析法得到海面高度反演值,频谱分析质量控制指标采用峰值信噪比和峰值最小值。然后,由反演值拟合潮波函数计算海面动态改正,详见文献[36]。本文中潮波函数考虑8个主要分潮的影响,拟合过程中以3倍拟合残差中误差(3m)为阈值剔除反演异常值,迭代拟合过程直至无异常值后计算顾及动态改正反演值。最后,在逆建模时间窗口内初始化待估参数,再由Matlab中非线性最小二乘解算函数lsqnonlin采用LM算法迭代估计未知参数,对各窗口输出的控制点参数进行B样条插值即可获得均匀的海面高度变化。
图5图5 数据处理流程
Fig.5 Data processing process
反演海面高度变化以天线相位中心为基准,对反演海面高度H与验潮站实测海面高度HTG之和取均值作为天线相位中心与验潮站零点之间的距离,则以验潮站零点为基准的反演海面高度H0为
(13)
式中,mean表示取均值。统一基准后,计算反演结果与实测海面高度之间的均方根误差(root-mean-square error,RMSE)和相关系数进行精度评估。
为突出采用多模多频SNR数据的优势,逐步增加不同GNSS系统和频率的SNR数据开展逆建模。以MAYG测站为例,其2023年7月的反演精度评估见表3。总体上看,随着所采用系统频率SNR数据的增加,逆建模滑动时间窗口内DSNR弧段数相应增加,RMSE在逐步减小,相关系数在逐步增大。因此,本文采用各测站全部多模多频SNR数据反演海面高度。此外,对MAYG测站DOY1逆建模时间窗口内的迭代次数和计算时间进行统计,仅采用GPS系统L1、L2、L5频SNR数据时,平均迭代次数为13.75次,平均计算时间为10.33 s;采用四系统的L1、L2、L5频SNR数据时,平均迭代次数为8.42次,平均计算时间约为21.08 s。采用多系统数据后较少的迭代次数便可收敛,但数据量的增加使得每次迭代计算耗时更长。
表3MAYG测站多模多频SNR数据的逆建模精度评估
Tab.3
GPSL1GPS-S1C248.300.843 75L1,L2GPS-S1C;GPS-S2X107.900.995 26L1,L2,L5GPS-S1C;GPS-S2X;GPS-S5X146.870.996 74GLO GPS GLO GAL GPSL1GPS-S1C;GLO-S1C;GLO-S1P720.180.968 56L1GPS-S1C;GLO-S1C;GLO-S1P;GAL-S1X914.630.982 97GPS GLO GAL BDSL1GPS-S1C;GLO-S1C;GLO-S1P;GA L-S1X;BDS-S2I149.410.992 95GPS-S1C;GPS-S2X;GLO-S1C;GLO-S1P;GLO-S2C;L1,L2GLO-S2P;GA L-S1X;GAL-S7X;BDS-S1X;BDS-S2I;453.530.999 14BDS-S6I;BDS-S7I;L1,L2,L5全部633.590.999 19新窗口打开| 下载CSV
将经典LSP频谱分析法、顾及动态改正LSP频谱分析法、逆建模法反演海面高度与验潮站实测海面高度进行对比,如图6所示,为体现反演结果的细节信息仅显示7 d的结果对比。由于在计算海面动态改正过程中加入了质量控制环节,测站环境较为复杂的MAYG和BRST测站的反演异常值大大减少,使得顾及动态改正海面高度与验潮站实测海面高度吻合度得以提升。以该结果为初值再次进行DSNR数据逆建模后,反演结果与验潮站实测结果几近重叠,为时间分布均匀的海面高度变化序列。
图6图6 GNSS-IR反演海面高度与验潮站实测海面高度对比
Fig.6 Comparison between sea surface height inversions based on GNSS-IR and sea surface height measurements from tide gauges
统计3个测站1 a的反演结果与实测海面高度之间的RMSE,如图7所示。BRST测站环境复杂、干扰因素较多,且本文采用的方位角范围较大,使得经典LSP频谱分析法反演结果精度较差。经质量控制并加入动态改正后,3个测站顾及动态改正LSP频谱分析法反演结果精度趋于稳定,MAYG、BRST和SC02站的RMSE分别为16.89、23.30和10.94 cm。此外,MAYG站平均每6 h的时间窗口内反演值约有56个,BRST站约70个。SC02站虽然可用SNR数据类型减少,但是可观测方位角较大,平均每6 h反演值约65个。顾及动态改正LSP频谱分析法反演结果保证了逆建模的时间窗口内有足够数量和精度的反演值。3个测站逆建模法反演结果精度如下:MAYG站RMSE为5.97 cm,相较顾及动态改正反演结果精度提升65%,相关系数为0.997 8;BRST站RMSE为8.78 cm,精度提升62%,相关系数为0.998 5;SC02站RMSE则达到2.38 cm,精度提升78%,相关系数达0.999 6,该测站观测条件较好,反演精度最高。在此类潮差较大、环境复杂程度不一的3个测站,逆建模法相较于顾及动态改正LSP频谱分析法反演精度提升60%~80%,单个反演值精度为2~9 cm。
图7图7 MAYG、BRST站2023年和SC02站2022年海面高度反演结果RMSE对比
Fig.7 RMSE of sea surface height inversions for MAYG and BRST stations in 2023, and SC02 station in 2022
图8(a)、(c)、(e)为范德卡斯特尔(Van de Casteele)图,其横坐标为反演海面高度与实测海面高度之间的偏差,纵坐标为验潮站实测海面高度,图形形状能有效揭示海面高度数据的误差特性[37]。图中逆建模反演海面高度偏差(深蓝)相较于顾及动态改正反演海面高度偏差(浅蓝)分布更加集中,出现了个别与海面高度相关的较大偏差。由于逆建模反演过程在滑动窗口中进行,若时间窗口中出现控制点初值精度较低、SNR观测值误差大或数据量少的情况,该窗口的控制点拟合精度将较低,插值后窗内结果呈现为整体性的偏离。但是该结果不对其他窗口的拟合造成影响,各窗口间拟合独立进行,保证了逆建模反演过程的稳定性,整体上看仍然能达到较高精度。图8(b)、(d)、(f)显示了逆建模反演海面高度与实测海面高度的相关性,且3个测站的回归分析结果特征相似,回归系数均略大于1,反映逆建模反演海面高度的尺度略小于实测海面高度,应是由GNSS信号对流层延迟等未改正的误差引起。
图8图8 反演海面高度的误差分析
Fig.8 The error analysis of sea surface height inversions
4 潮汐调和分析在开展潮汐调和分析之前,对逆建模反演海面高度和验潮站实测海面高度序列进行预处理。3个测站GNSS数据均存在短时缺失,由拟合的潮高函数进行插值补全逆建模反演海面高度,形成5 min或6 min间隔(与验潮站数据间隔相同)的时间序列。本文使用T_TIDE程序包,根据带有交点订正的潮高表达式依最小二乘法求解调和常数。
3个站点潮汐分析回报与原始序列的残差如图9所示,该残差体现站点海平面变化的非潮汐因素影响,反映站点海面对气象影响的静力学和动力学响应程度。逆建模反演海面高度和验潮站实测海面高度的潮汐调和分析残差总体上具有较好的一致性,但是逆建模反演海面高度的残差中呈现明显的高频信号,特别是海平面变化剧烈的BRST测站。该测站海域潮差高达约7.0 m,海面变化受气象影响也较为剧烈,非潮汐因素引起的海平面变化波动更大。对3个站点的实测海面高度和逆建模反演海面高度进行逐月潮汐分析,统计逐月拟合残差中误差的变化,并与ERA5提供的逐月风速数据进行对比,如图10所示。逐月拟合残差中误差的变化趋势与风速变化趋势较为吻合,呈现了气象因素与非潮汐水位的相关性;反演海面高度与实测海面高度的逐月拟合残差中误差变化高度一致,说明非潮汐水位对2个序列的影响高度一致,同时也说明了反演海面高度变化提取非潮汐水位变化的可靠性。
图9图9 潮汐分析回报与原始序列的残差
Fig.9 The errors between the prediction of tidal harmonic analysis and the original sequence
图10图10 逐月拟合残差中误差和风速对比
Fig.10 The monthly RMSE of the prediction of tidal harmonic analysis and the original sequence
表4统计了1 a数据的潮汐调和分析结果,3个测站的海面高度序列均分析得到59个分潮的调和常数,潮汐分析回报均能解释原始海面高度序列方差的98%以上。MAYG和SC02测站两个序列拟合残差中误差的差异为0.07和0.18 cm,BRST测站的相应差异也在1 cm左右,该测站反演的海面高度序列受非潮汐水位影响略大。
表4潮汐调和分析结果
Tab.4
MAYG验潮站实测5998.98.87逆建模反演5998.98.80BRST验潮站实测5998.916.57逆建模反演5998.717.63SC02验潮站实测5998.011.55逆建模反演5998.011.37新窗口打开| 下载CSV
按各分潮振幅贡献率排名,MAYG测站取前8个分潮的振幅、迟角对比见表5。除M2、S2、K1外,其他分潮振幅之间的差异在1 cm以内,分潮振幅差值的平均绝对误差(mean absolute error,MAE)为0.89 cm;除了分潮P1外,迟角之间的差异均在2°左右,迟角差值的MAE为2.05°。BRST测站潮差更大,取前12个分潮的振幅、迟角对比见表6。该测站除M2和M4外,其他分潮振幅之间的差异在1 cm以内,振幅差值的MAE为0.79 cm;除MU2、L2、M4外,其余迟角之间的差异均在2°左右,迟角差值的MAE为2.76°。SC02测站取前7个振幅、迟角对比见表7,除了分潮K1外,其他分潮振幅之间的差异在1 cm以内,振幅差值的MAE为0.40 cm;迟角之间的差异均在1°以内,迟角差值的MAE为0.26°。
表5MAYG测站实测海面高度和反演海面高度的分潮振幅、迟角
Tab.5
分潮名称振幅/cm迟角/(°)M2104.140101.3212.81926.34027.712-1.372S253.15851.6511.50765.86967.282-1.413N218.72018.5480.1726.9419.068-2.127K114.13712.9031.2344.2176.181-1.964K214.13014.0360.09461.53063.575-2.045O18.9538.1770.7766.9605.9830.977NU24.1723.9000.2725.6126.449-0.837P14.1033.8440.2595.442359.7525.690新窗口打开| 下载CSV
表6BRST测站实测海面高度和反演海面高度的分潮振幅、迟角
Tab.6
分潮名称振幅/cm迟角/(°)M2204.800202.1982.602105.745108.121-2.376S275.20774.3730.834145.625148.074-2.449N241.40640.6740.73287.88689.718-1.832K221.64121.1800.461142.788145.508-2.720MU28.7007.8810.819101.937105.252-3.3152N28.0927.2200.87277.74479.905-2.161NU27.7357.1360.59983.95386.322-2.369SSA7.4346.8410.59383.26284.007-0.745L26.7556.4280.327112.853119.561-6.708K16.5636.826-0.26373.20474.202-0.998O16.4666.558-0.092327.962328.277-0.315M45.8424.5091.33399.788106.965-7.177新窗口打开| 下载CSV
表7SC02测站实测海面高度和反演海面高度的分潮振幅、迟角
Tab.7
分潮名称振幅/cm迟角/(°)K175.62274.2721.350279.796279.1760.620M255.79755.4080.38910.17810.260-0.082O142.67242.2750.398258.082258.272-0.190P124.12123.7890.332279.097279.126-0.029S213.36513.1330.23234.91134.1830.729N211.96812.019-0.051343.233343.324-0.091Q17.2657.1990.066251.248251.345-0.098新窗口打开| 下载CSV
综上所述,逆建模反演海面高度的潮汐分析结果与验潮站实测结果对比,提取振幅精度在1 cm以内,提取迟角精度为0°~3°。考虑到测站潮差较大且GNSS-IR技术与验潮站观测环境不同,已然具有很好的一致性。目前研究认为GNSS-IR反演海面高度潮汐分析得到的分潮振幅差异或与卫星轨道周期有关,由表5—表7可知,振幅差异并未体现与卫星轨道周期有关的规律,本文采用多系统多频SNR数据,各卫星系统的轨道周期不同或能减轻单一轨道周期带来的系统性影响[26]。
5 结论和展望为评估多模多频SNR数据的逆建模反演海面高度变化的性能,对3个大潮差测站-MAYG(潮差约4.0 m)、BRST站(潮差约7.0 m)和SC02站(潮差约3.5 m)1 a的多模多频SNR数据进行处理和分析。根据顾及动态改正反演结果确定逆建模时间窗口内控制点参数初始值,获得了较为稳定的高精度海面高反演结果。结果表明,逆建模反演海面高度与验潮站实测结果对比,MAYG站的RMSE为5.97 cm,BRST站为8.78 cm,SC02站为2.38 cm。相较于顾及动态改正反演结果精度提升了60%~80%,单个反演值获得了2~9 cm的精度。多模多频SNR数据不仅增加了拟合时间窗口中观测值数量,还改善了观测值的时间分布,保证了反演过程的稳定性,整体上取得较高精度。
对逆建模反演海面高度和验潮站实测海面高度进行潮汐调和分析,结果表明,逆建模反演海面高度潮汐分析回报能解释原始海面高度序列方差的98%以上,年度、逐月拟合残差中误差高度一致,反映其非潮汐水位的提取具有较高可靠性;对提取的潮汐分潮振幅、迟角进行对比分析,分潮振幅精度优于1 cm,迟角精度为0°~3°。总体上看,逆建模反演海面高度和验潮站实测海面高度的潮汐分析结果具有高度一致性。
一般情况下,海面高度中的非潮汐影响通常为10~20 cm,而基于逆建模反演海面高度的精度能达到厘米级,明显小于海面高度中的非潮汐因素影响,因此,其可替代验潮站实测海面高度开展潮汐调和分析。然而,分析过程中发现逆建模反演海面高度与验潮站实测海面高度之间存在些许尺度的差异和整体的偏差,仍需进一步开展路径误差[38]、地面沉降[39]等带来的影响分析。此外,参照精密定位数据用作估计潮汐荷载位移[40],可进一步研究不同星座轨道几何结构对逆建模提取潮汐成分的影响。
来源:测绘学报