唐新杰 劉默然 喬建東 周晨

（武漢大學(xué)電子信息學(xué)院, 武漢 430072）

0 引 言

無線電波在大氣層中傳播時(shí)，由于大氣介質(zhì)不均勻分布的影響產(chǎn)生折射效應(yīng)，使得電波射線發(fā)生彎曲、傳播速度小于光速，導(dǎo)致雷達(dá)等測(cè)控系統(tǒng)在定位、導(dǎo)航時(shí)產(chǎn)生一定的偏差[1]。隨著技術(shù)的發(fā)展和各種應(yīng)用場(chǎng)景的需要，大氣電波折射誤差逐漸成為制約航天測(cè)定軌、空間目標(biāo)監(jiān)視等遠(yuǎn)程探測(cè)系統(tǒng)精度的重要因素，因此需要對(duì)測(cè)控系統(tǒng)探測(cè)目標(biāo)測(cè)量值的折射誤差予以修正。

現(xiàn)有的雷達(dá)折射誤差修正方法可以分為射線描跡法、近似修正法和新型修正法三種，其中射線描跡法因?yàn)槠渚雀?、?shí)用性強(qiáng)等優(yōu)點(diǎn)在研究和工程實(shí)踐中都有廣泛的應(yīng)用[2]。折射誤差修正精度不僅和選用的算法有關(guān)，還取決于選取的大氣折射率剖面的準(zhǔn)確度，學(xué)者們對(duì)此進(jìn)行了大量的研究，通過模型或一些探測(cè)方法對(duì)大氣折射率剖面進(jìn)行建模。王寧等[3]引入徑向基函數(shù)(radial basis function, RBF)神經(jīng)網(wǎng)絡(luò)算法反演大氣折射率，驗(yàn)證了利用微波輻射計(jì)對(duì)流層電波折射進(jìn)行修正的可行性；徐艷等[4]對(duì)中近程雷達(dá)應(yīng)用場(chǎng)景中高精度目標(biāo)定位誤差進(jìn)行了研究分析；張瑜等[5]針對(duì)下墊面復(fù)雜地區(qū)的三維大氣結(jié)構(gòu)，采用直接探測(cè)法和全國(guó)大氣剖面模型組合的方法，有效提高了電波射線折射率分布的準(zhǔn)確度；喬江等[6]對(duì)比研究了折射率剖面模型和映射函數(shù)在對(duì)流層折射誤差修正中的精度表現(xiàn)；李川川等[7]結(jié)合微波輻射計(jì)和GNSS接收機(jī)研制的電波折射誤差修正設(shè)備，實(shí)現(xiàn)了S波段雷達(dá)實(shí)時(shí)、高精度折射誤差修正；劉琨等[8]提出了適用于高軌道目標(biāo)折射誤差修正的空間投影法和自適應(yīng)網(wǎng)格法。

對(duì)于航天目標(biāo)的雷達(dá)探測(cè)，電波傳播主要受對(duì)流層和電離層兩方面的影響。而對(duì)于P波段雷達(dá)而言，由于其頻率相對(duì)較低，受電離層的影響更為顯著，也明顯比對(duì)流層影響更為嚴(yán)重。同時(shí)，低仰角相對(duì)于高仰角的誤差更大，而低仰角意味著水平距離更遠(yuǎn)，其環(huán)境參數(shù)與探測(cè)點(diǎn)的環(huán)境參數(shù)差別也更大，因此對(duì)低仰角的誤差修正就更為困難，這也一直是相關(guān)領(lǐng)域的難點(diǎn)之一。本文結(jié)合NRLMSISE-00大氣模型對(duì)海拔60 km以下的對(duì)流層折射率剖面進(jìn)行建模，并利用TEC數(shù)據(jù)同化的方式建立電離層現(xiàn)報(bào)模型，得到更精確的電離層實(shí)時(shí)信息，完成電離層折射率剖面的建模。在此基礎(chǔ)上，利用射線描跡法對(duì)P波段雷達(dá)探測(cè)數(shù)據(jù)進(jìn)行修正。與基于國(guó)際參考電離層（International Reference Ionosphere, IRI）模型的電離層折射率剖面相比，利用TEC數(shù)據(jù)同化現(xiàn)報(bào)模型得到的電離層折射率剖面可以極大提高修正精度。

1 建模方法

1.1 折射誤差計(jì)算方法

射線描跡法從斯奈爾定律和費(fèi)馬原理出發(fā)，基于數(shù)學(xué)幾何關(guān)系嚴(yán)格推導(dǎo)而來，具有較高的精度，因此成為工程中廣泛應(yīng)用的修正方法[9]。由于精度高，該方法也常被用來與其他的修正模型進(jìn)行對(duì)比驗(yàn)證。圖1給出了球面分層電波傳播軌跡示意圖，C為地心，O為測(cè)控站位置，T為目標(biāo)位置。

圖1 球面分層電波射線軌跡幾何圖形Fig.1 Geometry of spherical layered radio ray

高軌道遠(yuǎn)程目標(biāo)一般沉浸于電離層之中的空間環(huán)境，由斯奈爾定理可以推導(dǎo)出雷達(dá)探測(cè)獲得的目標(biāo)視在距離Re計(jì)算式為

式中：rI為電離層底的高度；r0為測(cè)控站處的地心半徑，r0=a+h0，a為地球半徑，h0為測(cè)控站海拔高度；n為電波射線軌跡上變化的折射指數(shù)；r為目標(biāo)至地心之間的距離；n0為地表處的折射指數(shù)；θ0為目標(biāo)視在仰角。

測(cè)控站與目標(biāo)之間地心張角φ的計(jì)算式為

目標(biāo)真實(shí)仰角α0的計(jì)算式為

目標(biāo)真實(shí)距離R0的計(jì)算式為

由上述式子計(jì)算得到距離折射誤差?R和仰角折射誤差ε分別為

1.2 基于電離層現(xiàn)報(bào)同化模型的電離層折射率剖面

數(shù)據(jù)同化按照理論基礎(chǔ)分類主要分為兩類：一類是基于統(tǒng)計(jì)學(xué)的估計(jì)理論，如最優(yōu)插值、卡爾曼濾波、集合卡爾曼等；一類是基于變分理論的方法[10]。卡爾曼濾波算法是常見的線性濾波和預(yù)測(cè)方法中的一種，是由匈牙利科學(xué)家Rudolph E.Kalman在1960年提出的。它是最早出現(xiàn)的順序同化算法，一般認(rèn)為它是順序同化算法的理論基礎(chǔ)。在近年的研究中，卡爾曼濾波算法越來越多地被應(yīng)用于電離層數(shù)據(jù)同化中。

根據(jù)卡爾曼濾波理論，對(duì)觀測(cè)資料的基本同化過程可表示為：

式中：為t+1時(shí)刻的背景場(chǎng)，由背景電離層模型給出；M為狀態(tài)轉(zhuǎn)換矩陣；為t時(shí)刻的分析場(chǎng)，即同化后得到的電子密度；為t+1時(shí)刻的背景場(chǎng)誤差協(xié)方差矩陣；為t時(shí)刻的分析場(chǎng)誤差協(xié)方差矩陣；Q為模型誤差方差矩陣；為t+1時(shí)刻的分析場(chǎng)；K為增益矩陣，起到觀測(cè)數(shù)據(jù)對(duì)背景場(chǎng)的調(diào)整作用；Y為觀測(cè)場(chǎng)；H為觀測(cè)算子，代表觀測(cè)量與模型參量之間的關(guān)系；R為觀測(cè)場(chǎng)誤差協(xié)方差矩陣；為t+1時(shí)刻的分析場(chǎng)誤差協(xié)方差矩陣。

本文采用目前被廣泛應(yīng)用于電離層研究的經(jīng)驗(yàn)?zāi)Ｐ?IRI模型作為背景模型，版本為IRI2016。IRI由國(guó)際無線電科學(xué)聯(lián)盟和空間研究委員會(huì)共同研發(fā)和維護(hù)，其數(shù)據(jù)來源于全球180多個(gè)電離層探測(cè)站的資料以及衛(wèi)星的資料，主要模擬了電子密度、離子成分、離子溫度、電子溫度等一系列電離層參量的全球分布[11]。背景場(chǎng)的經(jīng)緯度范圍選取以雷達(dá)測(cè)站為中心的探測(cè)覆蓋范圍及周邊區(qū)域，緯度為15°N~55°N，經(jīng)度為70°E~140°E，水平格點(diǎn)分辨率設(shè)定為1°×1°，高度范圍為60~1 500 km， 步長(zhǎng)為10 km，因此背景場(chǎng)誤差協(xié)方差矩陣Pb的大小為419 184 ×419 184。參與三維數(shù)據(jù)同化的觀測(cè)數(shù)據(jù)主要為中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)(Crustal Movement Observation Network of China, CMONOC)的261個(gè)GNSS監(jiān)測(cè)臺(tái)站提供的電離層TEC數(shù)據(jù)，GNSS接收機(jī)的采樣間隔為30 s。TEC是GNSS信號(hào)傳播路徑對(duì)路徑上電子密度的線積分，所以觀測(cè)算子H即為GNSS信號(hào)傳播路徑在電離層?xùn)鸥裰写┰降拈L(zhǎng)度。

基于統(tǒng)計(jì)的協(xié)方差模型在估計(jì)背景誤差協(xié)方差矩陣時(shí)會(huì)引入虛假相關(guān)，即在物理上不相關(guān)或空間上距離較遠(yuǎn)的狀態(tài)變量之間存在虛假相關(guān)。為了降低虛假相關(guān)對(duì)同化過程的影響，提升同化效率，本文應(yīng)用了協(xié)方差局地化(covariance localization，CL)方法[12]。CL方法采用一個(gè)基于距離的相關(guān)系數(shù)ρ與背景誤差協(xié)方差矩陣作舒爾積，來替代原有的背景協(xié)方差矩陣。相關(guān)函數(shù)采用常用的GC (Gaspari-Cogn)五階分段多項(xiàng)式函數(shù)，表達(dá)式如下：

式中：z為物理空間中兩個(gè)格點(diǎn)之間的距離；c為長(zhǎng)度尺度。

1.3 基于空間格點(diǎn)化獲取的射線路徑上折射率剖面

傳統(tǒng)的射線描跡法忽略水平方向的不均勻性，往往采用垂直方向上的空間環(huán)境剖面來代表不同仰角下觀測(cè)條件，即傳統(tǒng)的射線描跡法適用于下墊面平坦、大氣水平相對(duì)均勻的地區(qū)，而對(duì)于高目標(biāo)低仰角的目標(biāo)觀測(cè)而言，電波實(shí)際傳輸路徑較長(zhǎng)，跨越的區(qū)域范圍更大，經(jīng)歷的空間環(huán)境也更加多變復(fù)雜，因而不再滿足水平均勻的假設(shè)條件[13]。

考慮到大氣不均勻性產(chǎn)生的影響以及三維空間格點(diǎn)化大氣高度分布形式，根據(jù)射線傳播路徑獲取射線路徑上位于不同空間柵格對(duì)應(yīng)的空間環(huán)境參數(shù)值，進(jìn)而計(jì)算得到更加真實(shí)的射線傳播路徑上的折射率剖面[14]。

三維空間中的電波射線的路徑可以通過大地坐標(biāo)系（L，B，H）和站心地平直角坐標(biāo)系（E，N，U）描述，在站心地平直角坐標(biāo)系中以雷達(dá)測(cè)站位置作為站心P0建立坐標(biāo)系，電波射線軌跡以站心為起點(diǎn)向外延展到達(dá)目標(biāo)，電波射線路徑上的每個(gè)位置P可以用（EP，NP，UP）來表示，如圖2所示。

圖2 地心空間直角坐標(biāo)系和站心地平直角坐標(biāo)系Fig.2 Geocentric space and station center ground horizontal Cartesian coordinate system

NRLMSISE-00大氣模型是由美國(guó)海軍研究實(shí)驗(yàn)室設(shè)計(jì)研發(fā)的全球大氣經(jīng)驗(yàn)?zāi)Ｐ停枋隽藦牡孛娴綗釋痈叨确秶鷥?nèi)（0~1 000 km）的中性大氣密度、溫度等大氣物理性質(zhì)，是目前使用最多的大氣模型之一[15]。該模型在長(zhǎng)時(shí)間的觀測(cè)數(shù)據(jù)基礎(chǔ)上建立并不斷更新，主要數(shù)據(jù)源為火箭探測(cè)數(shù)據(jù)、衛(wèi)星遙感數(shù)據(jù)和非相干散射雷達(dá)數(shù)據(jù)等，通過采用低階球諧函數(shù)擬合大氣性質(zhì)隨經(jīng)緯度、年周期、半年周期、地方時(shí)的變化而建立的[16]。

本文使用NRLMSISE-00模型來描述和構(gòu)建雷達(dá)電波射線經(jīng)過的對(duì)流層大氣環(huán)境，通過公式(12)可以計(jì)算得到對(duì)流層的折射率剖面：

式中：P為壓強(qiáng)，hPa；T為大氣溫度，K；e為水汽壓，hPa。

NRLMSISE-00模型和IRI模型是基于大地坐標(biāo)系來計(jì)算獲取全球空間環(huán)境參數(shù)的，為了保持空間匹配，NRLMSISE-00模型也按照經(jīng)緯度1°×1°的柵格空間分辨率把雷達(dá)測(cè)站視野范圍內(nèi)的三維空間格點(diǎn)化處理，高度范圍為0~60 km，步長(zhǎng)為1 km。通過站心地平直角坐標(biāo)系和大地坐標(biāo)系的轉(zhuǎn)換，可以將電波射線路徑經(jīng)過的每個(gè)位置投影到通過模型構(gòu)建的空間柵格中，進(jìn)而獲取射線路徑上相關(guān)的空間環(huán)境參數(shù)P、T、e、ne關(guān)于高度的分布（如示意圖3所示），在此基礎(chǔ)上進(jìn)一步計(jì)算得到射線傳輸路徑上的折射率剖面。電離層折射指數(shù)n與空間環(huán)境中電子密度的關(guān)系可表示為

圖3 射線路徑上不同空間柵格中環(huán)境參數(shù)選取Fig.3 Selection of environmental parameters in different spatial grids on the ray path

式中：ne為電子密度，el/m3；f為電波頻率，Hz。

站心地平直角坐標(biāo)系轉(zhuǎn)換為空間直角坐標(biāo)系的公式為[17]

空間直角坐標(biāo)系轉(zhuǎn)換為大地坐標(biāo)系的公式為

式中：N為卯酉圈的半徑，為第一偏心率，為橢球長(zhǎng)半軸和短半軸。

2 結(jié)果分析

2.1 電子密度剖面同化模型效果比較

圖4給出了2015-03-17T9:00LT的數(shù)據(jù)同化實(shí)例，其中圖4(a)為背景模式的區(qū)域TEC分布結(jié)果，圖4(b)為進(jìn)行數(shù)據(jù)同化后的TEC分布結(jié)果。可以看出，數(shù)據(jù)同化方法實(shí)現(xiàn)了觀測(cè)數(shù)據(jù)和背景模型有效融合，數(shù)據(jù)同化后的區(qū)域TEC分布相比IRI模型直接獲取的結(jié)果展現(xiàn)了更為豐富的細(xì)節(jié)，更加符合實(shí)際。

圖4 2015-03-17T9:00LT數(shù)據(jù)同化前后區(qū)域TEC分布的對(duì)比Fig.4 Comparison of regional ionospheric TEC maps using data assimilation on 2015-03-17T9:00LT

電離層的折射率高度分布依賴于電子密度這一關(guān)鍵的環(huán)境參量，同時(shí)電子密度剖面可以有效地描述電離層垂直結(jié)構(gòu)的細(xì)節(jié)。為了驗(yàn)證三維同化模型數(shù)據(jù)同化的效果，選取處于同化區(qū)域內(nèi)的北京昌平站的電離層測(cè)高儀SAO格式觀測(cè)數(shù)據(jù)對(duì)同化前后的電子密度剖面進(jìn)行對(duì)比，SAO觀測(cè)數(shù)據(jù)來自空間環(huán)境地基綜合監(jiān)測(cè)網(wǎng)。Digisonde系列數(shù)字測(cè)高儀內(nèi)嵌的軟件，可自動(dòng)度量標(biāo)定電離圖得到電離層特征參數(shù)和電子密度剖面，生成SAO數(shù)據(jù)[18]。圖5給出了北京昌平站上空的電子密度剖面對(duì)比結(jié)果，可以直觀地看出，在大多數(shù)時(shí)刻同化后的電子密度剖面與真實(shí)的觀測(cè)值更為接近，電子密度峰值與測(cè)高儀探測(cè)結(jié)果相差更小，說明同化算法可以有效地把TEC觀測(cè)數(shù)據(jù)融合到背景場(chǎng)，改善提升背景場(chǎng)的數(shù)據(jù)精度。

圖5 2015-12-01北京昌平站上空的電子密度剖面同化結(jié)果對(duì)比Fig.5 Comparison of electron density profile assimilation results above Changping station, Beijing on 2015-12-01

2.2 電波折射誤差修正效果比較

為了驗(yàn)證本文提出的利用空間格點(diǎn)化柵格獲取電波射線路徑上折射率剖面以及結(jié)合IRI模型數(shù)據(jù)同化進(jìn)行電波折射誤差修正的效果，選取某臺(tái)工作于P波段的雷達(dá)在年積日（day of year, DOY）分別為71、72、73對(duì)高軌道衛(wèi)星目標(biāo)的三次跟蹤實(shí)測(cè)記錄作為實(shí)驗(yàn)數(shù)據(jù)，目標(biāo)在雷達(dá)跟蹤視野范圍內(nèi)仰角歷經(jīng)了從大到小的變化，俯仰范圍為4°~80°，南北跨越緯度約23°，東西跨越經(jīng)度約42°。使用GPS精密軌道數(shù)據(jù)作為測(cè)量基準(zhǔn)來確定和比較不同方法對(duì)各種距離仰角條件下目標(biāo)的折射誤差修正精度。

本文利用射線描跡法對(duì)比三種電離層折射率剖面構(gòu)建方法對(duì)修正結(jié)果的影響：1）垂直剖面法，即通過IRI模型構(gòu)建的水平均勻的電離層折射率剖面；2）射線路徑剖面法，即通過IRI模型構(gòu)建的隨傳播路徑變化的水平不均勻的電離層折射率剖面；3）數(shù)據(jù)同化法，即利用電離層同化現(xiàn)報(bào)模型構(gòu)建的電離層折射率剖面。基于以上的折射率剖面，利用射線描跡法可以計(jì)算得到折射誤差修正量和剩余殘差，雷達(dá)測(cè)距測(cè)角誤差分別為雷達(dá)探測(cè)視在距離和視在仰角與對(duì)應(yīng)目標(biāo)精密軌道值之差，剩余殘差為經(jīng)修正后的測(cè)距測(cè)角誤差。

圖6～8分別給出了三次觀測(cè)跟蹤中三種不同方法應(yīng)用于折射誤差修正中對(duì)距離誤差修正的表現(xiàn)?？梢钥闯?，隨著觀測(cè)仰角降低，三種方法計(jì)算的誤差修正量均隨之增大，在低仰角處基于模型數(shù)據(jù)同化獲取射線路徑剖面計(jì)算得到的距離折射誤差修正量與雷達(dá)測(cè)距折射誤差有更好的一致性，相應(yīng)地，修正殘差結(jié)果（子圖b）顯示利用數(shù)據(jù)同化法進(jìn)行距離折射誤差修正的殘差明顯更小，并且修正殘差伴隨觀測(cè)仰角變化的趨勢(shì)更加平緩，很大程度上消除了觀測(cè)仰角的影響，表明了數(shù)據(jù)同化方法的有效性。

圖6 DOY=71時(shí)三種方法距離折射誤差修正結(jié)果Fig.6 Correction results of distance refraction error by 3 methods when DOY=71

圖7 DOY=72時(shí)三種方法距離折射誤差修正結(jié)果Fig.7 Correction results of distance refraction error by 3 methods when DOY=72

圖8 DOY=73時(shí)三種方法距離折射誤差修正結(jié)果Fig.8 Correction results of distance refraction error by 3 methods when DOY=73

圖9～11分別給出了三次觀測(cè)跟蹤中三種不同方法應(yīng)用于折射誤差修正中對(duì)仰角誤差修正的表現(xiàn)。其低仰角情況下誤差更為顯著的結(jié)論對(duì)于仰角誤差而言仍然成立，最大仰角誤差可達(dá)0.4°。對(duì)比三種方法在仰角誤差修正中的差別可以發(fā)現(xiàn)，三種方法的修正結(jié)果依然是數(shù)據(jù)同化法更優(yōu)秀，射線路徑剖面法次之，垂直剖面法最差，但是不同于距離誤差，三種方法對(duì)于仰角誤差的修正效果差距并沒有那么明顯，或者說數(shù)據(jù)同化方法對(duì)仰角誤差的修正相對(duì)于其他兩種方法效果并不顯著，這個(gè)問題還有待更深入的研究。

圖9 DO Y=71時(shí)三種方法仰角折射誤差修正結(jié)果Fig.9 Correction results of elevation angle refraction error by 3 methods when DOY=71

圖10 DOY=72時(shí)三種方法仰角折射誤差修正結(jié)果Fig.10 Correction results of elevation angle refraction error by 3 methods when DOY=72

圖11 DOY=73時(shí)三種方法仰角折射誤差修正結(jié)果Fig.11 Correction results of elevation angle refraction error by 3 methods when DOY=73

為了量化對(duì)比三種方法的修正結(jié)果，對(duì)利用上述三種方法進(jìn)行雷達(dá)電波折射誤差修正的殘差進(jìn)行統(tǒng)計(jì)分析，結(jié)果見表1、2。其中，數(shù)據(jù)同化方法距離折射誤差修正平均殘差分別為15.43、16.97、13.64 m，相比于未進(jìn)行同化的射線路徑剖面方法的修正精度提升了41.2%、35.5%、44.6%，這表明，在三次雷達(dá)觀測(cè)修正中，通過數(shù)據(jù)同化方法進(jìn)行電波折射誤差修正能極大地提高修正精度。同時(shí)，數(shù)據(jù)同化法的修正結(jié)果均方根誤差(root mean square error,RMSE)也最小，這表明修正后的殘差更為“平坦”，也正如前文所說，對(duì)于低仰角的誤差修正得到了很大改善。而對(duì)于仰角誤差而言，如表2，其提升結(jié)果有限。

表1 不同修正方法距離折射誤差修正殘差統(tǒng)計(jì)分析Tab.1 Statistical analysis of the correction residuals of distance refraction error

表2 不同修正方法仰角折射誤差修正殘差統(tǒng)計(jì)分析Tab.2 Statistical analysis of the correction residuals of elevation refraction error

3 結(jié) 論

本文針對(duì)傳統(tǒng)折射誤差修正方法忽略大氣水平不均勻性而采用測(cè)站單點(diǎn)位置垂直方向的折射率剖面進(jìn)行電波折射誤差修正的局限性，提出了空間格點(diǎn)化大氣和電離層模型同時(shí)結(jié)合數(shù)據(jù)同化獲取電波射線路徑上電子密度剖面的方法。基于某臺(tái)P波段雷達(dá)的外場(chǎng)實(shí)測(cè)數(shù)據(jù)和精密軌道數(shù)據(jù)對(duì)三種方法應(yīng)用于折射誤差修正的精度進(jìn)行了比較分析，實(shí)驗(yàn)結(jié)果表明，通過空間格點(diǎn)化獲取射線路徑上的折射率剖面來進(jìn)行電波折射誤差修正相比于利用大氣和電離層模型直接獲取測(cè)站單點(diǎn)的折射率剖面來進(jìn)行修正有一定的改善，但仍有較大的修正剩余誤差，難以滿足工程應(yīng)用中的精度需求?；诳臻g格點(diǎn)化大氣和電離層模型結(jié)合數(shù)據(jù)同化的方法可以有效地增強(qiáng)電離層電子密度分布的準(zhǔn)確性，并顯著降低雷達(dá)測(cè)距測(cè)角的折射誤差修正殘差。研究結(jié)果驗(yàn)證了利用大氣和電離層模型數(shù)據(jù)同化方法提升電波折射誤差修正精度是一個(gè)可靠有效的途徑。