GNSS-IR测量水位的精度评估和站点对比：以中国南海北部和日本南部站点为例

叶脉，李琳琳,2*，彭冬菊，王培涛，邱强

（1.中山大学地球科学与工程学院，广东省地球动力作用与地质灾害重点实验室，广东珠海 519082；2.南方海洋科学与工程广东省实验室（珠海），广东珠海 519082；3.香港理工大学土地测量及地理资讯学系，香港 999077；4.国家海洋环境预报中心，北京 100081；5.中国科学院边缘海与大洋地质重点实验室, 中国科学院南海海洋研究所, 中国科学院南海生态环境工程创新研究院，广东广州 511458；6.南方海洋科学与工程广东省实验室（广州），广东广州 511458；7.中巴地球科学联合研究中心，巴基斯坦伊斯兰堡45320）

0 引言

在全球气候变化引起海平面上升的背景下，预计2050 年全球生活在沿海泛滥平原的人口可达到3.3～3.5 亿[1-2]。这些地势低洼的沿海地区将可能面临因海平面上升而逐渐加剧的海洋灾害侵袭，如极端台风引起的风暴潮、海域地震和火山活动触发的海啸灾害等[3-6]。我国以粤港澳大湾区为代表的东南沿海地区人口密集、经济发达，却长期遭受以风暴潮为主的多种海洋灾害威胁，极端灾害事件往往造成重大经济损失和人员伤亡[7]，例如2017 年在天文潮上涨时期登陆珠江口的台风“天鸽”引起的风暴潮和强降雨使澳门半岛一半以上的地区被淹没，最大淹没深度为3.1 m[8]。气候变化和沿海局地地表沉降引起的相对海平面上升会进一步加剧这些海洋灾害造成的威胁。全球气候变暖将导致到2100年时全球热带风暴的平均强度增加2%～11%，风暴中心100 km 内的降水率将增加20%左右[9]。假设2018 号台风“山竹”在2100 年（预计海平面相对现在上升0.9 m）发生并在极端天文高潮时登陆，则原本未受灾的澳门大学将被完全淹没，整个内港区域的淹没深度将达到2.5 m以上[10]；我国沿海相对海平面上升速度约为2.6 mm/yr，比全球平均上升速度高0.7 mm/yr[11]，海平面上升0.5 m 将使澳门遭受海啸的风险增加1倍[6]。因此，极端海平面监测对沿海地区受灾风险评估以及海洋灾害预警具有重要意义。

准确的海平面监测数据是建立准确的大洋潮汐模型、确定大洋环流以及中尺度气候模型的关键，对全球气候研究至关重要[12]。目前远海（与海岸的距离超过10～15 km）海平面监测主要依赖卫星测高，即测量相对于参考椭球的海平面变化（绝对海平面高度）。但在近岸地区，由于陆地污染或缺乏地球物理校正，卫星测高的结果并不可靠[13]。沿海地区海平面测量的主要设备是沿岸安装的潮位站，但潮位站监测存在一定的局限性。例如，在风暴潮、海啸等极端事件来临时，潮位站容易受到极端风浪破坏或断电而停止记录，例如在2005年飓风“卡特里娜”[14]、2007 年飓风“谷努”[15]、2011 年日本东北地震海啸事件[16]、2017 年台风“天鸽”期间均出现多个潮位站缺失记录的情况；受安装位置构造活动、人类活动诱发的地面沉降等影响，潮位站只能测量出相对海平面高度，需要和全球卫星导航系统（Global Navigation Satellite System，GNSS）台站绑定才能计算出绝对海平面高度[17-18]。因此，亟需发掘新的观测手段来提高沿海海平面观测能力。

自从研究人员发现全球定位系统（Global Position System，GPS）的直接信号和反射信号的干涉图样与其特征频率和GPS天线到反射面的距离有关后[19-20]，在前人的不断努力之下，全球卫星导航系统干涉反射计（Global Navigation Satellite System Interferometric Reflectometry，GNSS-IR）作为一种新兴的卫星遥感技术[21]，在测量积雪深度变化[22-24]、植被含水量[25]、土壤湿度[26-27]、水面高度[28-33]等方面有了许多研究和应用实例，具有广泛的应用前景。安装位置靠近岸边的GNSS台站可以通过GNSS-IR技术，利用信噪比（Signal to Noise Ratio，SNR）数据中记录的卫星直接信号与水面反射信号产生的干涉信号，提取出干涉信号的特征频率并计算反射高度（水面—天线的距离），从而获得水面高度的时间变化过程（见图1）。近岸GNSS 台站相对于传统潮位站具有以下优势：①台站可架设在岸边高地，与水面保持了一定的距离，被极端风浪破坏的风险极小；②台站能计算出绝对海平面高度[32]，并将全球海平面高度建立在同一参考系内，能同时测量地表沉降和构造活动；③现有的许多架设位置合适的测地GNSS 台站可满足反演要求，数据源多。在监测海平面高度方面，GNSS-IR 已被逐渐证实能以较高的精度监测短期潮位和长期水位变化[29-32]、反演风暴潮期间水位[34-36]和捕捉部分海啸信号[33]，从而在一定程度上和传统潮位站形成时空互补。已有研究显示近岸GNSS 站点的架设位置、设备条件和卫星信号接收系统等（见表1中相关信息）均会对GNSSIR 反演的效果造成不同程度的影响。为了量化评估影响GNSS-IR 反演效果的重要因素，本研究选取南海北部（见图2）和日本沿岸多个近岸GNSS站点，详细对比分析了GNSS台站架设和设备条件对反演结果的影响。我们首先针对香港HKQT 站点，采用GPS 与全球轨道导航卫星系统（Global Orbiting Navigation Satellite System，GLONASS）联合反演2016—2020 年的长期海平面变化，探究HKQT 站点是否具有长期监测海面的能力；然后选取我国南海北部的4个站点（广东省惠州市NHZH 站点、汕头市NSTO 站点、深圳市NSZN 站点和遮浪市NZLG 站点，见图2），将反演结果与HKQT 站点进行对比，分析不同设备条件对GNSS-IR 反演结果的影响；最后，选取日本南部的两个站点，首次使用潮位站缺失地区的GNSS 台站捕捉到完整的风暴潮波形以及反演风暴潮涨水前后几天的水位变化，进一步验证了架设位置合适的GNSS站点在缺失潮位站记录的地区具有提供珍贵的风暴潮期间水位记录、补充现有测潮网络的能力以及完善海洋灾害预警系统的潜力。

表1 前人GNSS-IR反演水位使用的站点信息Tab.1 Station information used by previous GNSS-IR inversion of water level

图1 GNSS-IR原理几何关系图（修改自PENG等[34]）Fig.1 Geometric diagram of GNSS-IR principle（modified from PENG et al.）

1 GNSS-IR 估测海面高度的原理及方法

图1 为GNSS-IR 原理的几何示意图。GNSS-IR仅需一根右旋圆极化（Right-Handed CircularPolarization，RHCP）的天线，同时接收卫星直接信号与海面的反射信号，对信噪比数据中二者的干涉信号的特征频率提取分析从而进行高度反演[28]。

本文使用科罗拉多大学Larson 教授团队开发并提供的开源模型gnssrefl（网址：https://github.com/kristinemlarson/gnssrefl），该模型可根据站点架设情况，通过设置坐标位置、反射高度、仰角范围、方位角范围等参数计算反射区域。在这些参数中，反射高度越大，反射区域半径越大；仰角越小，反射区域距离台站的水平距离越远；方位角决定接收反射信号的角度范围。gnssrefl 可从RINEX 数据和卫星轨道数据中提取出信噪比数据[37]，用于计算海面高度、积雪深度、土壤湿度等。根据仰角和直接信号与反射信号的路径差，由图1 的几何关系可以得到反射面高度的计算公式[34]：

式中：φ为直接信号与反射信号的相位差；λ为波长；H为天线到海面的距离为反射高度；θ为仰角（即直接信号的入射角）。将sinθ作为变量[38]，可以得到：

式中:f为特征震荡的频率。将式（2）简化后可得f与H的关系：

对于f的提取，首先将信噪比从对数尺度的分贝-赫兹转换为线性尺度的伏特/伏特，再利用低阶多项式去除直接信号的趋势项。将信噪比残差δ建模为[29]：

式中：λ为GNSS 信号的波长；A为振幅；φ为相位偏移。然后对残差序列利用Lomb-Scargle 周期图（Lomb-Scargle Periodogram，LSP）分析法求得f，从而计算出反射高度H，进一步计算潮位值。GNSS 4个系统（GPS、GLONASS、Galileo、BDS）不同波段的信号均可按照上述GNSS-IR 经典海面反演理论对海面高度进行估算[36]。

2 海平面高度反演实例与讨论

2.1 香港HKQT 站点在台风“天鸽”和“山竹”风暴潮期间潮位监测性能

HKQT 站台（22.291°N，114.213°E）位于香港鰂鱼涌附近（见图2），接收机型号为TRIMBLE NETR9，天线相位中心距离水面平均距离为6.4 m，属于香港卫星定位参考站网（Satellite Positioning Reference Station Network, SatRef），由香港地政总署（Survey and Mapping Office, SMO）运营和提供数据，附近2 m 处有Quarry Bay 潮位站提供实测潮位数据参考。本文选取2017 年8 月23 日对香港造成10.2 亿美元经济损失和近百人受伤的台风“天鸽”事件以及2018年9月16日登陆香港造成巨大经济损失和200 多人受伤的台风“山竹”事件[36]，对这两次极端风暴潮期间的海平面高度进行反演，并将风暴潮期间的反演结果与Quarry Bay潮位站的数据进行对比。

对于HKQT 站，选择仰角区间为4°～9°，方位角区间为-60°～105°，选用波段为GPS 的L1 和L5以及GLONASS 的L1 和L2，原始数据为5 s 采样率的接收机自主交换格式（Receiver Independent Exchange Format，RINEX）。不同频率波段测量得到的海平面高度会有差异，但这种差异比起误差本身来说极小，因此将3 个频率的测量结果相结合不会对最终结果产生太大影响，却可以极大地提高风暴潮期间的时间分辨率[34]。在LSP 分析过程中，为得到有效的反射高度，采取以下质量控制标准：选择起始仰角为4°～6°、终止仰角为7°～9°的卫星弧段；反射高度为2～8 m；峰值噪声比（Peak to Noise Ratio，PNR）大于2；GPS系统L1波段LSP频谱峰值>8 volt/volt，GPS 系统L5和GLONASS 系统L1、L2波段的频谱峰值>6 volt/volt。当风速＜17.5 m/s 时，利用信噪比分析法反演海面高度受到风速的影响不明显，而当风速＞18 m/s后，反射信号的功率会受到海面粗糙度的影响[39]，振幅值减小，波形也会变得较扁平[40]。在2017年台风“天鸽”和2018年台风“山竹”袭击香港时，最高风速分别高达28.3 m/s 和38 m/s，水面扰动较大，需要适当将峰值噪声比标准降低至1.5。反演结果见图3，其中彩色点为GNSSIR 计算出的由反射高度转化的相对海平面高度，黑色圆点（由于过度密集而看似曲线）为Quarry Bay潮位站实测值。图3a 为2017 年台风“天鸽”期间及前后7 d 内水位监测情况，风暴潮涨水过程中潮位上涨超过1 m，整个过程持续约7 h。我们将潮位站数据拟合为曲线，求得GNSS-IR 反演结果时间点的潮位站测量值，从而计算误差，整体结果相对于潮位站实测值的均方根误差（Root Mean Square Error，RMSE）为20.43 cm。2018 年台风“山竹”于2018 年9 月16 日袭击香港，从图3b 中可以看到，风暴潮造成涨水超过2 m，涨水过程持续近12 h，虽然峰值水位所在的2 d 内有一个6.6 h 的数据空档，但整体结果的RMSE 仅为15.54 cm。由于2017 年的数据质量较差，因此2018年台风“山竹”期间的反演结果优于2017年台风“天鸽”，但是两次事件的风暴潮波形都被GNSS-IR 很好地捕捉到。在GNSS-IR 反演结果中（见表2），2017 年台风“天鸽”期间GPS 系统L1波段测量值为152 个，L5 波段测量值为69 个，GLONASS 系统L1 波段测量值为118 个，L2 波段测量值为102 个，平均每天有63 个潮位值，峰值水位所在的2 d 内日均有57.5 个潮位值。2018 年台风“山竹”期间7 d 的GPS 系统L1 波段数据为202 个，L5 波段数据为98 个，GLONASS 系统L1 波段数据为146 个，L2 波段数据为150 个，平均每天有85 个潮位值，峰值水位所在的2 d 内有一个6.6 h 的数据空档，日均为61.5 个潮位值。在两次风暴潮事件期间，GPS 数据反演的潮位值个数分别占比50.11%和50.33%，GLONASS 数据反演的潮位值个数占比49.89%和49.67%。加入GLONASS 系统后，时间分辨率比单GPS 系统提高了近一倍。多模多频的GNSS-IR 比GPS-IR 有更多的信号种类和卫星弧段数量，有利于观测异常潮位的涨潮、峰值和落潮的完整过程，可以极大地完善风暴潮过程的监测。

表2 两次风暴潮期间HKQT站不同波段反演潮位值个数Tab.2 The number of tidal level values retrieved from different wavebands of HKQT station during two storm surges

图3 风暴潮期间HKQT站点反演结果Fig.3 Inversion results of HKQT station during storm surge

2.2 香港HKQT站点长期潮位监测能力评估

对2016—2020 年长期水位进行反演时，选取的质量检验标准与上文相同，应用大气折射和高度率进行校正[29]。这5 a 的GNSS 数据均有不同数量的缺失，会影响高度率校正的拟合计算。2016 年有3个空档，分别在年积日（Day of Year）第135 天缺失13 h，第152 天缺失53 h，第240 天缺失26 h；2017 年有22 个空档，年积日第10 天缺失22 h，第259 天缺失338 h，其余18 个空档分布较分散，均为5～8 h；2018 年有2 个空档，年积日第136 天缺失49 h，第259 天缺失6.6 h；2019 年仅有2 个5～8 h 的小空档，第342 天后数据都缺失，但因其在尾端因此不影响拟合；2020 年仅有1 个空档，第261 天缺失25 h。为避免因潮位站数据的质量问题造成偏差，因此将每日少于800 个数据点的天数去除，不参与对比。将5 a 的GNSS-IR 反演结果与由潮位站实测值计算出的日平均水位进行对比（见图4），从结果可以看出，2018 年、2019 年和2020 年的误差相对较低，RMSE分别为4.3 cm、4.3 cm 和4.7 cm。2016 年和2017 年的数据质量较差，误差相对较高，RMSE分别为6.5 cm和7.8 cm。5 a 的误差均未超过10 cm，其中2017 年

图4 2016—2020年日平均水位反演结果与潮位站实测值对比Fig.4 Comparison between the inversion results of daily mean water level and measured values of tide station from 2016 to 2020

8 月与2018 年9 月峰值处水位与潮位站记录相差较大，原因是风暴潮涨水期间GNSS-IR 反演的时间分辨率略有降低，导致平均水位低于潮位站结果。2017 年平均每天有61.7 个潮位值，其余4 a 的日平均数据点均为70～71 个。结果表明风暴潮期间与非极端天气条件下的时间分辨率差异较小，仅有6.81%和14.29%。

2.3 南海北部其他站点与HKQT站点监测性能对比

我们选取南海北部NHZH 站点反演的2017 年台风“天鸽”期间的水位结果与HKQT 站点进行对比。NHZH 站点（22.694°N，114.564°E）位于广东省惠州市，接收机型号为TPS NET-G3，GNSS 数据以及附近潮位站的数据由自然资源部海啸预警中心管理，GNSS 数据采样率为30 s。站点反射区见图2，选取的方位角范围为120°～360°，仰角范围为5°～15°，反射高度范围为4～10 m。数据中反射信号强度普遍极弱，质量检验过程中选取的峰值信噪比标准为2，频谱峰值标准为0。图5 为NHZH 站点的反演结果与潮位站数据对比，从中可以看出，NHZH站点的反射区优于HKQT 站，但反射信号功率极弱，难以筛选出有意义的反射高度，且由于该地为港口，反射信号会受到来往、停泊的船只干扰，导致反演结果中异常值较多；此外，由于接收机型号较老，仅有GPS 系统L1 波段的数据，导致时间分辨率与HKQT 站相差较大，7 d 的日均反演值仅为39 个，峰值水位所在的2 d 内日均仅有38.5 个潮位值。尽管如此，仍然可以从反演结果中捕捉到潮位变化的大致趋势，尤其是8 月24—26 日的反演结果与潮位站实测值基本吻合。由于NSZN 站点（22.587°N，114.273°E）、NZLG 站点（22.656°N，115.569°E）距离水面较高（与平均海平面的垂直距离分别为14.4 m 和18.4 m），而数据均为30 s 采样率，不满足奈奎斯特频率的要求[37]，无法反演出有效的海面高度；NSTO 站点（22.220°N，116.775°E）的数据质量极差，同样无法反演出有效海面高度，因此这3个站点的反演结果被舍弃。

图5 NHZH站点反演2017年台风“天鸽”期间水位结果Fig.5 Comparison between the inversion results of NHZH station and measured values of tide station

2.4 日本南部站点捕捉风暴潮信号实例

首先在内华达大地测量实验室（Nevada Geodetic Laboratory，NGL，网址：http://geodesy.unr.edu/index.php）数据库中对日本GNSS 站点进行筛选，挑选架设位置和反射区域能够用于GNSS-IR 反演的台站，最终在日本南部选取G212 站点与J425站点。根据日本近5 a 的台风数据集绘制出台风路径图，选取在日本境内造成明显风暴潮涨水且路径经过上述两个台站的台风事件，最终选择对2018年台风“谭美”期间的水位进行反演。

G212 站点与J425 站点均由日本国土地理院（Geospatial Information Authority of Japan，GSI）管理和提供数据，原始RINEX数据均为30 s采样率，只含GPS系统L1和L2波信号，其中G212站点有一个并行架设的潮位站NAHA，而J425站点附近150 km内没有公开的潮位站。G212 站点和J425 站点的地理位置、反射区以及台风“谭美”的移动路径见图6。G212 站点（26.213°N，127.665°E）位于日本那霸市，选取155°～305°的方位角，5°～15°的仰角，反射高度范围为4～9 m，质量检测过程设置峰值信噪比标准为2，GPS 系统L1 和L2 波段的频谱峰值标准均设置为0。J425 站点（34.472°N，134.314°E）位于日本小豆岛，方位角范围为160°～270°，仰角范围为5°～12°，反射高度范围为4 ～10 m，设置峰值信噪比标准为1.5，GPS 系统L1 的频谱峰值标准为4，L2 为0.8。G212站点接收到的反射信号功率极弱，需将频谱峰值标准设置为0，导致难以筛除异常值，但从图7 可以看出该站点的整体结果与潮位站实测值契合极好，完整地记录下9 月29 日持续近1 d 的风暴潮涨水过程；1 d 后，J425 站点在9 月30 日同样捕捉到风暴潮导致的水位变化过程，波形、涨水过程持续时间与G212 的反演结果极为接近，且波形更加完整。因此，J425站点在缺少潮位站的地区提供了重要的风暴潮记录，在空间上补充了潮位站未能覆盖的范围。

图6 选取的日本南部GNSS站点概况Fig.6 Overview of selected GNSS stations in southern Japan

图7 2018年台风“谭美”期间日本南部GNSS站点反演结果Fig.7 Inversion results of GNSS stations in southern Japan during 2018 Typhoon"Trami"

2.5 不同站点反演性能与自身架设条件的关系

表3 列出了本研究选取的4 个不同站点的架设条件和风暴潮期间反演结果的性能对比。GNSS-IR反演受到站点位置、设备条件、附近遮挡物等因素的影响。站点距离海面的水平距离越远，相同仰角的反射区与站点自身的水平距离越远；站点与海岸的水平距离越近，其反射区在海面的覆盖越大。站点位置决定了反演时可选取的方位角、仰角范围，从而影响可接收的卫星弧段数量，最终影响时间分辨率。例如NHZH、G212、J425 站点的位置相对于HKQT站点更靠近海面，视野更开阔，因此可用的方位角范围或仰角范围大于HKQT 站点。由于不同卫星系统开发和投入使用的时间不同（先后顺序为GPS、GLONASS、Galileo、BDS），不同波段的增设时间也有先后差异（先后顺序为L1、L2、L5）[41]，接收机型号同样会对反演结果产生较大影响。型号较老的接收机只能接收到GPS 系统信号，且接收到的波段数量有限，极大地影响反演时间分辨率。例如G212 站点和J425 站点虽然可用的仰角范围大于HKQT 站点，但接收的波段仅有2 个，最终时间分辨率与HKQT 站点大致持平；而NHZH 站点仅有1 个波段数据，仰角范围和方位角范围均远大于HKQT站，时间分辨率却仅有HKQT的一半左右。

表3 不同站点架设条件及风暴潮期间反演结果对比Tab.3 Comparison of erection conditions of different stations and inversion results during storm surge

反射区内的遮挡物会对反演的精度产生较大影响。NHZH 站点位于港口，其反射区内有来往船只，尤其是在风暴潮期间，大量船只被迫停靠在站点附近，极大干扰了卫星的反射信号，导致该站反演的误差超出正常站点两倍多。部分接收机的扼流线圈对反射信号功率的扼制较大，在反演时若选取较低的质量检测（Quality Check）标准，则难以剔除异常值，会对反演精度造成一定影响。

3 结论

海平面监测亟需精度高、范围广的监测新技术，其对气候、水文以及海洋灾害预警具有重要意义[24]。本文针对新兴的GNSS-IR 技术，通过多个实例量化分析近岸GNSS 站点的架设位置、设备条件和卫星信号接收系统等对GNSS-IR 反演的效果的影响。首先对中国香港HKQT 站点两次极端风暴潮事件期间以及2016—2020 年的数据进行GPS 与GLONASS 联合反演，数据结果表明，架设位置合适的站点（例如HKQT 站点）具备长期监测海平面变化趋势的能力，在反演过程加入GLONASS 系统后比单GPS 系统的时间分辨率提高近1 倍，风暴潮涨水期间GNSS-IR 的时间分辨率比平时仅有小幅度下降，两者分别为6.81%和14.29%；然后对2017 年台风“天鸽”期间广东省惠州市NHZH 站点的数据进行反演，与HKQT 站点的对比结果表明接收机接收的卫星信号波段数量会极大地影响反演时间分辨率，天线接收反射信号功率的强弱会对异常值的筛选产生较大影响；最后选取日本南部两个相隔数千里的G212 站点与J425 站点，对2018 年台风“谭美”期间的潮位进行反演，结果充分表明GNSS-IR能在潮位站记录的空白地区提供宝贵的风暴潮水位数据，可以为风暴潮建模研究提供关键的性能验证[42-43]。

结合本文的研究结果，我们对近岸GNSS 站点的架设提出以下参考建议：①设备应选取接收卫星信号波段尽可能多且没有特化对反射信号扼制效果的接收机；②卫星星座分布会导致其反射区出现一个空缺，架设站点时应在视野开阔的近岸，使空缺部分在陆地上，反射区覆盖在水面上，同时水面无凸出的礁石、船只等干扰物；③站点架设高度应结合架设位置、设备存储数据的采样率等因素综合考虑。架设位置与水面的水平距离较近时，高度选择相对自由，反之则需要提升站点高度以保证反射区能覆盖到水面。在奈奎斯特频率的限制下，设备存储数据的采样率为30 s时，站点与水面的垂直距离不应超过10～12 m，采样率为15 s时不应超过25 m，采样率为5 s 时不应超过80 m。具体的采样率和站点与水面垂直距离推荐范围的计算可参考Larson教授团队开发的网页程序（网址：https://gnssreflections.org/rzones）。