1 Hz GNSS電離層相位閃爍因子提取及在北極區(qū)域的驗(yàn)證

趙東升，李 旺，李宸棟，唐 旭，張克非

1. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116； 2. 寧波諾丁漢大學(xué)理工學(xué)院，浙江 寧波 315100； 3. 南京信息工程大學(xué)遙感與測(cè)繪學(xué)院，江蘇 南京 210044

全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)已經(jīng)成為監(jiān)測(cè)電離層的重要手段[1-4]。電離層閃爍是電離層中的不規(guī)則體，是近地磁赤道和地磁極區(qū)域頻發(fā)的一種天文災(zāi)害[5]，會(huì)對(duì)穿越其中的GNSS信號(hào)造成衍射和散射，進(jìn)而使得GNSS信號(hào)的振幅或者相位劇烈波動(dòng)，降低GNSS信號(hào)的信噪比甚至失鎖，給GNSS系統(tǒng)位置、導(dǎo)航與時(shí)間(PNT)服務(wù)的穩(wěn)定性帶來(lái)巨大挑戰(zhàn)[6-8]。電離層閃爍因子是電離層閃爍對(duì)GNSS信號(hào)影響的定量表征，可以反映出電離層閃爍的強(qiáng)弱，是實(shí)現(xiàn)電離層閃爍的監(jiān)測(cè)、建模、預(yù)報(bào)的基礎(chǔ)數(shù)據(jù)[9]，對(duì)于科學(xué)認(rèn)識(shí)電離層閃爍發(fā)生規(guī)律、修正或減弱電離層閃爍對(duì)GNSS的不利影響進(jìn)而實(shí)現(xiàn)高精度定位具有重要意義。

為了測(cè)量電離層閃爍對(duì)GNSS信號(hào)的影響，通常需要電離層閃爍監(jiān)測(cè)接收機(jī)(ionospheric scintillation monitoring receivers，ISMR)，它可以直接給出兩類閃爍因子：振幅閃爍因子(S4)和相位閃爍因子(σφ)。由于ISMR運(yùn)行在50 Hz的采樣頻率，使得ISMR可以直接跟蹤到由于電離層閃爍而產(chǎn)生的信號(hào)振幅和相位的變換，但是這也導(dǎo)致其需要較大的存儲(chǔ)空間，且價(jià)格較貴，限制了ISMR的布站數(shù)量。當(dāng)前布設(shè)ISMR的電離層閃爍監(jiān)測(cè)網(wǎng)絡(luò)主要有歐洲航天局的Monitor項(xiàng)目、巴西的CIGALA/CALIBRA、ICEA和LISN項(xiàng)目和加拿大高緯度北極電離層閃爍監(jiān)測(cè)網(wǎng)絡(luò)(CHAIN)，共約100個(gè)測(cè)站。相比于ISMR，采樣頻率通常在1 Hz及以下的接收機(jī)分布十分廣泛，基本可以實(shí)現(xiàn)對(duì)陸地及近海的全覆蓋，且這些測(cè)站的設(shè)立時(shí)間更久，可以提供更長(zhǎng)時(shí)間序列的觀測(cè)數(shù)據(jù)，有助于研究電離層閃爍變化的長(zhǎng)周期項(xiàng)，提高電離層閃爍監(jiān)測(cè)及建模的精度[10]，但在1 Hz及以下采樣頻率的觀測(cè)數(shù)據(jù)中，由于觀測(cè)數(shù)據(jù)時(shí)間分辨率的降低，電離層閃爍誤差極容易被其他誤差(尤其是接收機(jī)鐘差)所淹沒，為此需要專門的方法構(gòu)建基于低采樣頻率(1 Hz及以下)GNSS觀測(cè)數(shù)據(jù)的電離層閃爍因子。

地磁北極區(qū)域的電離層閃爍以相位閃爍為主[11]，因而本文將以相位閃爍因子作為主要研究對(duì)象。為構(gòu)建基于低采樣頻率觀測(cè)數(shù)據(jù)的相位閃爍因子，文獻(xiàn)[12]基于無(wú)幾何組合觀測(cè)值提出了總電子含量(TEC)變化率指數(shù)(ROTI)，可以較好減弱接收機(jī)鐘差的影響，已經(jīng)成為電離層閃爍監(jiān)測(cè)領(lǐng)域應(yīng)用最為廣泛的標(biāo)準(zhǔn)閃爍因子之一[12-14]，ROTI可以給出每個(gè)觀測(cè)歷元的閃爍因子。在此基礎(chǔ)上，文獻(xiàn)[15]通過對(duì)一段時(shí)間內(nèi)所有可視衛(wèi)星獲得的TEC值以衛(wèi)星高度角定權(quán)的方式進(jìn)行加權(quán)平均從而提出了另一種閃爍因子——TEC起伏(AATR)，它可以用來(lái)預(yù)報(bào)發(fā)生閃爍的時(shí)間段和區(qū)域，已經(jīng)被歐洲地球靜止導(dǎo)航重疊服務(wù)(EGNOS)選為表征電離層活動(dòng)的參數(shù)之一[14]。文獻(xiàn)[16]將無(wú)電離層組合殘差的標(biāo)準(zhǔn)差作為監(jiān)測(cè)電離層閃爍的因子并提出了對(duì)應(yīng)的修正無(wú)電離層組合中的接收機(jī)鐘差方法。文獻(xiàn)[17]利用1 Hz的斜路徑總電子含量和小波變換技術(shù)提出了一種類相位閃爍因子，該因子被認(rèn)為可以有效地預(yù)警歐洲部分北極區(qū)域的電離層閃爍事件。上述研究所構(gòu)建的閃爍因子均基于GNSS組合信號(hào)，其所獲得的閃爍因子為電離層閃爍對(duì)組合信號(hào)作用的結(jié)果，無(wú)法用于研究電離層閃爍對(duì)每一頻點(diǎn)信號(hào)的影響。

為構(gòu)建基于1 Hz GNSS載波相位觀測(cè)值的電離層相位閃爍因子，首先，本文將利用大地測(cè)量趨勢(shì)分離(geodetic detrending)、精密單點(diǎn)定位、無(wú)電離層組合、周跳探測(cè)與修復(fù)等技術(shù)剔除載波相位觀測(cè)值中的幾何距離以及固體潮、海潮、大氣潮、極潮、天線相位中心、相對(duì)論、相位纏繞、接收機(jī)鐘差和對(duì)流層延遲誤差，獲得載波相位觀測(cè)值的殘差；然后利用離散小波變換技術(shù)降低噪聲對(duì)電離層閃爍信號(hào)進(jìn)行提??；最后利用連續(xù)小波變換技術(shù)對(duì)去噪后的信號(hào)進(jìn)行時(shí)頻分析，對(duì)經(jīng)驗(yàn)頻帶內(nèi)的小波系數(shù)進(jìn)行逆變換，提取出電離層閃爍信號(hào)，并設(shè)定滑動(dòng)窗口，對(duì)窗口內(nèi)的電離層閃爍信號(hào)取標(biāo)準(zhǔn)差，從而完成每個(gè)頻點(diǎn)上電離層相位閃爍因子的構(gòu)建。本文將基于CHAIN提供的大量觀測(cè)數(shù)據(jù)，通過對(duì)比分析提出的電離層相位閃爍因子、ROTI與ISMR提供的相位閃爍因子的相關(guān)性和準(zhǔn)確性，驗(yàn)證其有效性。

1 閃爍因子的構(gòu)建方法

本文基于1 Hz GNSS觀測(cè)數(shù)據(jù)，致力于構(gòu)建可以提供每個(gè)頻點(diǎn)的電離層相位閃爍因子，詳細(xì)過程如下。

1.1 大地測(cè)量趨勢(shì)分離

大地測(cè)量趨勢(shì)分離的概念為利用公式、模型及公開參數(shù)等手段剔除GNSS載波相位中較為容易估計(jì)的誤差[16]。為提取GNSS載波相位觀測(cè)值中的電離層閃爍誤差，本文采用了如下方法進(jìn)行大地測(cè)量趨勢(shì)分離。

(1) 站星幾何距離改正。為改正每顆衛(wèi)星至測(cè)站之間的幾何距離，需準(zhǔn)確確定測(cè)站坐標(biāo)和每個(gè)觀測(cè)歷元的衛(wèi)星坐標(biāo)。測(cè)站坐標(biāo)利用CSRS-PPP在線精密單點(diǎn)定位軟件通過靜態(tài)解算方法獲得。衛(wèi)星坐標(biāo)通過國(guó)際GNSS服務(wù)組織(IGS)提供的精密星歷文件并采用二階拉格朗日內(nèi)插算法獲得。根據(jù)測(cè)站坐標(biāo)和每個(gè)觀測(cè)歷元的衛(wèi)星坐標(biāo)，通過歐幾里得度量獲得站星幾何距離改正參數(shù)。

(2) 固體潮誤差采用二階簡(jiǎn)化潮汐模型進(jìn)行改正，海潮、大氣潮和極潮誤差均采用國(guó)際地球自轉(zhuǎn)服務(wù)(IERS)推薦的慣用模型改正[18]。

(3) 測(cè)站的接收機(jī)天線相位中心(包括接收機(jī)天線參考點(diǎn)、接收機(jī)天線平均相位中心和接收機(jī)天線瞬時(shí)相位中心)和衛(wèi)星天線相位中心(包括衛(wèi)星天線平均相位中心和衛(wèi)星天線瞬時(shí)相位中心)利用IGS提供的天線相位中心模型(即igs14.atx)改正。

(4) 衛(wèi)星鐘差利用IGS提供的精密鐘差文件進(jìn)行改正。

(5) 相對(duì)論效應(yīng)和相位纏繞誤差采用文獻(xiàn)[19]所述相關(guān)方法予以改正。

(6) 對(duì)流層延遲分為干分量和濕分量?jī)刹糠?，其中干分量采用文獻(xiàn)[19]所述相關(guān)方法予以改正。由于難以利用模型對(duì)濕分量進(jìn)行精密估計(jì)，本文采用如下較為簡(jiǎn)單的方法估計(jì)：將濕分量的初值設(shè)為0.1 m，余下歷元的濕分量作為隨機(jī)游走參數(shù)進(jìn)行估計(jì)，對(duì)于未改正的對(duì)流層延遲濕分量作為精密單點(diǎn)定位中的參數(shù)進(jìn)行解算。

(1)

1.2 對(duì)流層濕延遲未能通過模型改正部分和測(cè)站接收機(jī)鐘差的改正

VIF=Ax-(LIF-DIF)

(2)

(3)

式中，VIF為殘差；A為系數(shù)矩陣；DIF為無(wú)電離層組合模型改正值；x為待解算參數(shù)包括測(cè)站坐標(biāo)(δX，δY，δZ)、天頂對(duì)流層濕延遲(δρz,wc)、接收機(jī)鐘差(cδtR)和同觀測(cè)歷元的m顆可視衛(wèi)星的模糊度(N)。該觀測(cè)方程采用卡爾曼濾波靜態(tài)解算方式進(jìn)行求解，濾波中的狀態(tài)轉(zhuǎn)移矩陣(Φ)為

Φ=diag[111011…1]

(4)

系統(tǒng)噪聲向量的協(xié)因數(shù)矩陣(Q)為

Q=diag[0009×101010-40…0]

(5)

為提高天頂對(duì)流層濕延遲和測(cè)站接收機(jī)鐘差的估計(jì)精度，本文采用正向和反向運(yùn)算相結(jié)合的方式，并將反向運(yùn)算的結(jié)果作為未能利用模型改正的天頂對(duì)流層濕延遲(δρz,wc)和測(cè)站接收機(jī)鐘差(c·δtR)的估計(jì)結(jié)果。值得指出的是在估計(jì)接收機(jī)鐘差和對(duì)流層濕延遲時(shí)，也可采用非組合精密單點(diǎn)定位技術(shù)[20]。利用Niell濕延遲投影函數(shù)，可將天頂對(duì)流層延遲(δρz,wc)換算到每一衛(wèi)星斜路徑上，進(jìn)而改正該衛(wèi)星觀測(cè)值中的對(duì)流層延遲誤差。如果接收機(jī)鐘比較穩(wěn)定，不存在鐘跳時(shí)，利用PPP方法獲得的接收機(jī)鐘差可直接用來(lái)改正觀測(cè)值中的鐘差；但對(duì)于大多數(shù)測(cè)地型接收機(jī)，接收機(jī)廠商會(huì)在鐘差漂移到某一閾值時(shí)，通過對(duì)其中插入鐘跳的方式，使得接收機(jī)內(nèi)部時(shí)鐘與衛(wèi)星鐘同步精度控制在一定范圍[21]。此時(shí)利用PPP方法無(wú)法準(zhǔn)確估計(jì)出接收機(jī)鐘差變化的細(xì)節(jié)部分，需要對(duì)鐘差進(jìn)行進(jìn)一步精細(xì)估計(jì)。

1.3 測(cè)站接收機(jī)鐘中存在鐘跳時(shí)的鐘差精細(xì)估計(jì)及周跳探測(cè)與修復(fù)

GNSS載波相位觀測(cè)值經(jīng)過上文所提出方法修正后，其殘差為

(6)

式中，δtR,cj表示未能經(jīng)過PPP方法修正的鐘差。為提取該殘差中的鐘差，本文利用兩個(gè)頻點(diǎn)的殘差組成無(wú)電離層組合觀測(cè)值，可表示為

(7)

該組合可以較為完好地消除電離層折射作用的影響，但無(wú)法消除電離層衍射與散射作用(即電離層閃爍)和模糊度的影響。由于模糊度在同一觀測(cè)弧段內(nèi)為同一值，因此可以通過歷元間做差的方法消除。忽略觀測(cè)噪聲的影響，歷元間做差之后的殘差可表示如下

(8)

為避免誤讀、誤會(huì)，釋義，便成了首要功課。“三雜”之雜，為多種多樣、豐富多彩意，決非雜亂之雜……大凡雜文家，一定會(huì)有相應(yīng)的知識(shí)儲(chǔ)備、文體實(shí)踐和人生閱歷——如此之“三”以“三雜”謂之，不算離奇、搞怪吧？

由式(8)可以看出，周跳會(huì)影響鐘差的估計(jì)精度，故而在估計(jì)鐘差之前需要進(jìn)行周跳探測(cè)與修復(fù)。本文首先利用式(6)給出的兩個(gè)頻點(diǎn)的殘差組成雙頻寬項(xiàng)組合和無(wú)幾何關(guān)系組合聯(lián)合進(jìn)行初次周跳探測(cè)，該方法可以修復(fù)大多數(shù)的周跳，但通常情況下，尤其是在電離層較為活躍的條件下，難以準(zhǔn)確修復(fù)全部周跳。為此需要對(duì)未準(zhǔn)確修復(fù)的周跳進(jìn)行進(jìn)一步探測(cè)和修復(fù)。本文將初步修正周跳后的每個(gè)頻點(diǎn)的殘差組成的無(wú)電離層組合觀測(cè)值并進(jìn)行歷元間差分，由于該組合觀測(cè)值可以極大地降低電離層的影響，從而提供精度較高的檢測(cè)值，能夠較為準(zhǔn)確的檢測(cè)出小周跳的存在。但只采用這一組合，無(wú)法確定周跳的大小，導(dǎo)致無(wú)法修復(fù)探測(cè)到的周跳。為此本文將探測(cè)發(fā)生周跳的歷元作為新弧段的開始，從而在該觀測(cè)弧段內(nèi)可以為鐘差的估計(jì)提供無(wú)周跳影響的歷元間差分無(wú)電離層組合觀測(cè)值。需要指出的是單獨(dú)用無(wú)電離層組合觀測(cè)值進(jìn)行周跳探測(cè)會(huì)存在不敏感周跳組問題[22]，仍然難以探測(cè)全部周跳，但因?yàn)楸疚牟捎昧诵〔ㄗ儞Q技術(shù)，少量周跳不會(huì)影響最終閃爍因子的估計(jì)。

為避免多路徑效應(yīng)對(duì)估計(jì)接收機(jī)鐘差的影響，本文設(shè)置衛(wèi)星的截止高度角為20°。為綜合利用每個(gè)歷元內(nèi)截止高度角以上的所有衛(wèi)星觀測(cè)值，本文采用衛(wèi)星高度角定權(quán)的方法對(duì)各衛(wèi)星的歷元間差分無(wú)電離層組合觀測(cè)值進(jìn)行加權(quán)平均，其權(quán)重(P)為

P=sin2E

(9)

加權(quán)平均后的歷元間差分后的鐘差可表示為

(10)

式中，m表示歷元k中衛(wèi)星的數(shù)量；s表示衛(wèi)星的序號(hào)。歷元k時(shí)刻的鐘差可以表示為從初始?xì)v元至歷元k的歷元間差分后鐘差的數(shù)值積分，即

(11)

式中，n表示在初始?xì)v元到歷元k的積分變量；dn為積分步長(zhǎng)，通常等于觀測(cè)的歷元間隔。本文中初始的歷元鐘差設(shè)為0 s。

1.4 離散小波變換去噪

GNSS載波相位觀測(cè)值經(jīng)過上文修正后，其殘差可表示為

(12)

(1) 分解。選用合適的小波基對(duì)數(shù)據(jù)進(jìn)行多級(jí)分解，獲得小波系數(shù)。去噪的小波基需要正交或雙正交，本文選用較為常用的多貝西極限正交相位小波，消失矩設(shè)為2。分解層數(shù)(level)的確定方法如下

level=max(5,[log2n])

(13)

式中，[·]表示向下取整運(yùn)算；n為數(shù)據(jù)長(zhǎng)度。

(3) 重構(gòu)。利用閾值處理后的多級(jí)小波系數(shù)進(jìn)行小波重構(gòu)，獲取目標(biāo)信號(hào)。

1.5 電離層閃爍信號(hào)提取及閃爍因子構(gòu)建

GNSS載波相位觀測(cè)值經(jīng)過1.4所提出方法進(jìn)一步修正后，觀測(cè)噪聲已經(jīng)被大大削弱了，其殘差可以表示為

(14)

(15)

2 數(shù)據(jù)簡(jiǎn)介

為確定構(gòu)建所提出閃爍因子的相關(guān)經(jīng)驗(yàn)參數(shù)及驗(yàn)證其對(duì)電離層閃爍的探測(cè)效果，本文選取了加拿大高緯度北極電離層閃爍監(jiān)測(cè)網(wǎng)絡(luò)(CHAIN)中的11測(cè)站[23]，它們的分布如圖1所示，均位于極隙區(qū)至極光區(qū)內(nèi)，其中chuc測(cè)站具有并址磁測(cè)站，可以提供時(shí)間分辨率為1 min的地磁場(chǎng)強(qiáng)度信息。每個(gè)測(cè)站均配備Septentrio PolaRxS Pro電離層閃爍監(jiān)測(cè)接收機(jī)和Septentrio的PolaN GG天線。CHAIN中的接收機(jī)被設(shè)定為接收雙頻GPS信號(hào)，且可以利用采樣頻率為50 Hz的觀測(cè)數(shù)據(jù)直接產(chǎn)生電離層相位閃爍因子σφ，本文將其視為評(píng)價(jià)1 Hz數(shù)據(jù)閃爍因子的參考值。CHAIN同時(shí)可以提供1 Hz的GPS觀測(cè)數(shù)據(jù)，本文利用其來(lái)研究所提出電離層相位閃爍因子的有效性。

本文選擇以上11個(gè)測(cè)站于2009年DOY 230至DOY 260、2013年DOY 120至DOY 150、2015年DOY 90至DOY 110、2017年DOY 100至DOY 130、2018年DOY 220至DOY 250和2020年DOY 80至DOY 112，共6組188 d的觀測(cè)數(shù)據(jù)作為研究對(duì)象，每組數(shù)據(jù)中均包括強(qiáng)弱不同的地磁活動(dòng)，利用每組數(shù)據(jù)獲得的結(jié)論相似，為避免重復(fù)，這里重點(diǎn)介紹最后一組的地磁強(qiáng)弱情況。電離層閃爍與地磁活動(dòng)相關(guān)[24]，因此本文選擇Kp、Dst和ASY 3個(gè)參數(shù)來(lái)反映的地磁活動(dòng)情況(如圖2所示)。Kp指數(shù)表征全球地磁風(fēng)暴的幅度，其中Kp≥5表示發(fā)生地磁風(fēng)暴。Dst指數(shù)表示地球表面南北極與極赤道的軸向?qū)ΨQ擾動(dòng)磁場(chǎng)。Dst中主要擾動(dòng)為負(fù)值，即地磁場(chǎng)的減少，主要產(chǎn)生于磁層中的電子環(huán)流；而Dst的正向擾動(dòng)主要是太陽(yáng)風(fēng)作用在磁層造成的。ASY表征中緯度地區(qū)地磁擾動(dòng)，其有兩個(gè)分量，即水平分量(偶極子-極方向，H)和垂直分量(東-西方向，D)，其中ASY-H的變化趨勢(shì)被證實(shí)與AE指數(shù)具有強(qiáng)相關(guān)性。AE指數(shù)通常被用來(lái)表征極光區(qū)電子活躍程度，但由于近期的AE指數(shù)無(wú)法獲得，故本文用ASY指數(shù)近似替代AE指數(shù)來(lái)表征研究區(qū)域的電子活躍程度。由圖2可以看出，2020年DOY 111具有較強(qiáng)的地磁風(fēng)暴，該日的觀測(cè)數(shù)據(jù)將是本文重點(diǎn)研究對(duì)象。

3 利用小波變換提取電離層閃爍的經(jīng)驗(yàn)參數(shù)確定

如上文所述，本文利用連續(xù)小波變換的方法將電離層閃爍信號(hào)與電離層折射信號(hào)相分離，而恰當(dāng)?shù)倪x取小波變換的小波基和小波變換相關(guān)參數(shù)會(huì)提高該分離效果，進(jìn)而提高所構(gòu)建相位閃爍因子的準(zhǔn)確性。下面將對(duì)電離層閃爍時(shí)頻分析的小波基及對(duì)稱參數(shù) (symmetry parameter) 和時(shí)間帶寬積(time-bandwidth product)的選取進(jìn)行詳細(xì)闡述。

圖1 本文選用的11個(gè)測(cè)站的分布Fig.1 Distribution of the 11 stations selected in this paper

圖2 2020年年積日80至112日的地磁活躍程度Fig.2 Degree of geomagnetic activity on DOY 80 to DOY 112 of 2020

3.1 小波基的選取

可進(jìn)行時(shí)頻分析的連續(xù)小波變換通常有3種小波基，即廣義Morse小波(generalized morse wavelets)、解析Morlet小波(analytic morlet wavelet)和Bump小波。本文采用受電離層閃爍影響較為嚴(yán)重的觀測(cè)弧段(arcc測(cè)站PRN07衛(wèi)星2020年DOY 111 15:00至19:00的觀測(cè)數(shù)據(jù))，對(duì)上述3種小波的時(shí)頻分析能力進(jìn)行比較，方法為：首先以1 Hz為采樣頻率、以該弧段的觀測(cè)歷元數(shù)14 400為信號(hào)長(zhǎng)度(其中Morse小波的對(duì)稱參數(shù)和時(shí)間帶寬分別設(shè)置為3和60)分別建立以上3種小波的小波濾波器庫(kù)，并利用Matlab提供的頻響特性函數(shù)(freqz)繪制幅頻響應(yīng)圖(圖3)；然后利用建立的濾波器庫(kù)對(duì)該弧段觀測(cè)數(shù)據(jù)進(jìn)行小波變換，繪制時(shí)頻譜圖(圖4(a)—圖4(c)，觀察其對(duì)電離層閃爍信息的響應(yīng)能力；接著以0.1 Hz至0.4 Hz為特征頻率區(qū)間，將對(duì)應(yīng)的小波系數(shù)進(jìn)行逆變換，確定每種小波基對(duì)應(yīng)的電離層相位閃爍因子(圖4(d)—圖4(f)，依據(jù)其與參考值(σφ)的契合程度，尤其是在發(fā)生電離層閃爍時(shí)段(σφ≥0.2 rad)的契合程度，選取恰當(dāng)?shù)男〔ɑ?/p>

從幅頻響應(yīng)圖(圖3)可以看出，Bump小波的頻率響應(yīng)區(qū)間最小(如圖3(a)所示)，而電離層閃爍通?？缭降念l率區(qū)間比較大，使得利用Bump提取的電離層閃爍信號(hào)強(qiáng)度會(huì)明顯低于電離層閃爍接收機(jī)高頻觀測(cè)數(shù)據(jù)給出的相位閃爍因子(如圖4(a)所示)，導(dǎo)致最終確定的閃爍因子遠(yuǎn)低于參考值(如圖4(d)所示)。如圖3(b)和圖3(c)所示，Morlet小波和Morse小波在全頻率域均有響應(yīng)，但Morse小波在對(duì)高頻信號(hào)響應(yīng)更佳。Morlet小波和Morse小波均能較好地反映出電離層閃爍信號(hào)的時(shí)頻信息(如圖4(b)和圖4(c)所示)，但利用Morlet小波變換后所得到的電離層相位閃爍因子的幅值明顯大于參考值(如圖4(e)所示)，進(jìn)而導(dǎo)致過多的誤報(bào)，而利用Morse小波獲得相位閃爍因子與參考值基本相當(dāng)(如圖4(f)所示)，且Morse小波可以通過調(diào)整其對(duì)稱參數(shù)和時(shí)間帶寬積參數(shù)而進(jìn)一步提高所構(gòu)建閃爍因子的準(zhǔn)確性，所以本文選擇Morse小波作為連續(xù)小波變換的小波基。

圖3 3種小波基的幅頻響應(yīng)圖Fig.3 Amplitude-frequency response diagrams of the three wavelets

圖4 arcc測(cè)站PRN07衛(wèi)星2020年DOY 111 15時(shí)至19時(shí)觀測(cè)弧段的時(shí)頻譜圖和相應(yīng)閃爍因子Fig.4 The time-frequency spectrum and corresponding scintillation index based on the observation arc from 15h to 19h of PRN07 collected on DOY 111 of 2020 at arcc station

3.2 對(duì)稱參數(shù)和時(shí)間帶寬積的選取

對(duì)稱參數(shù)和時(shí)間帶寬積是利用Morse小波基時(shí)頻分析的重要參數(shù)。Morse小波的對(duì)稱參數(shù)通過解調(diào)偏斜(demodulate skewness)控制著小波在時(shí)間上的對(duì)稱性；時(shí)間帶寬積與小波的持續(xù)時(shí)間成正比，決定時(shí)域小波的中心窗口在峰值頻率處可以容納振蕩的個(gè)數(shù)。二者共同決定著小波的形狀，進(jìn)而影響小波變換的效果，因此恰當(dāng)?shù)倪x擇小波參數(shù)對(duì)于電離層閃爍信息的準(zhǔn)確提取至關(guān)重要。

為確定最優(yōu)Morse小波參數(shù)，進(jìn)而準(zhǔn)確提取1 Hz GNSS觀測(cè)值中的電離層閃爍信息，首先分別利用有電離層閃爍(arcc測(cè)站2020年DOY 111 PRN07衛(wèi)星16:00至18:00)和無(wú)電離層閃爍(arcc測(cè)站2020年DOY 111 PRN02衛(wèi)星00:00至03:30)弧段以及兩組距離較遠(yuǎn)的參數(shù)組合作特例分析，以確定最優(yōu)參數(shù)組合的所在區(qū)間。分析方法如下：利用對(duì)稱參數(shù)3和6、時(shí)間帶寬積60和20組成兩個(gè)參數(shù)組合(3,60)和(6,20)，分別對(duì)以上兩種弧段的觀測(cè)數(shù)據(jù)進(jìn)行連續(xù)小波變換，并對(duì)特征頻率區(qū)間(0.1～0.4 Hz)的小波參數(shù)進(jìn)行連續(xù)小波逆變換，構(gòu)建相應(yīng)的相位閃爍因子；以參考值為依據(jù)對(duì)不同參數(shù)組合構(gòu)建的閃爍因子進(jìn)行3方面評(píng)估：相關(guān)性、線性擬合的斜率和殘差的均方根誤差(RMS)。相關(guān)性越高，說(shuō)明二者的相似程度越高；線性擬合的斜率越接近于1，殘差的均方根誤差越小說(shuō)明構(gòu)建的閃爍因子越接近參考值。

特例分析的結(jié)果如圖5所示。由圖5(a)—圖5(d)可以看出利用相同參數(shù)組合作用于閃爍強(qiáng)度不同的數(shù)據(jù)獲得的閃爍指數(shù)與參考值的符合程度不同，利用不同參數(shù)作用于同一弧段數(shù)據(jù)所得閃爍指數(shù)也不相同。為評(píng)價(jià)這兩種參數(shù)組合參數(shù)表現(xiàn)優(yōu)劣，進(jìn)一步對(duì)由小波變換構(gòu)建的閃爍因子與參考值做線性擬合并進(jìn)行相關(guān)性分析和殘差的均方根誤差分析，其結(jié)果如圖5(e)—圖5(l)所示?？梢钥闯?，在有電離層閃爍發(fā)生時(shí)段，本文所提出方法構(gòu)建的閃爍因子與參考值的相關(guān)性更高，其擬合斜率更接近于1，但由于在強(qiáng)閃爍發(fā)生時(shí)，閃爍幅值較大，導(dǎo)致殘差的均方根誤差較大；而在閃爍強(qiáng)度小于0.05 rad時(shí)，由于受到未能剔除的觀測(cè)噪聲影響，估計(jì)的閃爍指數(shù)呈現(xiàn)出較強(qiáng)的隨機(jī)性，進(jìn)而降低了其與參考值的相關(guān)性。利用參數(shù)組合(3,60)獲得的閃爍指數(shù)的擬合斜率小于1而參數(shù)組合(6,20)獲得擬合斜率大于1，因此通過以上特例分析可以看出，小波參數(shù)的選擇會(huì)影響所估計(jì)的閃爍指數(shù)，最佳的參數(shù)組合應(yīng)該位于對(duì)稱參數(shù)3～6、時(shí)間帶寬積20～60。

在大體確定最佳參數(shù)組合潛在范圍的基礎(chǔ)上，本文進(jìn)一步利用全部11測(cè)站188 d的觀測(cè)數(shù)據(jù)選取最佳的參數(shù)組合，方法為：分別以1和5為步長(zhǎng)選取位于3至6區(qū)間和20至60區(qū)間的對(duì)稱參數(shù)和時(shí)間帶寬積，組成多種參數(shù)組合；基于以上數(shù)據(jù)，利用多種組合參數(shù)的Morse小波構(gòu)建相應(yīng)的相位閃爍因子，并比較其與參考值線性擬合的斜率、相關(guān)性和殘差的均方根誤差。圖6和圖7分別統(tǒng)計(jì)了由188 d的數(shù)據(jù)獲得的上述3種指標(biāo)在有閃爍發(fā)生時(shí)段和無(wú)閃爍發(fā)生時(shí)段的變化情況，并以誤差線的形式展示每種指標(biāo)在對(duì)應(yīng)參數(shù)組合下波動(dòng)的標(biāo)準(zhǔn)差?？梢钥闯觯陨戏治隹梢源_定獲得較好結(jié)果的參數(shù)組合所在的大致范圍。在有閃爍發(fā)生時(shí)段，對(duì)稱系數(shù)為3時(shí)可以給出更為優(yōu)良的結(jié)果，其與時(shí)間帶寬積45、50、55和60組成的參數(shù)組合具有相似且較為優(yōu)良的表現(xiàn)，均可以使線性擬合斜率趨近于1，相關(guān)性達(dá)到92%而殘差的均方根誤差降至0.03 rad；而在無(wú)閃爍發(fā)生時(shí)段，依舊是對(duì)稱系數(shù)3表現(xiàn)更佳，其與時(shí)間帶寬積30、35、40、45組成的參數(shù)組合可以獲得較好的結(jié)果，但其相關(guān)性要略低于有閃爍發(fā)生時(shí)段，筆者認(rèn)為這主要是由于觀測(cè)噪聲在無(wú)閃爍發(fā)生時(shí)段帶來(lái)的影響導(dǎo)致估計(jì)的閃爍因子呈現(xiàn)一定的隨機(jī)性導(dǎo)致。

為了從圖6和圖7給出的較優(yōu)參數(shù)組合范圍中確定最佳參數(shù)組合，本文進(jìn)一步分析了每個(gè)測(cè)站的最優(yōu)和次優(yōu)參數(shù)組合，并計(jì)算了最優(yōu)和次優(yōu)參數(shù)組合在線性擬合斜率、相關(guān)性和殘差的均方根誤差的指標(biāo)差，計(jì)算方法如下

線性擬合斜率指標(biāo)差=

(16)

相關(guān)性指標(biāo)差=最優(yōu)參數(shù)組合相關(guān)性-

次優(yōu)參數(shù)組合相關(guān)性

(17)

圖5 利用arcc測(cè)站2020年DOY 111的數(shù)據(jù)對(duì)組合參數(shù)(3,60)和(6,20)特例分析的結(jié)果Fig.5 Special case analysis of two parameter groups (3,60) and (6,20) observation arcs collected on DOY 111 of 2020 at arcc station

4 與ROTI指數(shù)的比較

為驗(yàn)證所提出的閃爍因子σφf(shuō),wavelet的有效性，本文利用全部11測(cè)站6組共188 d的觀測(cè)數(shù)據(jù)，以CHAIN網(wǎng)絡(luò)高頻采樣數(shù)據(jù)輸出的相位閃爍因子σφ為參考值，對(duì)所提出的閃爍因子和ROTI指數(shù)進(jìn)行對(duì)比。由于每組數(shù)據(jù)獲得的結(jié)果類似，為避免重復(fù)，這里重點(diǎn)介紹由第6組(即2020年DOY 80至DOY 112)觀測(cè)數(shù)據(jù)獲得的結(jié)果。ROTI是當(dāng)前最為常用的基于低采樣頻率觀測(cè)數(shù)據(jù)的閃爍因子[25]，被諸多學(xué)者用來(lái)監(jiān)測(cè)電離層閃爍[26-32]。ROTI表示斜路徑總電子含量(STEC)的變化率，其計(jì)算方法可參考文獻(xiàn)[12]。ROTI中的STEC直接由雙頻載波無(wú)幾何關(guān)系組合獲得，因此ROTI容易受到周跳的影響。當(dāng)前已有方法難以探測(cè)并修復(fù)全部周跳，尤其對(duì)于低衛(wèi)星高度角或有閃爍發(fā)生的弧段，為此本文將截止高度角設(shè)為20°并采用文中所述方法，以最大限度減弱周跳對(duì)ROTI的影響。值得一提的是，本文所提出的閃爍因子構(gòu)建方法中由于引入了小波變換技術(shù)，可以較為容易地克服未正確探測(cè)及修復(fù)的周跳。

圖6 有閃爍發(fā)生時(shí)段，所構(gòu)建閃爍因子與參考值之間的統(tǒng)計(jì)分析Fig.6 Statistical analysis of the reference and the proposed scintillation index

圖7 無(wú)閃爍發(fā)生時(shí)段，所構(gòu)建閃爍因子與參考值之間的統(tǒng)計(jì)分析Fig.7 Statistical analysis of the reference and the proposed scintillation index

表1 有閃爍發(fā)生時(shí)段，各測(cè)站最優(yōu)與次優(yōu)參數(shù)組合及其指標(biāo)差

表2 無(wú)閃爍發(fā)生時(shí)段，各測(cè)站最優(yōu)與次優(yōu)參數(shù)組合及其指標(biāo)差

4.1 電離層閃爍的監(jiān)測(cè)效果

為驗(yàn)證所提出閃爍因子對(duì)電離層閃爍的監(jiān)測(cè)效果，本文利用chuc測(cè)站2020年DOY 110和DOY 111的觀測(cè)數(shù)據(jù)，對(duì)比分析了σφf(shuō),wavelet、σφ和ROTI在每小時(shí)時(shí)間內(nèi)所探測(cè)出的電離層閃爍歷元的數(shù)量。選用該測(cè)站的原因是因?yàn)閏huc測(cè)站有并址磁測(cè)站，而電離層閃爍與本地的地磁場(chǎng)活動(dòng)直接相關(guān)，因而可以利用磁測(cè)站提供的地磁場(chǎng)數(shù)據(jù)，驗(yàn)證所探測(cè)出受電離層閃爍影響歷元數(shù)量的正確性。需要指出的是σφf(shuō),wavelet與ROTI的時(shí)間分辨率與GPS數(shù)據(jù)的采樣頻率相同(本文中為1 s)，而σφ的時(shí)間分辨率為60 s，為方便對(duì)比，本文將σφf(shuō),wavelet與ROTI按照σφ中對(duì)應(yīng)的時(shí)刻進(jìn)行了重采樣。本文采用經(jīng)驗(yàn)值的方式確定電離層閃爍的閾值，即σφf(shuō),wavelet與σφ為0.2 rad、ROTI為0.1 TECU/min[33]。當(dāng)閃爍因子的幅值大于閾值時(shí)即判定發(fā)生電離層閃爍。

圖8展示了chuc測(cè)站在第110和111天的地磁場(chǎng)變化情況及每小時(shí)受電離層閃爍影響的歷元數(shù)量。由圖8(a)可以看出DOY 111的地磁活動(dòng)要強(qiáng)于DOY 110。由圖8(b)和8(c)可以看出3種閃爍因子所探測(cè)出發(fā)生電離層閃爍的歷元均分布在地磁場(chǎng)較為活躍的時(shí)段，這在一定程度上再次證明了電離層閃爍與本地的地磁場(chǎng)活動(dòng)直接相關(guān)的結(jié)論以及3種閃爍因子對(duì)電離層閃爍監(jiān)測(cè)的有效性；σφf(shuō),wavelet與ROTI所探測(cè)出的發(fā)生電離層閃爍的歷元數(shù)量與參考值所給出的數(shù)量基本相同(尤其是L2載波)，說(shuō)明σφf(shuō),wavelet與ROTI均可以較為準(zhǔn)確地監(jiān)測(cè)電離層閃爍；盡管ROTI無(wú)法給出每個(gè)頻點(diǎn)上的閃爍數(shù)量，但是其探測(cè)結(jié)果與σφf(shuō),wavelet在L1、L2頻點(diǎn)探測(cè)結(jié)果均相同，這主要是由于在高緯度地區(qū)，電離層閃爍對(duì)GPS的L1、L2頻點(diǎn)的影響幾乎相同[34]。綜上所述，本文所提出的閃爍因子σφf(shuō),wavelet和ROTI因子均能準(zhǔn)確地探測(cè)出電離層閃爍的發(fā)生，可以用于電離層閃爍的定性監(jiān)測(cè)研究。

4.2 電離層閃爍的幅值估計(jì)準(zhǔn)確度

由4.1中的試驗(yàn)可以證明σφf(shuō),wavelet和ROTI均可以較為準(zhǔn)確地探測(cè)出電離層閃爍發(fā)生，但無(wú)法證明二者估計(jì)的電離層閃爍幅值的準(zhǔn)確度。本文將相關(guān)性作為為衡量電離層閃爍因子估計(jì)的幅值準(zhǔn)確度的指標(biāo)。在σφf(shuō),wavelet和ROTI均能準(zhǔn)確探測(cè)出電離層閃爍發(fā)生的前提下，二者與σφ的相關(guān)性越高，說(shuō)明其對(duì)電離層閃爍估計(jì)的幅值越接近σφ。本文將全部11個(gè)測(cè)站188 d的觀測(cè)數(shù)據(jù)分成有閃爍發(fā)生(σφ≥0.2 rad)和無(wú)閃爍發(fā)生(σφ<0.2 rad)兩種類別，分別研究由每顆衛(wèi)星獲得的σφf(shuō),wavelet、ROTI與σφ的相關(guān)性。

圖9和10分別展示了利用2020年DOY 80至DOY 112數(shù)據(jù)獲得的有閃爍和無(wú)閃爍發(fā)生時(shí)段的相關(guān)性結(jié)果，其中，空白區(qū)域表示該衛(wèi)星不可用或無(wú)閃爍發(fā)生。因在相同緯度附近測(cè)站的結(jié)果類似和篇幅限制，本文僅展示了arcc、sarc、repc和chuc這4個(gè)測(cè)站的結(jié)果。由圖9可以看出，σφf(shuō),wavelet與σφ的相關(guān)性在全部有閃爍發(fā)生的弧段均在0.8以上，且高于ROTI與σφ的相關(guān)性，說(shuō)明在有閃爍發(fā)生時(shí)，σφf(shuō),wavelet的幅值更加接近σφ。由圖10可以看出，在無(wú)閃爍發(fā)生的時(shí)段，緯度較高的測(cè)站(arcc和sarc)所給出的σφf(shuō),wavelet、ROTI與σφ的相關(guān)性均要高于緯度較低的測(cè)站，這主要是因?yàn)閍rcc和sarc測(cè)站位于極隙區(qū)，而其余兩個(gè)測(cè)站分布在極光區(qū)，極隙區(qū)的磁場(chǎng)線呈現(xiàn)開放狀態(tài)，導(dǎo)致該區(qū)域的電離層更容易受到太陽(yáng)風(fēng)的影響而變得活躍，進(jìn)而導(dǎo)致電離層閃爍的發(fā)生頻率更高。盡管此時(shí)極隙區(qū)和極光區(qū)閃爍均沒有超過閾值，但是極隙區(qū)的閃爍幅度要高于極光區(qū)的閃爍，使得由1 Hz采樣頻率觀測(cè)數(shù)據(jù)估算出來(lái)的閃爍因子受觀測(cè)噪聲的影響更小，進(jìn)而提高了σφf(shuō),wavelet、ROTI與σφ的相關(guān)性。圖10也可以看出，對(duì)于極隙區(qū)測(cè)站(arcc和sarc)，σφf(shuō),wavelet的相關(guān)性明顯優(yōu)于ROTI，但該圖難以反映二者的相關(guān)性在其他區(qū)域測(cè)站(如repc、chuc)的情況，需進(jìn)一步統(tǒng)計(jì)分析。

圖8 σφf(shuō),wavelet、ROTI及σφ對(duì)chuc測(cè)站電離層閃爍的探測(cè)結(jié)果對(duì)比Fig.8 Ionospheric scintillation detection results provided by σφf(shuō),wavelet、 ROTI and σφ utilizing the observations from chuc station

為進(jìn)一步比較σφf(shuō),wavelet與ROTI在無(wú)閃爍發(fā)生時(shí)段的電離層閃爍的幅值估計(jì)準(zhǔn)確度，本文對(duì)全部的σφf(shuō),wavelet、ROTI與σφ的相關(guān)性按照下式進(jìn)行如下統(tǒng)計(jì)分析

(18)

式中，corr表示相關(guān)性；s表示測(cè)站衛(wèi)星；NS表示測(cè)站在一天內(nèi)全部可視衛(wèi)星數(shù)量；DOY表示年積日；P表示σφf(shuō),wavelet與σφ的相關(guān)性大于ROTI與σφ的相關(guān)性的概率。式(18)表示的分析方法為：將每個(gè)測(cè)站每天數(shù)據(jù)獲得的σφf(shuō),wavelet與σφ的相關(guān)性減去ROTI與σφ的相關(guān)性，找出差值大于0的衛(wèi)星的數(shù)量，并除以全部的衛(wèi)星數(shù)量，若最終的結(jié)果大于50%，則說(shuō)明在該測(cè)站該天，σφf(shuō),wavelet比ROTI有更大的可能性獲得相關(guān)性更高的閃爍因子。全部11個(gè)測(cè)站于2020年DOY 80至DOY 112觀測(cè)數(shù)據(jù)的統(tǒng)計(jì)結(jié)果如圖11所示，相比于ROTI，σφf(shuō),wavelet有72.9%的概率獲得與σφ相關(guān)性更高的閃爍因子，即給出的閃爍幅值更為準(zhǔn)確。

圖9 有閃爍發(fā)生時(shí)段σφf(shuō),wavelet、ROTI與σφ(參考值)的相關(guān)性比較Fig.9 Comparison of the correlation provided by σφf(shuō),wavelet、 ROTI with σφ (the reference)

5 結(jié) 論

本文以1 Hz采樣頻率GNSS載波相位觀測(cè)值為基礎(chǔ)，通過大地測(cè)量趨勢(shì)分離、精密單點(diǎn)定位、無(wú)電離層組合估計(jì)接收機(jī)鐘差、周跳探測(cè)與修復(fù)、離散小波去噪等技術(shù)，剔除了載波相位觀測(cè)值中的衛(wèi)地幾何距離以及固體潮、海潮、大氣潮、極潮、天線相位中心、相對(duì)論、相位纏繞、對(duì)流層延遲、接收機(jī)鐘差、周跳和觀測(cè)噪聲誤差，并詳細(xì)研究了利用連續(xù)小波變換技術(shù)提取電離層閃爍信號(hào)中的最優(yōu)小波基及相關(guān)參數(shù)的經(jīng)驗(yàn)值，進(jìn)而構(gòu)建了基于1 Hz觀測(cè)數(shù)據(jù)的電離層相位閃爍因子。為了驗(yàn)證所提出閃爍因子的有效性和準(zhǔn)確性，本文以高采樣頻率數(shù)據(jù)的相位閃爍因子為參考值，采用加拿大高緯度北極電離層閃爍監(jiān)測(cè)網(wǎng)絡(luò)(CHAIN)中11個(gè)測(cè)站188 d的1 Hz觀測(cè)值對(duì)所提出閃爍因子的有效性進(jìn)行驗(yàn)證，并與已經(jīng)廣泛使用的ROTI因子進(jìn)行對(duì)比，結(jié)果表明在電離層閃爍發(fā)生時(shí)段，本文所提出的閃爍因子與ROTI均能有效地探測(cè)到電離層閃爍的發(fā)生，但本文提出的閃爍因子所給出的電離層閃爍的幅值更接近參考值，準(zhǔn)確度更高，尤其在有強(qiáng)電離層閃爍發(fā)生的時(shí)段。

文獻(xiàn)[35]指出接收機(jī)相位偏差存在受溫度影響的短期變化，而本文為了簡(jiǎn)便，認(rèn)為相位偏差在同一觀測(cè)弧段不變化，未來(lái)的研究若能考慮到短期的溫度變化，或可在一定程度上提高本文所構(gòu)建相位閃爍因子的對(duì)電離層閃爍的監(jiān)測(cè)精度。需要指出的是，構(gòu)建該閃爍因子的原理也可用于測(cè)地型接收機(jī)的觀測(cè)值，由于篇幅限制，筆者將在未來(lái)通過試驗(yàn)驗(yàn)證其可靠性，以拓展該閃爍因子的應(yīng)用場(chǎng)景。

圖10 無(wú)閃爍發(fā)生時(shí)段，σφf(shuō),wavelet、ROTI與σφ(參考值)的相關(guān)性比較Fig.10 Comparison of the correlation provided by σφf(shuō),wavelet、ROTI with σφ (the reference)

圖11 σφf(shuō),wavelet與σφ的相關(guān)性大于ROTI與σφ的相關(guān)性概率統(tǒng)計(jì)Fig.11 Statistics of the probability that the correlation provided by σφf(shuō),wavelet and σφ is greater than that provided by ROTI and σφ