GNSS信噪比在火山灰監測中的應用

張勤耘，金雙根

(1.上海大學 通信與信息工程學院，上海 200444；2.中國科學院 上海天文臺，上海 200030)

0 引言

火山噴發時，會噴射出大量的巖石塊、火山灰、火山氣體等物質。這些物質會嚴重影響建筑物的完整性，破壞陸地或空中交通安全，扼制植物的生長，影響人類的健康生活和安全。因此監測并預報火山灰十分重要。世界各地的火山活動可通過各種各樣的地基或空基設備進行監測，但這2類設備各有優勢與缺點：地基設備可以提供相對持續的高精度時空覆蓋度，但是密度不足夠大、覆蓋度有限；相反地，衛星系統可以確保測量的全球空間密度以及所研究區域的概要圖，但是空間覆蓋率較低，精度不高。為了更好地降低火山煙流的危害，相關的監測預報系統和設備如意大利地震臺網中心(Istituto Nazionale di Geofisica e Vulcanologia，INGV)、激光雷達(light detection and ranging，LiDAR)、多普勒雷達(doppler radar，DR)等已經逐漸得以完善[1-3]。除此之外，火山煙流也可能影響全球衛星導航系統(global navigation satellite system，GNSS)的衛星信號。

文獻[4]利用無電離層延遲測量(ionosphere-free linear combination，LC)的異常，采用地震層析成像算法將折射率異常和推算的火山灰的溫度映射在時空圖中，分析了2000-08的日本三宅島火山(Miyake Jima)相關情況，表明熱火山灰造成了全球定位系統(global position system，GPS)信號路徑的延遲[4]。文獻[5]采用穩健型統計學方法，以2007-09-04和2007-09-05的意大利埃特納火山(Etna)為例，探索大氣中火山煙流是否會影響GPS非差擬合相位殘差(post-fit phase residuals，PPR)，研究結果表明，火山煙流明顯影響了GPS PPR值，且是日本三宅島火山所得殘差值的1/10，猜測殘差值與火山爆發強度存在一定的相關性[5]。文獻[6]以運動學定位方法為基礎，以2009年阿拉斯加州的里道特火山(Redoubt)為例，獲取全球衛星導航系統(global navigation satellite system，GNSS)數據的相位殘差，將大規模、短周期的定位偏移與穿過稠密的火山煙流的“衛星-測站對”的相位延遲聯系起來，并提供了一種測量火山煙流擴散方向的方法即沿著衛星在天空圖中的運行軌跡表示出相位殘差，由此猜想利用衛星與地表測站之間的GPS信號的相位殘差可以監測火山，同時也指出由于火山煙流的濃度可能不足夠高而不能對衛星信號造成顯著影響，即這種方法并不適用于所有的火山煙流[6]。上述方法都是使用標準的地學定位處理技術，并將火山灰信號作為一種非典型的大氣效應，需要分析GPS載波相位觀測范圍，還要評估站心坐標、軌道、時鐘、大氣延遲以及相位模糊度等影響。而且火山灰與GPS接收機隨時間變化，“衛星-測站對”的路徑并不一定就通過火山灰；需要利用最小二乘估計法解算，其他模型誤差都會影響待估參數。文獻[7]提出利用GPS信號信噪比數據(signal-to-noise ratio，SNR)估計2008年阿拉斯加州奧克莫克火山(Okmok)和2009年里道特火山(Redoubt)，研究表明在火山爆發后的一小段時間內，SNR會在有一定的突降之后恢復到之前的變化趨勢上[7]。文獻[8]使用了非差擬合相位殘差方法PPR和SNR探測方法監測2012-08-06爆發的新西蘭湯加里羅火山(Tongariro)，進一步表明使用GPS數據監測火山灰的可能性，有助于提高對火山的認識[8]。

但以上研究仍存在諸多問題，如不是所有GPS SNR觀測都能監測火山灰[9-10]、信噪比受電離層閃爍的影響，以及火山灰傳播速度不清楚。本文利用GPS SNR數據監測2016-10-08的日本阿蘇山火山灰變化，給出GPS SNR異常值探測火山煙流的結果，并估計火山灰平均傳播速度。

1 理論和方法

GPS信號是廣大用戶發送的用于導航定位的調制波，它包含處于L波段的載波、測距碼、導航電文。常用載波有2種：L1和L2波段的載波頻率分別為fL1=154f0=1 575.42 MHz、fL2=120f0=1 227.60 MHz；對應的波長為λ1=19.032 cm、λ2=24.420 cm。現代化后，又增加了L5載波，頻率為fL5=115f0=1 176.45 MHz、λ5=25.480 cm。測距碼包括C/A碼和P碼，前者只能調制在L1上，碼率為1.023 MHz，后者的碼率為10.23 MHz。

另外，GPS接收機通常記錄SNR數據，這些數據是載波與噪聲密度比(C/N0)，即信號功率與噪聲功率的比值。SNR與C/N0之間的關系可以表示為SNR=(C/N0)/B，以分貝(dB)為單位；其中B為帶寬。相對應載波L1和L2，信噪比類型主要有S1和S2，此外三頻觀測還有S5等。SNR數據不提供衛星發射機與接收機之間的距離信息，因此對定位解決方案沒有直接貢獻，也是由于這個原因，SNR數據通常被地理物理學家和大地測量工作者忽略。然而SNR數據可以直接測得信號阻塞。當GPS信號經過火山灰時，部分信號會發生衰減或者散射，這就意味著到達GPS接收機的信號功率會降低。相反地，大氣中的水汽并不直接影響信號功率[7]。這種對水不敏感的特性正是建立L波段雷達觀測灰云的原因。PPR方法既受水汽、溫度的影響，同時也受火山灰的影響；而SNR僅受火山灰的影響：這也是利用GPS信噪比SNR探測火山灰的優點之一[10]。

信噪比SNR數據值非常依賴高度角。當衛星高度角從5°到65°變化時，SNR的值會緩慢增加，這主要是歸功于天線增益的類型。L1 SNR數據在高度角大于25°時有高頻噪聲(系統噪聲或隨機噪聲)。相反地，海拔更高的L2 SNR數據的高頻噪聲的級別更小。這種差異是由于C/A碼比L2P或L2C短而更易遭受交流干擾。然而，L2P卻具有比C/A或L2C更低的SNR值，這是因為接收機在檢索過程中不能使用加密碼。注意只有在2005年以后發射的衛星上才可以獲得L2C碼。此外除非用戶要求記錄L2C碼，否則不會記錄。

當高度角小于30°時，SNR數據的振蕩是由地面反射導致的。反射信號比直接信號經過更長的距離并且與直接信號進行干涉，這將導致需要觀測調制。高度角低于30°的SNR數據的頻率與土壤濕度、雪深和海平面高度相關。這里，只考慮具有更高的高度角的SNR數據作為探測火山灰的興趣因子。這主要是因為高度角高的L2P SNR數據具有更好的精確度，所以用它去分離火山灰的影響。

在利用GPS信噪比SNR數據探測火山煙流時，首先要對SNR數據進行預處理，舍棄高度角低于30°的數據，因為低于30°的SNR數據極容易受到地面反射的影響，與土壤濕度、雪深、海平面高度等因素相關[11-15]。通過多項式擬合得到擬合信噪比值，并與實際信噪比值相減，得到信噪比差異值。分析SNR異常值隨時間變化的情況，進而利用SNR數據監測火山灰變化。此外，還要注意可以取火山爆發前30 min的SNR數據的標準差作為性能控制尺度，探測值需要是標準差的2.5倍；若探測值不滿足這個條件，則該衛星記錄的SNR數據不可用。

此外，在條件允許的情況下，需要驗證信噪比是否受電離層閃爍的影響，即需要一些接收機數據作為對比，以證明SNR數據的改變不是由電離層閃爍產生的；要想更充分地分析利用SNR數據監測火山灰的優劣性，還需要對比雷達掃描結果，以確定灰塵粒子尺寸和密度對于使用SNR數據分析火山灰的限制。

實際上火山爆發時，噴發出的物質除了呈固態小顆粒的火山灰，還有含有水蒸氣和硫化物的火山氣體。在海拔較高的地方，火山灰和水蒸氣會分離。這是因為在海拔較高的位置，水蒸氣依然可以保持充足的密度和體積，但是火山灰卻會沉降并向水平方向的擴散，這也導致在水平方向和垂直方向上離火山越遠火山灰密度越低；所以距離火山較近的測站受火山灰影響大，而其他距離較遠的測站受火山灰影響微弱，即火山灰密度在信噪比SNR的衰減現象中扮演著十分重要的角色。

2 結果與分析

阿蘇山(Mount Aso)是日本最大的活火山，也是世界大型的活火山之一，位于九州島熊本縣東北部，最高點海拔達到1 592 m。2016-10-08，當地時間凌晨1：46，日本氣象廳(Japan Meterological Agency，JMA)報道阿蘇火山中岳火山口發生爆發式噴發。由于當日環境比較惡劣，影響了視覺觀察，但是衛星影像顯示火山爆發當日火山灰到達了11 km的高度。火山沉積物向西北側面擴散了1.6 km，向東南側擴散了1 km，在東北側也比較多。

我們利用的是日本國土地理院(Geographical Information Authority of Japan，GSI)所建立的GNSS地球觀測網(GNSS Earth observation network system，GEONET)數據，這些GNSS數據的采樣間隔為30或1 s，采樣間隔依賴于接收機所處位置的環境、GNSS接收機的類型和遙測系統。根據文獻[6]的研究結果，15到30 s的標準采樣率就足以分析經過煙流的信號；本次研究所使用的數據間隔是30 s。

本文選取阿蘇山火山30 km范圍以內的測站，其測站位置與阿蘇山的關系如圖1所示，其中較粗曲線部分表示第281天16：00到18：00的衛星軌跡，較細曲線部分表示觀測窗口全時段內的衛星軌跡。

針對各個測站提取火山爆發時刻前后16：00—18：00(UTC)時間段內的GNSS SNR數據，分別分析SNR與時間以及衛星高度角之間的變化關系，探索SNR與時間以及衛星高度角之間的相關性，并通過二項式擬合得到對應于各個測站的SNR異常值，然后統計信噪比異常發生的概率。經過對上述9個測站的試驗與分析，最終呈現PRN32衛星與阿蘇山10 km范圍以內的3個測站即0701、0703、0704之間的“測站-衛星對”的結果，其中0703測站的情況如圖2所示。

從圖2可知：對于0703測站，第279天、第280天、第282天、第283天所對應的信噪比SNR隨高度角的增大而逐漸增大，而變化趨勢則逐漸減小，沒有任何變化異常；而對應于第281天在衛星高度角40°左右時的信噪比SNR明顯下降3 dB左右(而對于0701測站和0704測站，第279天、第280天、第281天、第282天、第283天，SNR沒有任何異常變化現象，且各天SNR的變化趨勢同一，變化曲線基本同一)。故得到初步結果：火山灰僅對0703-PRN32的“測站-衛星對”有影響，這與火山灰的擴散方向和衛星測站對之間的路徑是否穿過火山灰等因素有關。

由于接收機觀測數據中并未記錄到對應于L1的信噪比數據S1，但記錄了與L2和L5對應的信噪比數據S2和S5；所以對于有明顯下降現象的測站0703，除了考慮上述S2外，又考慮了S5(見圖3)。綜合圖2與圖3可以看出，同0703測站的S2數據相比，整體來看S5的值要比S2的值高15 dB左右，而且在高度角40°左右時，第281天的S5也有相對明顯的異常下降現象。

為進一步研究火山灰與信噪比異常變化之間的關系，可以利用4階多項式擬合得到期望SNR值，并與實際SNR值相減，得到SNR差異值。如圖4為第281天針對0703測站的PRN32衛星的觀測值信噪比差異值隨高度角變化的情況。

由圖4可知：S2的RMSE為0.683 5，S5的RMSE為0.393 7；在火山爆發前，SNR差異值穩定在[—1，1]這個區間；在火山爆發后，SNR異常差異值振蕩現象大，其中對應于S2的最大差異絕對值為3 dB，對應于S5的最大差異絕對值為1.9 dB。

為進一步研究信噪比與火山灰之間的關系，可以從“信噪比的下降現象是否在火山爆發后發生”這一方面考慮(如圖5所示)。

由圖5可知，信噪比異常下降現象的確發生在火山爆發時刻之后2 min左右，而且異常變化間隔大概在3 min左右，變化間隔較短。這可能一方面是因為在火山爆發瞬間噴發出大量高濃度的火山煙流，而煙流快速上升的同時火山灰也逐漸向側方向擴散，火山煙流的濃度降低；另一方面因為衛星與測站的相對位置也在發生變化，衛星的直射信號路徑隨著時間可能穿過濃度較低的火山煙流或者不穿過火山煙流。

經過上述針對各測站的衛星高度角和時間與信噪比之間的變化關系的“控制變量”對比研究分析，即針對同一測站同一衛星研究不同天的SNR變化分析，以及針對不同測站同一衛星同一天的SNR變化分析，排除了由于測站自身因素、測站周邊地況、衛星自身因素以及其他天氣情況等因素，得到0703測站所對應的PRN32衛星在第281天(火山爆發當天)的SNR異常的確是由于火山灰引起。

通過大量數據綜合考慮0703測站所有衛星的信噪比差異值(實際信噪比值與擬合信噪比值相減得到)，并統計各信噪比差異值的概率情況如圖6所示。

由圖6可知：信噪比差異值出現的概率小于0.12，而差異值在[—3，—1.75]區間內的概率小于0.02；故可以判斷信噪比異常現象的出現概率是極其小的，若沒有外在因素的影響，信噪比異常現象不會出現，即火山灰的確對GPS信號造成了影響。阿蘇山周邊9個測站對應的PRN32信噪比SNR2的時間序列變化分析結果如圖7所示。

由圖7可知：對于PRN32衛星，整體來看，除了0703測站對應的信噪比值在16：50左右有一個突然的降低之外，各測站的信噪比SNR沒有驟降而是隨時間的增長而逐漸降低而且降低的幅度越來越大；在17：20之前，除了0703測站對應的信噪比差異值在16：50左右最大驟降幅度達到了3 dB Hz，其他測站對應的信噪比差異值絕大多數在[—1，1]之間；而在17：20之后，所有測站的信噪比差異值振蕩幅度增大，而考慮到在高度角小于30°時測站接收機受多路徑效應的影響，這種大幅度振蕩現象也可以得到解釋。圖8再次驗證在“衛星-測站對”的信號路徑未經過火山灰時，信噪比差異值絕對值小于1 dB，也從另一個角度反映了測站到火山的距離并不是唯一影響信噪比差異值的因素。

綜合火山爆發后火山灰到達衛星測站之間的連線(下文將此相交點稱為交點)所經過的距離，以及火山爆發到信噪比發生變化時所經過的時間，大致可以估計出火山灰傳播的平均速度。當信噪比發生變化時，已知變化時間dt，衛星高度角ele，測站與火山口之間的距離s，估計的交點到火山口之間連線與鉛垂線之間的夾角β(β可以是一個區間角)，可以求出衛星-測站連線與交點-火山口連線之間的夾角α，最后可以根據正弦定理求得交點與火山口之間的距離ρ，進而求出火山灰傳播的平均速度v(30 m/s左右)。

3 結束語

本文利用GEONET的SNR數據探測了2016-10-08日本阿蘇山火山灰變化，驗證利用SNR數據的變化監測火山灰的可行性。在阿蘇山火山爆發期間，火山北側的0703測站所接收到的PRN32衛星的SNR數據受到較大影響，信噪比的值在較短時間內會有明顯的降低，在16：50左右最大驟降幅度達到了3 dB；其他衛星幾乎沒有影響，而其他附近測站SNR數據幾乎也沒有明顯異常現象，主要是這些測站與衛星之間信號沒有穿過一定濃度的火山灰。根據火山灰到達衛星測站之間的距離，以及火山爆發到信噪比發生變化時所經過的時間，可以估計出火山灰傳播平均速度約為30 m/s。將來可以利用格洛納斯衛星導航系統(global navigation satellite system,GLONASS)、伽利略衛星導航系統(Galileo satellite navigation system，Galileo)和中國北斗衛星導航系統(BeiDou navigation satellite system，BDS)的SNR數據對火山灰進行探測，并與GPS數據的結果進行對比，進一步研究利用GNSS SNR數據探測火山灰的可行性。此外，還可以進一步考慮火山煙流的成分對GNSS SNR數據的具體影響。

[1] SCOLLO S，PRESTIFILIPPO M，SPATA G，et al.Monitoring and forecasting Etna volcanic plumes[J].Natural Hazards & Earth System Sciences，2009，9(5)：11993.

[2] SCOLLO S，PRESTIFILIPPO M，PECORA E，et al.Height estimation of the 2011—2013 Etna lava fountains[J].Annals of Geophysics，2014，57(2)：1-25.

[3] GOUHIER M，HARRIS A，CALVARI S，et al.Erratum to：Lava discharge during Etna’s January 2011 fire fountain tracked using MSG-SEVIRI[J].Bulletin of Volcanology，2012，74(5)：1261-1261.

[4] HOULIé N，BRIOLE P，NERCESSIAN A，et al.Sounding the plume of the 18 August 2000 eruption of Miyakejima volcano(Japan)using GPS[J].Nature，2005，414(5)：583-583.

[5] ARANZULLA M，CANNAVF，SCOLLO S，et al.Volcanic ash detection by GPS signal[J].GPS Solutions，2013，17(4)：485-497.

[6] GRAPENTHIN R，FREYMUELLER J T，KAUFMAN A M.Geodetic observations during the 2009 eruption of Redoubt Volcano，Alaska[J].Journal of Volcanology & Geothermal Research，2013，259(6)：115-132.

[7] LARSON K M.A new way to detect volcanic plumes[J].Geophysical Research Letters，2013，40(11)：2657-2660.

[8] 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[J].Journal of Volcanology & Geothermal Research，2014(286)：387-396.

[9] ARANZULLA M，CANNAVO F，Scollo S.Detection of Volcanic Plumes by GPS：the 23 November 2013 Episode on Mt.Etna[J].Annals of geophysics(Annali di geofisica)，2014，57(Fast Track 2).

[10] 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[J].Earth，Planets and Space，2015，67(1)：157.

[11] LARSON K M，BRAUN J J，SMALL E E，et al.GPS multipath and its relation to near-surface soil moisture content[J].IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing，2010，3(1)：91-99.

[12] NIEVINSKI F G，LARSON K M.Inverse modeling of GPS multipath for snow depth estimation—part I：formulation and simulations[J].IEEE Transactions on Geoscience & Remote Sensing，2014，52(10)：6555-6563.

[13] NIEVINSKI F G，LARSON K M.Inverse modeling of GPS multipath for snow depth estimation—part II：application and validation[J].IEEE Transactions on Geoscience & Remote Sensing，2014，52(10)：6564-6573.

[14] JIN S G，NAJIBI N.Sensing snow height and surface temperature variations in Greenland from GPS reflected signals[J].Advances in Space Research，2014，53(11)：1623-1633.

[15] JIN S G，QIAN X D，KUTOGLU H.Snow depth variations estimated from GPS-reflectometry：a case study in Alaska from L2P SNR data[J].Remote Sensing，2016，8(1)：63.