馮智杰,張雙成,3,周昕,陳雄川
( 1. 長安大學地質工程與測繪學院, 西安 710054;2. 自然資源部陜西西安地裂縫與地面沉降野外科學觀測研究站, 西安 710054;3. 地理信息工程國家重點實驗室, 西安 710054 )
GNSS 基準站網(wǎng)及其相關數(shù)據(jù)處理系統(tǒng)能夠為研究人員提供迅速、可靠、有效的信息服務,廣泛地滿足基礎測繪、環(huán)境監(jiān)測、災害預警、氣象預報等信息需求.從20 世紀90 年代初期開始,包括我國在內的世界上多個國家、地區(qū)及相關國際組織相繼建立了高精度GPS 連續(xù)觀測網(wǎng)絡[1-2].隨著網(wǎng)絡技術與通信技術的不斷發(fā)展,GNSS 技術與其結合發(fā)展出具有網(wǎng)絡實時差分定位(real time kinematic,RTK)功能的衛(wèi)星導航定位連續(xù)運行參考站(continuously operating reference stations,CORS)系統(tǒng),其憑借快速高效、高精度、高可靠性、網(wǎng)絡化、自動化和智能化等優(yōu)點,給包括測繪行業(yè)在內的多個領域帶來巨大影響.
CORS 站作為一種地面GPS 接收機,可提供對地球表面位置和運動的精確和連續(xù)測量.Rwabudandi等[3]利用2015—2018 年印尼地理空間信息局提供的10 個CORS 站數(shù)據(jù)對東爪哇北部地區(qū)的地殼形變進行分析,結果表明CORS 站數(shù)據(jù)能夠較好地解譯出站點在地球表面位置的水平移動以及豎直方向上的沉降運動.Cie?ko 等[4]利用3 個相互間隔適當?shù)腃ORS 站作為流動站估計GPS 位置,并利用卡爾曼濾波處理數(shù)據(jù),驗證了基于CORS 站的測試區(qū)域可有效獲得GPS 位置.徐創(chuàng)富等[5]結合臺州15 個CORS站2017—2019 年的時間序列,利用頻譜分析及小波變換分析CORS 站周期特征,得到其在垂直方向相較于水平方向變化明顯很多,且更具有周期特征.李婉秋等[6]提出一種多通道奇異譜分析(multi-channel singular spectrum analysis,MASS)方法,該方法能有效提取CORS 站網(wǎng)大地高時序形變信息,其利用38 個CORS 站進行時序形變信息提取,92%的測站時序形變序列有不同程度的改善.廖超明等[7]開展沿海區(qū)域CORS 站穩(wěn)定性實驗,借助于GAMIT/GLOBK進行穩(wěn)定性分析,得到CORS 站對于解譯站點形變時間序列具有較好的穩(wěn)定性.以上實驗均在一定程度上驗證了CORS 站數(shù)據(jù)用于解譯站點形變序列的穩(wěn)定性和準確性.
然而,CORS 站分布在一個廣泛開闊的區(qū)域,擁有一個可用于跟蹤環(huán)境參數(shù)變化的參考點網(wǎng)絡,其不僅能夠測量地表位置和運動的變化,還可通過將CORS 數(shù)據(jù)與遙感數(shù)據(jù)相結合,綜合監(jiān)測地球表面環(huán)境參數(shù)變化情況.CORS 站所接收的衛(wèi)星反射信號中存在信噪比(signal to noise ratio,SNR)數(shù)據(jù),而SNR數(shù)據(jù)可解譯出土壤濕度及積雪深度等地表環(huán)境參數(shù).2008 年,Larson 等[8]首次提出了基于大地測量型GNSS 接收機的SNR 數(shù)據(jù)估計土壤濕度算法,其將直射信號與反射信號分離后發(fā)現(xiàn)反射信號振幅波動與土壤濕度變化存在良好一致性,驗證了利用衛(wèi)星SNR 數(shù)據(jù)解譯土壤濕度的可行性.Chew 等[9]分別研究了衛(wèi)星SNR 數(shù)據(jù)的相位偏移量及振幅對土壤濕度變化敏感性,實驗得到相位是裸露土壤條件下反映土壤濕度變化的最佳參數(shù).Ran 等[10]提出一種圓弧編輯方法,只保留為余弦波形的DSNR 數(shù)據(jù),驗證了其能夠提高在起伏地面的土壤濕度解譯精度.Li 等[11]利用Lomb-Scargle 頻譜分析等多種方法對PBO 網(wǎng)站中的CORS 站進行雪深反演實驗,其結果表明:利用SNR 數(shù)據(jù)解譯得到的積雪深度與實際測量值的相關性達到0.9 以上,能夠很好地反演積雪深度的趨勢信息.Hu 等[12]利用相同的CORS 站進行實驗,并提出一種具有自適應高通濾波器特性的變分模態(tài)分解算法,實驗結果與原位積雪深度高度吻合.An 等[13]提出了一種改進的利用GNSS 觀測數(shù)據(jù)估算積雪深度的方法,其在LS 譜分析前加入小波分解算法,能夠有效分離噪聲功率,提高了積雪深度變化劇烈時的反演精度.衛(wèi)星信號穿過對流層時,會受到大氣的影響出現(xiàn)路徑彎曲和速度衰減,從而使觀測值中含有一定的大氣延遲[14],這些延遲均會被CORS 站接收并保存于數(shù)據(jù)中,為大氣水汽探測技術提供了重要的數(shù)據(jù)支撐.施闖等[15]對中國-中南半島地區(qū)地的測站數(shù)據(jù)進行水汽反演,研究了該地區(qū)大氣水汽平均含量、年周期振幅和半年周期振幅等氣候特征.劉夢杰等[16]提出一種臺站處分布式北斗/GNSS 實時水汽解算模式,并在中國三個典型區(qū)域的測站進行實驗,其結果與實時水汽產(chǎn)品對比相關性達到0.95.上述實驗均已驗證利用CORS 站監(jiān)測地表環(huán)境參數(shù)提供了一種高效的數(shù)據(jù)收集方式,將CORS 數(shù)據(jù)與全球導航衛(wèi)星系統(tǒng)干涉反射測量(GNSS interferometric reflectometry,GNSS-IR)技術相結合,研究人員可以更全面地監(jiān)測地表環(huán)境參數(shù)變化.
為進一步推進CORS 站用于地表環(huán)境參數(shù)監(jiān)測,本文提出一種集“積雪深度、土壤濕度、大氣水汽及地表形變”的地表環(huán)境多參數(shù)綜合監(jiān)測體系,加強CORS 站在地表環(huán)境監(jiān)測方面功能.文章基于齊齊哈爾CORS 站BFQE 進行實驗驗證,采用GNSS-IR 技術及GNSS-PWV(precipitable water vapor,PWV)技術,對CORS 站的接收數(shù)據(jù)進行地表環(huán)境參數(shù)解譯監(jiān)測,獲得特定時間的土壤含水率、積雪深度以及全年的大氣水汽解譯結果,最后在討論章節(jié)利用GAMIT基線解算獲得其地表形變解譯結果,并在各解譯參數(shù)結果中結合降雨數(shù)據(jù)分析其與降雨量之間的響應關系,充分發(fā)揮CORS 站在地表環(huán)境綜合監(jiān)測方面的作用.
CORS 站所接收的GNSS 衛(wèi)星反射信號能夠反應反射面的物理特征的變化,而在實際環(huán)境中,空氣中的水汽含量[17],土壤中的水分含量[18],地表的積雪深度[19]以及細微形變[20],均會造成GNSS 衛(wèi)星信號的改變.本文利用GNSS 衛(wèi)星信號可監(jiān)測反射面物理參數(shù)變化的特點,提出了一種多參數(shù)的地表環(huán)境綜合監(jiān)測體系,包括土壤濕度、積雪深度和大氣水汽多個地表環(huán)境參數(shù),如圖1 所示.首先,利用云服務器遠程獲取CORS 站所接收的GNSS 觀測數(shù)據(jù),包括觀測文件、廣播星歷、精密星歷,其次對GNSS 觀測數(shù)據(jù)進行綜合處理分析,對重采樣的SNR 數(shù)據(jù)采用非線性最小二乘及Lomb-Scargle 譜分析等GNSS-IR 技術解譯特定時間段的淺層土壤濕度及地表積雪深度,然后通過聯(lián)測遠距離IGS 站采用GNSS-PWV 技術獲取測站大氣水汽時間序列.最后,利用降雨數(shù)據(jù)進一步對解譯的地表環(huán)境參數(shù)結果進行相關性分析,獲得各參數(shù)之間的響應關系,為CORS 站用于地表環(huán)境參數(shù)監(jiān)測方面提供案例支撐.
圖1 地表環(huán)境參數(shù)綜合監(jiān)測
GNSS-IR 技術利用衛(wèi)星信號在GNSS 接收機天線處疊加產(chǎn)生的干涉現(xiàn)象來提取地表反射面的物理特征.與傳統(tǒng)的衛(wèi)星遙感手段比較,其無需單獨制造特定的發(fā)射機,大大降低了成本,且擁有及其豐富的免費導航衛(wèi)星信號源,利用L 波段進行觀測,受大霧和雨雪等氣候影響較小,對地表參數(shù)的變化十分敏感.對于在典型的GPS 應用中,反射信號的多路徑效應被認為是誤差的來源,而在2009 年,Larson 發(fā)現(xiàn)SNR 數(shù)據(jù)與積雪深度具有相關性并利用實驗證明其可用于解譯積雪深度[21].SNR 是衡量GNSS 信號強度的量化指標,其主要受接收機天線增益、接收機噪聲和多路徑效應等因素影響,多路徑效應影響尤為顯著.因此,可通過分析SNR 中的干涉信號來獲取地表環(huán)境參數(shù),這極大的節(jié)約了監(jiān)測成本,實現(xiàn)了“一機多用”的成效[22].圖2 為利用GNSS-IR 解譯積雪深度及土壤濕度示意圖.
圖2GNSS-IR 解譯積雪深度及土壤濕度示意圖
圖2 中H為接收機天線相位中心到反射面的高度,θ為衛(wèi)星高度角,紅色虛線直反射信號間的路徑差為2Hsin θ,接收機天線所接收到的信號是由直射信號的反射疊加產(chǎn)生的干涉信號,假設平面僅發(fā)生一次鏡面反射的簡化模型,SNR 可表示為[23]
式中:Ad和Am分別為直射信號和反射信號的振幅;ψ為直反射信號之間的相位差.
利用低階多項式擬合SNR 得到衛(wèi)星信號的直接分量,減去直接分量后的SNR 即為反射分量,可表示為
式中: SNRm為反射分量;λ為衛(wèi)星信號載波波長;φ為反射分量相位.反射分量是以衛(wèi)星高度角的正弦值為自變量的擬余弦函數(shù).對其利用Lomb-Scargle 譜分析方法可獲得峰值頻率,結合頻率和波長關系式可得反射高
式中:h為反射高;f為利用Lomb-Scargle 譜分析方法獲得的峰值頻率.
將接收機高度天線高度H和反射高h帶入GNSS-IR 積雪深度測量的幾何模型中即可得到積雪深度結果?h
利用GNSS-IR 解譯土壤濕度原理為:淺層土壤濕度的變化會引起土壤中介電常數(shù)發(fā)生改變,從而導致GNSS 衛(wèi)星反射信號的SNR 數(shù)據(jù)產(chǎn)生波動,該波動可被接收機記錄儲存并加以利用;對低高度角的SNR 數(shù)據(jù)進行分解、分析處理即可解譯得到土壤濕度結果.
式(2)中SNRm為反射信號SNR 數(shù)據(jù),θ為衛(wèi)星高度角,兩者均為觀測數(shù)據(jù),對其利用非線性最小二乘即可求解得到相位觀測量φ,而相位是土壤濕度的最佳指示器,可用一個簡單的線性函數(shù)來解譯土壤濕度
式中: SMCt為衛(wèi)星解譯的土壤濕度值;S為相位φ與實測土壤濕度之間斜率;SMCresid為相位φ與實測土壤濕度之間截距.
在進行GNSS-IR 解譯土壤濕度時,需要確定GNSS 衛(wèi)星的反射區(qū)域面積,即GNSS 反射遙感的空間分辨率.根據(jù)慧根斯—菲涅爾原理,鏡面反射產(chǎn)生的反射區(qū)域被稱為菲涅爾反射區(qū),其大小可以看作GNSS 反射遙感的有效探測區(qū)域[24],菲涅爾反射區(qū)域是一系列橢圓,橢圓的長半軸a和短半軸b分別為[25]:
圖3 為GNSS-PWV 技術解譯大氣水汽示意圖,其原理是通過聯(lián)測遠距離IGS 測站利用精密單點定位(precise point positioning,PPP)解算GNSS 觀測數(shù)據(jù),可得到由對流層靜力學延遲(zenith hydrostatic delay,ZHD)和對流層濕延遲(zenith wet delay,ZWD)組成的高精度對流層天頂總延遲(zenith total delay,ZTD)[26].
圖3 GNSS-PWV 解譯大氣水汽示意圖
其中,ZHD 可根據(jù)Saastamonien 模型計算,公式為
式中,Pc、φc、Hc分別為測站氣壓、緯度、海拔.
利用ZTD 與ZHD 相減可得到ZWD
而PWV與ZWD存在線性關系,公式為
式中:K為大氣水汽轉換系數(shù);Rv為水汽的比氣體常數(shù);Tm為加權平均溫度;k3與均為大氣折射率實驗常數(shù).
齊齊哈爾市的BFQE 測站地處47.376269°N,123.923742°E,海拔約為148m.如圖4 所示,實驗區(qū)域位于一處農(nóng)田,所處位置地勢平坦,四周無高大建筑阻擋,GNSS 反射信號接收較為良好;四周僅有少數(shù)低矮農(nóng)作物,GNSS 信號受植被影響較少.
圖4 BFQE 測站環(huán)境圖
黑龍江省齊齊哈爾市春季干旱多風,夏季炎熱多雨,秋季短暫,冬季干冷漫長.本文收集了2019—2021 年BFQE 測站附近的氣象數(shù)據(jù),結果如圖5 所示.齊齊哈爾市冬季溫度低且時間較長,低于0°的天數(shù)約占全年天數(shù)的1/3.低溫為積雪提供了優(yōu)質的儲存環(huán)境,積雪深度在短時間內較為穩(wěn)定,給GNSS 雪深測量創(chuàng)造了良好條件.同時夏季的齊齊哈爾市降雨充足,相對濕度較高,受降雨影響,土壤水分含量變化較為明顯,適合于GNSS 土壤濕度變化的探測.
圖5 2019-2021 年氣象序列圖
根據(jù)歷史數(shù)據(jù)顯示,2021 年年積日第1—61 天和第310—365 天均存在積雪情況,這段時期可作為監(jiān)測積雪深度的時間段;而第89—300 天,測站地面則沒有觀測到積雪現(xiàn)象,這一時間段可用于監(jiān)測土壤濕度的變化.
圖6 展示了衛(wèi)星反演雪深與實際雪深的對比圖及相關性分析圖.在驗證2021 年存在積雪的日期時,使用10°~25°及10°~30°的低高度角數(shù)據(jù)進行反演.結果表明,這兩組數(shù)據(jù)的相關性均為顯著,分別達到0.91 和0.88,RMSE 分別為2.06cm 和2.11cm.無論是10°~25°還是10°~30°的結果,在雪深的變化趨勢上均與實際雪深相似,特別是在第312—313 天雪深急劇上升的階段,反演數(shù)據(jù)都能很好地反映出其上身趨勢.然而,也需要指出,在一些天數(shù)上,反演數(shù)據(jù)與實際雪深存在一定差異.例如在第344—355 天,實際雪深約為10cm,而兩組數(shù)據(jù)的反演雪深都偏低,約為8cm.
圖6 2021年雪深反演
冬季的降雨量會部分轉換為積雪,從而增加積雪深度.圖7 顯示了利用10°~25°衛(wèi)星數(shù)據(jù)反演得到的積雪深度與實測值的對比圖,藍色柱狀圖代表降雨量.在第1—61 天以及第310—365 天中,共出現(xiàn)5 次明顯的降雨事件,對應時間的積雪深度均呈現(xiàn)上漲趨勢.衛(wèi)星反演結果與實際積雪深度在降雨后幾乎同步上漲.在第312 天及313 天,分別出現(xiàn)了15.7mm 和7.1mm 的降雨,衛(wèi)星反演的積雪深度與實際積雪深度同時達到峰值,分別為20cm 和18cm.
圖7 2021年積雪與降雨關系圖
綜上所述,利用GNSS-IR 技術解譯得到的積雪深度結果雖然在數(shù)值上與實測積雪深度存在細微差別,但整體趨勢上與實際情況相符.同時,該技術能夠很好地反應降雨后的積雪深度變化,并與降雨存在較好的響應關系.
利用GNSS-IR 解譯土壤濕度前,必須首先確定實驗區(qū)域的第一菲涅爾反射區(qū),即此次實驗中GNSS-IR 解譯土壤濕度的有效區(qū)域.圖8 展示了利用CORS 站BFQE 數(shù)據(jù)繪制的不同低高度角(5°~30°)下各衛(wèi)星的第一菲涅耳反射區(qū).該圖結合谷歌地圖繪制而成,圖中標簽位置為CORS 站BFQE 所在位置,不同顏色的橢圓代表了不同的衛(wèi)星.
圖8 BFQE測站菲涅耳反射區(qū)
圖9 展示了利用4 顆GPS 衛(wèi)星進行土壤濕度反演的結果及相關性分析.衛(wèi)星反演得到的土壤濕度結果與實際測量的土壤含水率在整體趨勢上呈現(xiàn)較好的一致性,隨著實測土壤含水率的增加和減少,相關系數(shù)均達到顯著水平,分別為0.70、0.74、0.71、0.81.然而,在特定時期,衛(wèi)星反演的土壤濕度結果與實測值存在一定差異,甚至出現(xiàn)相反的趨勢.例如,在第119—130 天期間,展示的衛(wèi)星反演結果與實際測量值對比均存在一定的誤差;在第206 天,實際測量值下降至極小值點,而5、6 號衛(wèi)星反演結果的趨勢與其相反,呈現(xiàn)極大值點,但10、31 號衛(wèi)星反演結果與實際測量值相似.因此,基于單顆衛(wèi)星的反演結果時常會受到環(huán)境因素或其他因素的影響而出現(xiàn)一定的誤差,不能完全詳細反演土壤濕度的細微變化.
圖9 2021年土壤濕度反演
利用GNSS-IR 技術解譯得到的土壤濕度反映了CORS 站周圍淺層土壤的濕度情況,這些數(shù)據(jù)距離地表僅有5~15cm 的深度,因此受降雨過程中雨水入滲的影響較大,土壤濕度含量與降雨之間存在著密切的關聯(lián).為探討之間的存在的響應關系,依據(jù)相關系數(shù)將多顆衛(wèi)星反演土壤濕度數(shù)據(jù)進行加權融合處理得到的結果與降雨數(shù)據(jù)結合得到圖10.可以發(fā)現(xiàn)土壤濕度數(shù)據(jù)隨著降雨的出現(xiàn)而產(chǎn)生增長的趨勢,其中,第182 天出現(xiàn)了78.2mm 的大幅度降雨,土壤濕度從前一天的11.87%增長至16.46%,而衛(wèi)星反演數(shù)據(jù)雖存在一定的延遲現(xiàn)象,但同樣出現(xiàn)增長情況,從13.51%增長至頂峰的17.82%.到第220 天,降雨量達到52.2mm,實測土壤濕度數(shù)據(jù)到達峰值23.22%,而衛(wèi)星反演數(shù)據(jù)的峰值19.26%延遲到第226 天出現(xiàn).多星融合反演結果相比于單星反演結果在土壤濕度的趨勢上更加擬合實測值,相關性也更高,達到0.89.與單星反演的平均相關性相比,多星融合反演的相關性提高了18.7%.
圖10 2021年土壤濕度與降雨關系圖
綜合來看,GNSS-IR 技術具備解譯CORS 站周圍淺層土壤濕度的能力.該技術解譯得到的土壤濕度結果與實測土壤濕度的趨勢相一致,并且在反應降雨對土壤濕度的影響關系方面表現(xiàn)較為出色.
圖11 展示了結合降雨數(shù)據(jù)繪制的大氣水汽含量圖,其中2021 年BFQE 測站的大氣水汽呈現(xiàn)倒V 字趨勢.在夏季時段,大氣水汽含量豐富,數(shù)值保持較高水平,在第196 天達到峰值55.63mm.而在春冬季節(jié),大氣水汽含量開始降低,保持在小于20mm 的低水平,第289 天降至最低值3.53mm.
圖11 2021年大氣水汽反演
降雨是指大氣中冷凝的水汽下降到地球表面的天氣現(xiàn)象,大氣中水汽的含量直接影響降雨量.由圖11可知,降雨出現(xiàn)前,大氣水汽含量呈上升趨勢,而降雨后則出現(xiàn)驟減現(xiàn)象.例如,在第165、181、220 及235 天,大氣水汽持續(xù)積累,導致降雨現(xiàn)象的出現(xiàn).其中,第220 天大氣水汽從前一天的27.44mm 增長至43.21mm,增長率達到57%,同日出現(xiàn)52.2mm 的降雨,降雨后又迅速減至29.89mm.
據(jù)此可見,利用CORS 站接收的GNSS 數(shù)據(jù)解譯的大氣水汽結果與實際情況相符,并且能較好地反應大氣水汽與降雨之間存在的響應關系.
為進一步完善CORS 站綜合監(jiān)測能力,文章將CORS 站測量地表位置變化的基本用途進行實驗并與降雨數(shù)據(jù)結合分析獲得結果.圖12 為利用GAMIT基線解算獲得的BFQE 測站2021 年在北(north,N)、東(east,E)、天頂(up,U)方向上的形變數(shù)據(jù)以及形變速率圖.由圖12(a)可知,CORS 站BFQE 形變量在2021 年全年時間內較為穩(wěn)定,N、E、U 方向全年累計形變量分別為-12mm、9mm、7mm.其中,N 方向和E 方向在全年形變波動較小,基本均在±0.01m內,而U 方向形變波動相比于N,E 方向較大,在±0.02m 內波動,且在夏季波動尤為明顯.為進一步分析其波動時段,文章計算了測站U 方向的形變速率,得到形變速率圖,圖中紅色曲線為擬合得到的形變速率趨勢項,可以看出在年積日第131—275 天,形變速率呈顯為上凸的趨勢,明顯大于其余時段.
圖12 2021年測站形變反演
雨水入滲過程會直接導致土壤含水率的變化,從而影響土地土質的強度變化,直接導致了地表穩(wěn)定性的變化.因此,地表形變與降雨過程中的降雨量關系密切.圖13 為結合降雨數(shù)據(jù)繪制的CORS 站BFQE在U 方向上的形變速率圖.在第172 天及220 天,降雨量分別達到峰值78.2mm 及次峰值52.2mm,形變速率趨勢在次峰值降雨的后一天達到頂點.
圖13 形變速率與降雨關系圖
文章統(tǒng)計了形變速率位于前15%的天數(shù)分布,均在降雨豐富的第131—275 天,同時值得注意的是在第107—133 天時間段,無降雨情況,地表形變速率表現(xiàn)為0.03m/d,遠低于降雨時段.由此可得出降雨是影響測站周圍形變量改變的一個重要因素.
CORS 站作為一種地面GPS 接收機,能夠快速接收衛(wèi)星數(shù)據(jù),本文將CORS 站接收數(shù)據(jù)與GNSS技術的有效解譯能力相結合,用于地表環(huán)境參數(shù)的解譯.通過利用GNSS-IR 技術解譯得到的土壤濕度和積雪深度,與實測值的相關系數(shù)均表現(xiàn)出強相關性,分別達到0.89 和0.90,并且在降雨天氣出現(xiàn)時呈現(xiàn)上升趨勢.利用GNSS-PWV 技術解譯得到的大氣水汽含量的趨勢與實際情況相符,始終在降雨之前開始積累,在降雨后開始下降,并且與降雨之間存在良好的響應關系.此外,地表形變在U 方向上跟隨降雨量及土壤濕度的變化而出現(xiàn)波動.
綜上所述,將CORS 站用于地表環(huán)境參數(shù)綜合監(jiān)測,既能有效解譯各參數(shù)值,獲得其變化趨勢,同時能較好地反應各參數(shù)之間的響應關系.這種綜合分析的方法為研究人員提供了一種便捷而有效的途徑,有助于更全面地理解地表環(huán)境參數(shù)之間的關聯(lián),進一步推進相關領域的研究工作.