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）。