張勤耘，金雙根

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

0 引言

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

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

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

1 理論和方法

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

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

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

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

由圖4可知：S2的RMSE為0.683 5，S5的RMSE為0.393 7；在火山爆發(fā)前，SNR差異值穩(wěn)定在[—1，1]這個區(qū)間；在火山爆發(fā)后，SNR異常差異值振蕩現(xiàn)象大，其中對應(yīng)于S2的最大差異絕對值為3 dB，對應(yīng)于S5的最大差異絕對值為1.9 dB。

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

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

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

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

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

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

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

3 結(jié)束語

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

[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，F(xiàn)REYMUELLER 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.