孫佳欣 彭懷云 李清亮

（中國電波傳播研究所 電波環(huán)境特性及?；夹g(shù)重點實驗室，青島 266107）

引 言

低電離層(60～120 km高度)是甚低頻(very low frequency, VLF)電磁波地-電離層波導(dǎo)傳播模型的上邊界[1]，其電子濃度、電子溫度、碰撞頻率等參量的分布、擾動變化是影響VLF遠距離通信效能、導(dǎo)航定位精度的重要因素[2]. 因此，低電離層參量的分布和擾動特性研究是VLF電波傳播研究的重要方向之一.

1990年，Inan[3]分析了南極Palmer站接收的NAA臺VLF傳播數(shù)據(jù)，發(fā)現(xiàn)該VLF信號幅度被NAU臺發(fā)射的VLF信號調(diào)制，實驗證實了VLF波能顯著改變低電離層電子溫度和濃度的理論預(yù)言.Inan等人[4]、Rodriguez等人[5]分析了NSS、NAU、NAA三個VLF臺發(fā)射的強VLF電磁波對夜間低電離層電子溫度和電子密度的影響，發(fā)現(xiàn)在70～90 km高度電子溫度升高2.5倍，電子密度下降26%. 1997年，劉立波等人[6]在Rodriguez的基礎(chǔ)上，分析了電離層中氧氣、氮氣密度等背景參量對電子溫度、濃度擾動的影響，發(fā)現(xiàn)中性大氣數(shù)密度越低或溫度越高，電子參數(shù)變化越大. 2007年，Peter等人[7]分析了閃電激發(fā)的VLF電波對電子密度的擾動，發(fā)現(xiàn)85 km處D層電子密度提高了約15%. 2013年，Graf等人[8]計算了VLF波引起的低電離層擾動范圍和跨低電離層衰減的估計值. 2017年，馮陽[9]推導(dǎo)了全波分析法下的VLF波傳播模型，并對不同輻射源參數(shù)、地磁場參數(shù)、電離層參數(shù)下的VLF波傳輸特性進行了分析.但他研究的重點是電離層中的VLF波傳播情況，并沒有對VLF波引起的電離層擾動情況進行分析.2018年，ATICI等人[10]利用Glukhov-Pasko-Inan電離層化學(xué)經(jīng)驗?zāi)Ｐ陀嬎懔死纂娂ぐl(fā)的VLF電磁波對低電離層電子密度的影響，發(fā)現(xiàn)由于負離子與氧、氮原子的反應(yīng)，在場強為10 V/m的情況下，低電離層電子濃度提高了5.5倍，與無線電波對電離層的擾動相比，雷電引起的電子濃度擾動幅度更大. 同年，成琦等人[11]在電離層自恰加熱模型的基礎(chǔ)上，引入電子復(fù)合和吸附效應(yīng)，對高頻電波引起的電離層擾動進行仿真. 成琦使用的電離層化學(xué)反應(yīng)模型和擾動模型與本文相似，不同之處在于，本文使用的是VLF波擾動電離層. 由于低電離層電子濃度、碰撞頻率參數(shù)的高度分布在一個波長范圍內(nèi)劇烈變化，還需考慮各層反射、透射波的相互疊加情況.

本文首先通過全波解反射系數(shù)和波模振幅遞歸算法，引入不同入射角的平面波，分析了各層中上行和下行波的能量密度分布規(guī)律. 然后基于Glukhov-Pasko-Inan低電離層化學(xué)反應(yīng)模型，利用梯度下降法[12]，分析了不同入射角平面波對電子溫度、濃度改變的影響規(guī)律. 本研究可為細致分析偶極子源(例如VLF臺、雷電)輻射的不同入射角VLF波對低電離層擾動的貢獻以及擾動現(xiàn)象的分析、認識提供思路.

1 計算方法

1.1 平面波入射的全波法

建立電離層水平分層均勻模型，x為正北方向，y為正東方向，z為垂直向上方向. 從地面至電離層上邊界(700 km)，按高度分為N層，每層的邊界為zi(i=0,1,2,···,N)，z0代表地面，zN為電離層上邊界.由于VLF電波頻率遠高于離子回旋頻率，離子影響可忽略. 電離層相對磁導(dǎo)率為1，時諧因子為 e?jωt.θ為入射角，為波矢量與x方向夾角.l、m、n為地磁場與x、y、z軸夾角的余弦，l=cosφBcosθB,m=sinφBcosθB,n=sinθB, θB為 磁傾角，φB為地磁場水平分量與正北方向夾角.

電離層中每層的電極化率矩陣由該層電子濃度、碰撞頻率決定，每層中電極化率矩陣[13]表示為

全波法將各分層中水平電場、磁場進行二維空間域傅里葉變換，并代入Maxwell方程，消去和，可得無源層中水平電、磁場的矩陣形式：

式中：

基于各層分界面水平電、磁場連續(xù)的邊界條件，利用反射系數(shù)遞歸算法，可獲得各分層界面對上行、下行波的反射系數(shù)矩陣當(dāng)?shù)?層(為空氣層)中入射波為TM平面波時入射磁場為水平方向，且垂直于入射面，其二維空間域傅里葉變換表示為

式中：a為 磁場幅度，單位為A/m. 入射波的波模振幅表示為

同理，當(dāng)入射波為TE波時入射電場為水平方向，且垂直于入射面，其二維空間域傅里葉變換表示為

式中，β為電場幅度，單位為V/m. 對應(yīng)的上行波模振幅為

在上述平面波引入的基礎(chǔ)上，利用波模振幅遞歸算法[14]，獲取各層中的上行、下行波模振幅. 再依據(jù)式(4)、Maxwell方程和傅里葉逆變換計算各層中的上下行4個波模的電場和磁場. 進而基于坡印廷定理計算出各波模的平均能流密度

1.2 低電離層參數(shù)擾動的梯度下降法

在低電離層中VLF電波能量主要被電子吸收，其他粒子的吸收可忽略[15]，電離層各分層中單位時間、單位體積內(nèi)吸收的能量

是4種波模吸收能量之和，單位為W/(m3·s?1)；(m=1, 2, 3, 4)為4種波模的吸收指數(shù)，

式中：c為光速；Im(·)表示取虛部.

各分層中電子溫度變化由電子能量方程確定：

式中：k為玻爾茲曼常數(shù)；Li為各層電子能量損失項.

能量損失主要來源于電子與其他粒子的碰撞、中性粒子的能級激發(fā)等，其計算模型為Stubbe[16]能量損失模型.

通過式(12)可以得到電子溫度隨時間變化速率dTe/dt，選取時間步長Δt，用差分法將電子溫度的變化過程離散化，相鄰兩時刻電子溫度之間關(guān)系可表示為

式中：Te0為前一時刻的電子溫度；Te為當(dāng)前時刻電子溫度. 電子溫度的變化改變了電離層中電子與其他帶電粒子的化學(xué)反應(yīng)過程，使其濃度發(fā)生變化. 各粒子的濃度變化速率可通過Glukhov-Pasko -Inan低電離層化學(xué)反應(yīng)模型[17]表達：

式中：Ne、N?、N+、Nx+分別表示電子、負離子、正離子、團簇離子的濃度；I0是電離層自然狀態(tài)下的電子產(chǎn)生率，白天太陽是主要的電離源，在夜間主要由月亮反射的太陽光、人工源、其他自然因素或宇宙背景下的電離源產(chǎn)生；γ是電子從負離子上的脫離率；β是電子附著到中性粒子的效率；αd是電子與正離子復(fù)合系數(shù)；是電子與團簇粒子的結(jié)合率；αi是負離子與正離子和水合正離子之間的復(fù)合系數(shù)；B是初級正離子轉(zhuǎn)化為水和離子的轉(zhuǎn)化率. 電子溫度的擾動改變了帶電粒子間的反應(yīng)系數(shù)，進而對各粒子成分濃度造成影響. 與式(13)同理，相鄰兩時刻的各粒子濃度關(guān)系可表示為

將初始電子溫度、濃度等參量帶入式(1)～(11)，用全波法計算電波能流分布、低電離層能量吸收，進而得到時間步長Δt后的電子溫度、電子和離子濃度. 如此循環(huán)迭代，直到電子溫度和帶電粒子濃度不再發(fā)生變化，即可實現(xiàn)VLF電磁波對低電離層擾動的計算仿真.

為加快迭代模型的收斂速度，Δt的選擇采取梯度下降法，將電子濃度的變化率作為損失函數(shù)，使Δt隨電子濃度變化率改變：

式中，ζ為學(xué)習(xí)率，取ζ=1×10?9.

2 計算分析與討論

2.1 擾動結(jié)果對比

文中使用的電離層中性大氣參數(shù)來源于MSISE-90模型. 正離子密度、電子溫度、80 km高度以上的電子數(shù)密度來源于IRI-2016或文獻[4]中使用過的經(jīng)驗?zāi)Ｐ? 在無特殊說明的情況下，平面波入射方向為垂直向上，即θ=0°、φ=0°，地磁場方向為θB=60°、φB= 0°，大小為50μT. 入射波頻率設(shè)置為f= 10 kHz，輸入能流為8×10?6J/(m3·s?1).

在上述條件下，分別用Rodriguez[18]方法和本文方法對電離層擾動情況進行仿真. 圖1(a)、(b)分別給出了全波法得到的電場、磁場分量二維仿真圖. 由于本節(jié)引入的是TM平面波垂直入射電離層，電場、磁場垂直分量很小. 隨著透入電離層的高度增加，電波方向發(fā)生改變，80 km以上高度出現(xiàn)了電場z分量.圖2為兩種方法得到的電離層能量吸收隨高度的變化. 在85 km以下高度，能量吸收差別很?。?5 km以上區(qū)域，兩種方法得到的能流吸收峰值都在90 km高度附近. 不同的是，由于全波法考慮到電離層各分層中反射、透射波相互疊加的情況，在進入電離層能流相同的情況下，電離層能流吸收要大于Rodriguez方法得到的能流吸收.

圖1 電磁場分量隨高度的變化Fig. 1 Variation of electric field and magnetic field component with height

圖2 電離層能流吸收對比Fig. 2 Comparison of energy flow absorption in the ionosphere

圖3給出了Rodriguez方法和全波計算方法得到的電子溫度擾動對比. 與圖2中的能流吸收情況相對應(yīng)，Rodriguez方法中電子溫度的變化主要集中在70～95 km高度范圍內(nèi)，電子溫度在85 km和90 km高度附近存在兩個峰值點，最高可以達到290 K. 用全波法得到的仿真結(jié)果中由于電離層能流吸收集中在90 km高度附近，電子溫度的擾動范圍主要集中在75～100 km且僅在90 km處有一個峰值點，最高可達320 K. 可以看出，波模的反射迭加不僅提高了電離層的能流吸收，電子溫度擾動幅度也隨之增大.

圖3 電子溫度擾動對比Fig. 3 Comparison of electronic temperature disturbance

電子溫度的擾動導(dǎo)致電子濃度發(fā)生變化，圖4給出了電子濃度與初始時刻的比值，受背景電子濃度分布的影響，擾動的幅度和高度與電子溫度存在一定差異. Rodriguez方法仿真的結(jié)果中，電子濃度在83 km處損耗最大，降低了約12%. 用全波法得到的電子濃度模型中，電子在85 km損耗最大，由于85 km處的背景電子濃度較高，濃度損失反而不如Rodriguez方法的結(jié)果明顯，降低了約10%. 在93 km以上區(qū)域，由于電子溫度的提高降低了電子與正離子的復(fù)合效率，該區(qū)域內(nèi)的電子濃度反而有所提高.

圖4 電子濃度擾動對比Fig. 4 Comparison of electron concentration disturbance

2.2 不同入射波、地磁場條件下的擾動結(jié)果對比

2.2.1 入射角變化對電離層擾動的影響

將磁傾角設(shè)為60°，地磁場水平分量指向正北方向，在入射波水平方向的分量與正北方向夾角φ為180°的情況下，本節(jié)對比了入射角 θ對電離層擾動的影響. 圖5(a)、(b)分別給出了各高度下電場、磁場矢量隨入射角的變化，可以看出，電場的衰減主要發(fā)生在90 km以上高度，磁場的衰減主要發(fā)生在85 km以上高度. 圖5(c)為不同入射角下的電波能流對比，可以看出，90 km高度以下區(qū)域，入射角θ = 30°的情況下電波能流密度最大，反射透射波的疊加使80～90 km高度范圍內(nèi)的能流密度大于入射波能流密度.90 km高度以上區(qū)域，入射角θ=0°時，能流密度最大，這是因為VLF波垂直入射電離層時穿透性最強，能流密度隨高度衰減最慢.

圖6給出了不同入射角對應(yīng)的電離層能流吸收情況，電離層能流吸收主要集中在80 km以上區(qū)域，隨著入射角增大，能流吸收峰值所在高度逐漸從92 km向85 km偏移，由于電波在反射面附近能流吸收最大，進一步證明了改變?nèi)肷浣菚绊戨姴ǖ姆瓷洹⑽崭叨? 入射角越大的情況下，電波穿透性越弱，反射吸收面高度越低.

圖5 電離層參數(shù)隨入射角的變化Fig. 5 Variation of ionospheric parameters with incident angles

圖6 不同入射角下的電離層能流吸收Fig. 6 Ionospheric energy flow absorption at different incident angles

與電離層能流吸收相對應(yīng)，圖7給出了電子溫度的擾動情況，電子溫度變化主要集中在75～95 km高度范圍內(nèi). 由式(12)可知，電子溫度的變化除受電離層能量變化影響外，還與背景電子濃度有關(guān)，電子濃度越低的區(qū)域電子溫度變化越明顯. 在入射角為30°和60°的情況下，電子溫度擾動峰值都在85 km高度附近，θ= 30°時擾動幅度大一些，溫度最高可達550 K，由于90 km高度電子濃度比85 km大得多，導(dǎo)致θ= 0時電子溫度的擾動峰值很小，溫度最高僅為300 K.

圖8給出了不同入射角下的電子濃度與初始時刻的比值，93 km以下區(qū)域，入射角 θ越大，電子濃度主要擾動區(qū)域所在高度越低，電子濃度損失越大. θ=60°情況下，電子濃度擾動幅度最大，83 km附近電子濃度僅為初始時刻濃度的80%，θ= 30°情況下，電子濃度擾動幅度最小，85 km附近電子濃度僅損失15%. 93 km以上區(qū)域，電子濃度因電子與正離子復(fù)合效率降低而提高.

圖7 不同入射角下的電子溫度Fig. 7 Electron temperatures at different incident angles

圖8 不同入射角下的電子濃度與初始時刻的比值Fig. 8 Ratio of electron concentration to initial time at different incident angles

2.2.2 入射波方位角對電離層擾動的影響

在入射角θ= 30°，地磁場方向為θB=60°、φB= 0°的情況下，圖9(a)、(b)分別給出了不同方位角下的電場與磁場分布. 電場的衰減主要發(fā)生在90 km以上高度，磁場的衰減主要發(fā)生在85 km以上高度. 在本文給出的初始條件下，場強隨方位角增加而增大.圖9(c)給出了不同入射波方位角φ下電離層能流吸收隨高度的變化. 電離層對電波能流的吸收主要集中在75 km以上區(qū)域，改變方位角會影響電離層的能流吸收幅度，φ值越大，吸收峰值越大. 峰值點所在高度不受φ值變化的影響，集中在90 km高度附近，說明方位角變化不會改變電波的反射吸收高度.

圖10給出了不同φ值下的電子溫度擾動情況，溫度的擾動主要發(fā)生在75 km以上的電離層區(qū)域.受背景電子濃度影響，電子溫度的峰值點出現(xiàn)在87 km高度附近且位置不隨φ值變化. 在φ為135°的情況下電子溫度提升最明顯，最高可達350 K；φ為45°時電子溫度擾動幅度最小，其溫度峰值僅為300 K.

圖9 電離層參數(shù)隨方位角的變化Fig. 9 Variation of ionospheric parameters with azimuth angles

圖11給出了不同φ值下的電子濃度擾動情況.可以看出，電子濃度擾動主要集中在85 km附近，φ為135°情況下，85 km處電子濃度損失最大，濃度僅有初始時刻的87%；φ取45°情況下，電子濃度擾動幅度最低，85 km處電子濃度僅損失10%.

圖10 不同方位角下的電子溫度Fig. 10 Electronic temperature at different azimuths

圖11 不同方位角下的電子濃度與初始時刻的比值Fig. 11 Ratio of electron concentration to initial time at different azimuths

2.2.3 磁傾角對電離層擾動的影響

受地理位置影響，不同緯度地區(qū)磁傾角存在較大差異，在入射波垂直向上的情況下，本節(jié)分析了磁傾角θB分別為30°、45°、60°時，基于全波法的電離層擾動情況. 圖12(a)、(b)分別給出了不同磁傾角下的電場強度與磁場強度對比. 在場強衰減比較大的90 km以上高度，二者夾角越大，電場與磁場的衰減越小. 圖12(c)給出了θB不同的情況下電波的能流密度. 我們發(fā)現(xiàn)，θB的變化對能流密度的影響主要集中在80 km以上的電離層區(qū)域. 在80～90 km高度范圍內(nèi)，能流密度隨θB的增大而降低，特別是θB= 30°的情況下，反射、透射波相互疊加使該區(qū)域內(nèi)能流密度提升明顯. 90 km以上區(qū)域，θB越大電波能流衰減越慢.

圖13給出了不同θB下 的電離層能流吸收情況，電波能流的吸收主要集中在80～100 km高度范圍內(nèi)，三種夾角下的能流吸收峰值點都在90 km高度附近. θB越 大，電離層能流吸收峰值越小，這是因為當(dāng)來波方向與地磁場夾角越接近90°，波的衰減越大，電離層能流吸收越大.

圖12 電離層參數(shù)隨磁傾角的變化Fig. 12 Variation of ionospheric parameters with magnetic dip angles

圖13 不同磁傾角下的電離層能流吸收Fig. 13 Ionospheric energy flow absorption at different magnetic dip angles

圖14給出了不同θB下的電子溫度模型，與圖8中的電離層能流吸收模型相對應(yīng)，電子溫度的擾動范圍主要集中在80～95 km高度內(nèi)，θB越小，電子溫度提升越明顯，θB= 30°的情況下，電子溫度最高可達600 K. 不同磁傾角下，電子溫度極大值點所在高度也有差異，增大θB，溫度擾動的峰值點所在高度逐漸從85 km向90 km偏移.

圖14 不同磁傾角下的電子溫度Fig. 14 Electron temperature at different magnetic inclinations

電離層電子溫度的變化會引起電子與其他粒子的附著、分離效率發(fā)生改變，圖15給出了不同θB下的電子濃度變化. 對比可知，80～92 km高度內(nèi)，電子溫度的提高增強了中性粒子對電子的吸附效率，導(dǎo)致電子濃度降低，負離子濃度提高，且θB越小，電子損耗越明顯，θB=30°的情況下，85 km高度處電子濃度只有初始時刻的80%. 92 km以上高度，由于電子溫度的增加抑制了電子與正離子的復(fù)合效率，電子濃度由降低轉(zhuǎn)變?yōu)樘岣?

圖15 不同磁傾角下的電子濃度與初始時刻的比值Fig. 15 Ratio of electron concentration to initial time at different magnetic inclinations

2.3 不同電離層背景參量下的擾動結(jié)果對比

由于中性成分參數(shù)、化學(xué)反應(yīng)系數(shù)等參量控制著電離層碰撞頻率、能量損失和化學(xué)反應(yīng)效率的變化，改變電離層背景參量亦會對電離層擾動情況產(chǎn)生影響. 本節(jié)分析了在入射波方向θ= 30°、φ= 180°，頻率f= 10 kHz，地磁場方向θB= 60°、φB= 0°，強度為5 0μT情況下，不同背景參量下的電離層擾動情況.

2.3.1 化學(xué)反應(yīng)系數(shù)對電離層擾動的影響

100 km以下的低電離層區(qū)域，中性成分主要為氮氣和氧氣，中性氣體對電子的吸附是引起電子損耗的重要因素. 圖16給出了中性粒子的電子附著率β對電子濃度的影響. 讀圖可知，隨著電子附著到中性粒子的效率提高，電子濃度向減小的趨勢發(fā)展. 電子濃度的變化進一步影響到電子溫度的擾動情況，由圖17可知，隨著電子濃度的降低，電子溫度在75～90 km高度范圍內(nèi)的擾動幅度逐漸提高.

圖16 不同β下的電子濃度與初始時刻的比值Fig. 16 Ratio of electron concentration to initial time at different values of β

圖17 不同β下的電子溫度Fig. 17 Electron temperature at different values of β

2.3.2 中性成分濃度對電離層擾動的影響

圖18給出了中性成分濃度Nn分別為原值N0的0.6倍、1倍、1.4倍情況下，電子溫度的擾動情況.結(jié)合式(12)和文獻[15]給出的能量損失模型可知，Nn越大，能量損失項Li越大，電子溫度的擾動幅度越小. 中性成分濃度的變化不會影響電子溫度峰值點所在高度，電子溫度的峰值點都出現(xiàn)在88 km附近.

圖18 不同中性成分濃度下的電子溫度Fig. 18 Electron temperature at different concentrations of neutral components

圖19給出了不同Nn下的電子濃度擾動情況，Nn的增大同時提高了中性粒子對電子的吸附效率和電子從負離子上的脫離效率，但是前者的提升更為明顯. 與圖16類似，電子濃度隨Nn的提高向變小的趨勢發(fā)展. 高度越低，背景電子濃度越小，電子濃度的變化越明顯.

圖19 不同中性成分濃度下的電子濃度與初始時刻的比值Fig. 19 Ratio of electron concentration to initial time at different concentrations of neutral components

2.3.3 中性成分溫度對電離層擾動的影響

圖20給出了中性成分溫度Tn分別為原值T0的0.8倍、1倍、1.2倍情況下，電子溫度的擾動情況.與圖18相反，Tn越大電子溫度的擾動幅度越大.Tn變?yōu)樵?.8倍的情況下，電子溫度最大值約為325 K；Tn為原值1.2倍的情況下，電子溫度最大值可提高到425 K.

圖20 不同中性成分溫度下的電子溫度Fig. 20 Electron temperature at different neutral component temperatures

圖21給出了不同Tn下的電子濃度擾動情況，可以看出在92 km以下高度，由于Tn的增大提高了中性成分對電子的吸附效率，電子濃度隨Tn的增大向變小的趨勢發(fā)展. 92 km以上高度，正離子濃度的增大使電子與正離子的復(fù)合效率成為影響電子濃度的重要因素，由于電子與正離子的復(fù)合效率隨Tn的增大而減小，負離子上的電子脫離率隨負離子濃度的提高而增大，電子濃度隨Tn的增大向變大的趨勢發(fā)展.

圖21 不同中性成分濃度下的電子溫度與初始時刻的比值Fig. 21 Ratio of electron concentration to initial time at different neutral component temperatures

3 結(jié) 論

通過全波法和波模振幅遞歸算法，結(jié)合Glukhov-Pasko-Inan低電離層化學(xué)反應(yīng)模型，對VLF波影響下的電離層電子濃度、溫度的擾動情況進行仿真. 主要成果與結(jié)論如下：

1)本文使用的全波計算方法與其他文獻方法的差別在于全波法考慮到上、下行波的疊加情況. 與Rodriguez方法對比發(fā)現(xiàn)，在進入電離層能流相同的情況下，兩種方法得到的電離層能流吸收都集中在75～100 km高度范圍內(nèi)，峰值點都出現(xiàn)在90 km高度附近. 用全波法算出的電離層能流吸收大于Rodriguez方法得到的能流吸收，進而造成電子溫度、濃度擾動的不同.

2)對比了不同電波入射角、方位角、磁傾角下的電離層擾動情況，發(fā)現(xiàn)增大平面波入射角會導(dǎo)致電波能流的反射吸收高度發(fā)生改變，入射角越大，電波的反射吸收高度越低. 方位角、磁傾角的變化不會改變電波反射高度，但會對電離層能流吸收幅度造成影響，入射波與地磁場夾角越接近90°，電離層能流吸收的峰值越高. 由于電離層中電子的擾動情況還與電離層背景參量、化學(xué)反應(yīng)過程有關(guān)，在多種因素調(diào)控下，電子溫度、濃度的擾動幅度和擾動范圍存在明顯差異.

3)對比了中性成分濃度、溫度等背景參量對電離層擾動的影響. 發(fā)現(xiàn)上述參量主要通過改變電離層能量損失項L控制電子溫度的擾動幅度. 背景參量的變化還會改變離子間的化學(xué)反應(yīng)效率，進而影響電子濃度的擾動情況. 不同之處在于，中性成分濃度越高，電子溫度擾動幅度越小，電子濃度的損耗越大. 中性成分溫度越高，電子溫度擾動幅度越大，92 km以下高度電子濃度損耗越大，92 km以上高度電子濃度損耗越小.