趙珍珠 劉江凡 蒲玉蓉 張金生 席曉莉*

（1.西安理工大學(xué),西安 710048；2.火箭軍工程大學(xué),西安 710025）

引言

羅蘭C 系統(tǒng)是陸基無線電導(dǎo)航授時系統(tǒng)，具備較強(qiáng)的抗干擾能力，可彌補(bǔ)GNSS 的脆弱性，是GNSS 重要的補(bǔ)充和備份[1-2].羅蘭C 導(dǎo)航系統(tǒng)工作在長波頻段，電波以地波和天波兩種方式傳播.地波沿地球表面繞射傳播至接收點(diǎn)，是當(dāng)前羅蘭C 接收機(jī)授時導(dǎo)航中主要應(yīng)用的信號，其覆蓋范圍海上約為1 500 km，陸地只有1 200 km 或更小.而天波經(jīng)電離層反射到達(dá)接收點(diǎn)，傳播距離較遠(yuǎn).

電離層復(fù)雜的時空變化是影響低頻天波傳播特性的關(guān)鍵因素，制約著利用天波進(jìn)行導(dǎo)航定位的應(yīng)用實(shí)現(xiàn).1980 年，陜西天文臺的苗永瑞對天波授時的可行性進(jìn)行了充分討論[3]，并給出一跳天波場強(qiáng)、相位預(yù)測模型.隨后大量關(guān)于羅蘭信號天波傳播的研究開始展開，包括地球模型[4]、季節(jié)變化[5-6]、太陽活動[7]等對低頻天波傳播影響因素的研究，以及基于實(shí)測數(shù)據(jù)對電離層分布模型的反演[8-9]等.2004 年，潘威炎出版了專著《長波超長波極長波傳播》[10].近些年，隨著電磁理論基礎(chǔ)的研究深入，Cummer 率先提出在長波傳播研究中可使用時域有限差分法(finite-difference time-domain，F(xiàn)DTD)[11]，開啟了長波傳播研究的新思路.隨后，文獻(xiàn)[12-16]總結(jié)了FDTD在長波信號傳播中的應(yīng)用并從計算方法、參數(shù)選擇、實(shí)測改進(jìn)等方面對其進(jìn)行改進(jìn)創(chuàng)新.2019 年，周麗麗等人將 FDTD 方法應(yīng)用于羅蘭系統(tǒng)天波傳播研究中[17]，隨后仿真計算了地-電離層波導(dǎo)中的多跳天波，對多跳天波帶來的多徑時延進(jìn)行了分析預(yù)測[18].然而，在這些研究中總是將電離層考慮為完全反射面或指數(shù)形式，未能考慮其自身復(fù)雜性(時空變化、多層反射等)對低頻天波信號的影響.

國際參考電離層(international reference ionosphere,IRI)模型能夠提供任意位置、任意時刻的電子密度分布.同時，對電波傳播影響較大的不僅有電子密度分布，還有碰撞頻率，一般研究中考慮為常數(shù)或者指數(shù)分布.而結(jié)合IRI 模型和NRLMSISE-00 大氣模型，即可得到任意位置、任意時間的各粒子的濃度分布信息，從而計算得到碰撞頻率分布[19].本文基于這兩種模型以及雙線性變換FDTD(bilinear transform FDTD,BT-FDTD)[20]方法，仿真計算了羅蘭C 天波經(jīng)電離層反射的反射信號，分析了反射信號幅度、時延、多徑隨入射角度、一天中的時間以及季節(jié)的變化.仿真結(jié)果表明，對于TE 波，反射信號強(qiáng)度隨入射角度的增大先減小后增大；一天中，其他條件不變，經(jīng)電離層反射的羅蘭C 天波信號幅度相差最大可達(dá)32.22 dB，時延差可達(dá)69.03 μs；選取每個月同一天同一時刻的電離層模型，其他條件不變，信號幅度差可達(dá)20.97 dB，時延差可達(dá)6.49 μs.低頻天波經(jīng)電離層反射信號一般包含一個完整的羅蘭C 信號，在晝夜過渡時會出現(xiàn)兩個信號.本文研究可為羅蘭C 采集信號分析及接收機(jī)的研制提供理論支持.

1 理論計算模型

電離層是影響羅蘭C 天波信號特性的關(guān)鍵因素.本文結(jié)合IRI 模型和NRLMSISE-00 大氣模型來構(gòu)建電波傳播信道模型，將電離層人為分層后，低頻天波傳播問題則可看作分層半空間介質(zhì)中平面波的傳播問題.因此，采用分層半空間中平面波傳播的準(zhǔn)一維BT-FDTD 方法計算電波經(jīng)電離層反射的反射信號.下面對模型及計算方法進(jìn)行介紹.

1.1 IRI 模型

IRI 是國際無線電科學(xué)聯(lián)合會(International Union of Radio Science,URSI)根據(jù)全球地面觀測站所得的大量電離層觀測數(shù)據(jù)和多年來電離層模型的理論研究成果，在全球范圍內(nèi)普遍適用.本文采用目前最新版本IRI-2016 模型.時間采用世界時，地理位置信息為：108.999 855°E，34.259 563 4°N.圖1 為2020-04-02的電子密度分布，圖2 為取每月2 日06:00 得到的電子密度分布年變化.由圖1 可以看到，電離層電子密度在1 d 內(nèi)呈現(xiàn)較大的變化，每立方米相差達(dá)到3 個量級.以所選仿真日期當(dāng)天來看，在05:00 達(dá)到最大值，17:00 達(dá)到最小值，在11:00、23:00 時間段變化劇烈.由圖2 可以看到，電子密度年變化在最底層較為明顯，整體變化范圍沒有一天中的電子密度變化范圍大.

圖1 2020-04-02 電子密度分布Fig.1 The electron density distribution on April 2,2020

圖2 2020 年電離層電子密度隨月份的變化Fig.2 Variation of ionospheric electron density with month in 2020

1.2 碰撞頻率模型

電離層碰撞頻率的計算公式為[19]

通過IRI-2016 模型，可得到電子密度和電子溫度數(shù)據(jù)；結(jié)合NRLMSISE-00 大氣模型，可得到氮?dú)?、氧氣和氧原子的密度，最后即可得到電離層碰撞頻率.

圖3 為2020-04-02 的碰撞頻率分布.圖4 為取每個月2 日06:00 得到的碰撞頻率分布年變化.結(jié)合圖1～4 可以看出，雖然電子密度模型晝夜變化較大，但是碰撞頻率變化不大.

圖3 2020-04-02 碰撞頻率分布Fig.3 Collision frequency distribution on April 2,2020

圖4 2020 年碰撞頻率隨月份的變化Fig.4 Changes in the frequency of collisions by month in 2020

1.3 BT-FDTD 方法

麥克斯韋旋度方程為

本構(gòu)方程

非磁化電離層介質(zhì)中的復(fù)介電常數(shù)

式中：ωp為電離層電子振蕩頻率；ν為電離層碰撞頻率；ω為電磁波頻率.

將式(5)代入到式(4)中，有

采用雙線性變換，令

將式(7)、式(8)代入到式(6)中，得

其他場量之間的關(guān)系同F(xiàn)DTD 方法.

1.4 斜入射TE 波的準(zhǔn)一維BT-FDTD

麥克斯韋方程的TE 波頻域形式為

式(10)中前兩式的時域形式為

離散式(11)可得到由Hz→Ex和Hz→Ey的時域步進(jìn)計算式.

根據(jù)分層半空間特點(diǎn)，推導(dǎo)式(10)中第三式的時域形式.設(shè)入射波為 exp(-jkxx-jkyy)，將式(10)中的第二式對x求導(dǎo)，并利用平面波算子替換關(guān)系?/?x→-jkx=-jk0sin θ，得到

由于波矢量的切向分量連續(xù)，則kx=k0sin θ 在各個分層均適用，將式(12)代入式(10)中的第三式并記第l層的介電系數(shù)為 εlr，則有

則式(13)變?yōu)?/p>

過渡到時域形式為

對離散式(16)進(jìn)行FDTD 可得到Ex→的時域遞推計算式.對于半空間平面波斜入射TE 情形，時域方程為式(11)和式(16)，可按照1.3 節(jié)中的BT-FDTD方法進(jìn)行計算.這些方程類似于一維波方程，因而稱為分層半空間中斜入射平面波傳播的TE 準(zhǔn)一維波方程[21].

最后，推導(dǎo)TE 波的準(zhǔn)一維方程中由磁場輔助量到磁場的遞推公式.

由式(14)得

將式(5)代入到式(17)中，并采用雙線性變換得

1.5 方法驗(yàn)證

在1.4 節(jié)中，對分層半空間中斜入射TE 波的準(zhǔn)一維BT-FDTD 方法進(jìn)行了公式推導(dǎo).本節(jié)通過一個計算算例，與混合矩陣法(hybrid matrix method,HMM)[22-23]進(jìn)行比較，驗(yàn)證算法的正確性.仿真計算模型如圖5 所示，電離層高度為60～160 km，每10 m取一個值，共10 000 層；準(zhǔn)一維BT-FDTD 方法總網(wǎng)格大小為 200×10 300，200～10 200 層為電離層，計算網(wǎng)格大小為1 0 m×10 m，時間步長為1.67×10-8s，吸收邊界80 層，源位于(100,120)網(wǎng)格處，入射角為θ，采樣點(diǎn)位于(100，200).

圖5 仿真計算模型Fig.5 Simulation calculation model

算例模型參數(shù)為：電離層模型選取2020-04-02 T06:00 的電子密度分布模型和碰撞頻率模型，入射角度為 θ=45°.兩種方法計算得到的反射系數(shù)隨入射波頻率的變化如圖6 所示.可以看到，電波斜入射時，兩種方法得到的反射系數(shù)完全吻合，所采用的準(zhǔn)一維BT-FDTD 方法用于計算分層半空間中低頻波的反射信號是正確的.

圖6 2020-04-02T06:00 反射系數(shù)隨入射波頻率的變化Fig.6 Change of reflection coefficient with the frequency of incident wave at 6:00 on April 2,2020

2 信號反射特性分析

采用上述模型及方法，對羅蘭C 信號經(jīng)電離層反射傳播中的信號反射特性進(jìn)行分析.電離層范圍選取及準(zhǔn)一維BT-FDTD 網(wǎng)格參數(shù)同1.5 節(jié).

2.1 反射系數(shù)隨入射角度的變化

計算圖1～4 中06:00、18:00 模型參數(shù)下反射系數(shù)隨入射角度的變化，結(jié)果如圖7 所示.可以看到夜間(18:00)電離層反射系數(shù)幅度大于白天(06:00)，天波信號更強(qiáng)；隨著入射角的增大，反射系數(shù)先減小再增大；不管白天還是夜間，入射角度大時，反射信號強(qiáng)度強(qiáng)于小角度入射.因此，天地波復(fù)合定位中，近距離時地波信號較強(qiáng)，采用地波定位；當(dāng)距離較遠(yuǎn)時，采用天波定位.

圖7 2020-04-02 反射系數(shù)隨入射角的變化Fig.7 Change of reflection coefficient with incident angle on April 2,2020

2.2 隨時間的變化

選取θ=60°計算一天、一年當(dāng)中反射信號的幅度和時延.此處選取反射系數(shù)中100 kHz 頻點(diǎn)的反射信號幅度；因直接采用相位計算時延存在相位模糊度問題，此處采用FFT/IFFT 頻譜相除法[24]對羅蘭C 信號過電離層的時延進(jìn)行計算.一天中反射信號幅度、時延計算結(jié)果如表1 所示.可以看到：一天之中，其他條件不變，僅通過電離層反射，幅度差可達(dá)32.22 dB，時延差可達(dá)69.03 μs；在晝夜過渡期間差異較大，其他時間較為平穩(wěn).一年中每個月同一天同一時刻反射信號幅度、時延計算結(jié)果如表2 所示.可以看到，一年之中，其他條件不變，僅通過電離層反射，幅度差可達(dá)20.97 dB，時延差可達(dá)6.49 μs.

表1 2020-04-02 各個時刻反射信號幅度和時延Tab.1 Reflected signal amplitude and time delay at various times on April 2,2020

表2 2020 年每月2 日06:00 反射信號幅度和時延Tab.2 Reflected signal amplitude and time delay at 6:00 on the 2nd of each month in 2020

2.3 反射信號多徑分析

在對2020-04-02 一天中信號時延計算中，發(fā)現(xiàn)當(dāng)出現(xiàn)晝夜交替電離層電子密度變化劇烈期間，會出現(xiàn)不止一個羅蘭C 信號，如圖8 所示，尤其23:00多徑信號幅值較大.因此，當(dāng)接收機(jī)工作在這個時段時，要考慮這個多徑信號的影響；在對實(shí)際天波采集信號中的多徑信號來源進(jìn)行分析時，也可結(jié)合實(shí)際電離層模型對反射信號進(jìn)行分析，考慮采集時間段內(nèi)來自電離層本身的多徑.

圖8 2020-04-02 信號FFT/IFFT 頻譜相除法分析結(jié)果Fig.8 Signal analysis result by FFT/IFFT spectrum division method on April 2,2020

3 結(jié)論

本文針對低頻天波經(jīng)電離層反射傳播反射信號特性進(jìn)行分析.一天中，其他條件不變，經(jīng)電離層反射的羅蘭C 低頻天波信號幅度相差最大可達(dá)32.22 dB，時延差可達(dá)69.03 μs；一年之中，其他條件不變，僅通過電離層反射，幅度差可達(dá)20.97 dB，時延差可達(dá)6.49 μs.低頻天波經(jīng)電離層反射信號一般包含一個完整的羅蘭C 信號，在晝夜過渡時會出現(xiàn)兩個信號.本文的分析對羅蘭C 實(shí)際采集信號分析以及接收機(jī)的研制具有積極的意義.