蔚娜,柳文,魯轉俠,馮靜,楊龍泉,郭文玲

?

F1層未充分發(fā)展時的電離層剖面反演

蔚娜,柳文,魯轉俠,馮靜,楊龍泉,郭文玲

中國電波傳播研究所，青島266107

摘要垂測電離圖反演對研究電離層結構及運動、電離層波傳播、空間氣象應用等具有重要意義，一直以來受到十分廣泛的重視.F1層未充分發(fā)展的回波描跡是垂測電離圖中觀測到的最多的一種情況，此時垂測描跡上表現(xiàn)為F1層到F2層的平滑過渡，而不是F1層充分發(fā)展時出現(xiàn)的尖點，然而現(xiàn)有電離層模型以及反演算法更多的是針對F1層充分發(fā)展情況，即模型剖面在F1層峰高處近似為拋物形狀，此處電子濃度剖面梯度為無窮大，這不符合F1層未充分發(fā)展時的剖面形態(tài)，即F1層剖面未到達拋物頂點處就已經達到該層電子濃度的最大值然后到達F2層，剖面梯度為一個有限的數(shù)值.因此，本文針對垂直探測中經常出現(xiàn)的F1層未充分發(fā)展情況，引入F1層模型設定臨頻這一參數(shù)，建立了F1層未充分發(fā)展時的基于移位切比雪夫多項式的電離層模型，并從電離層剖面的光滑性考慮，在所建模型的基礎上，提出了約束優(yōu)化F1層參數(shù)、F2層參數(shù)的垂測電離圖反演方法.然后，通過仿真分析，從理論上驗證了所建電離層模型和反演算法的合理性，并通過反演電子濃度剖面合成的垂測、斜測描跡與實測數(shù)據(jù)的對比，對反演算法的有效性進行了進一步的試驗驗證.

關鍵詞F1層；未充分發(fā)展；反演；電離層剖面；垂直探測

1引言

電離層是由太陽高能電磁輻射、宇宙射線和沉降粒子作用于地球的高層大氣使之電離而產生的，處于離地面以上約60～2000 km之間，是近地空間環(huán)境的一個重要組成部分.電離層由于含有足夠多的自由電子，會對無線電波的傳播產生一定的影響，其中影響最嚴重的是短波信號(Felgate和Golley，1971)，這些影響攜帶了電離層結構組成及當時狀態(tài)的特征信息，因此，人們常利用短波信號對電離層的特性進行實時監(jiān)測，以全面了解電離層的結構形態(tài)及變化、擾動等.

電離層短波垂直探測(簡稱垂測)(Hunsucker，1991)是一種歷史最悠久、至今仍然很重要的常規(guī)電離層探測手段，通過垂直向上發(fā)射不同頻率的電磁波，測量電磁波從電離層折回的時間能夠獲得反映電離層虛高與頻率關系的垂測電離圖.利用對垂測電離圖自動度量(Reinisch和Huang，1983；Galkin等，1996；Gilbert和Smith，1988；柳文等，2009；Scotto和Pezzopane，2004)獲得的垂測描跡可以反演電離層電子濃度剖面(電離層高度與等離子體頻率或電子濃度的對應關系)，它能夠提供詳細的電離層結構和變化的信息，對研究電離層結構、電離層波傳播問題、空間氣象應用等具有重要意義，因此，垂測電離圖的反演一直以來受到人們十分廣泛的重視.

目前，垂測電離圖反演方法可以歸納為以下兩種：①直接計算法，該方法直接利用實測虛高計算對應頻率的真實反射高度(簡稱真高)，主要有分片法(Titheridge，1959)、單一多項式法(Titheridge，1961)、重疊多項式法(Titheridge，1967)等，這些方法中有的是通過建立真高和虛高的聯(lián)立方程組，直接根據(jù)實測虛高求解出真高，如分片法和單一多項式法，有的是從較低的頻率到較高的頻率，通過計算逐步確定各個頻率對應的真高，如重疊多項式法，分片法也可以采用此方法，基于直接計算法的思想，Titheridge開發(fā)了垂測電離圖自動反演的POLAN程序(Titheridge，1985)；②模式法(Dyson和Bennett，1988；丁宗華等，2007；Jiang等，2014；Reinisch和Huang，1983；Scotto，2009；鄭傳青，1992)，該方法假設電離層剖面可用某種模型表征，通過尋找使基于該模型合成的垂測描跡與實測描跡在某種意義上最佳吻合的模型參數(shù)來確定電離層剖面，相較于直接計算法，模式法對于電離圖質量要求不那么苛刻，并且反演效果也不錯，應用較為普遍，國內外學者利用不同的探測數(shù)據(jù)發(fā)展了多種基于模式法的電離層參數(shù)反演方法.

基于模式法思想，Reinisch和Huang公開了一種基于移位切比雪夫多項式模型反演電離層剖面的方法(Reinisch和Huang，1983)，并開發(fā)了垂測電離圖自動反演NHPC程序應用于電離圖實時自動度量(ARTIST)系統(tǒng)中.該方法將F層建模為移位切比雪夫多項式，求解滿足計算虛高與實測虛高在最小二乘意義上最佳吻合的多項式系數(shù)，從而確定電離層剖面，后來，針對回波描跡中明顯有F1層描跡的情況，做了模型上的一些改進，即將F1層和F2層分別用不同的移位切比雪夫多項式表示(Reinisch等，1988)，但該方法更適用于F1層充分發(fā)展的情況，對于F1層未充分發(fā)展的情況，基于反演結果合成的垂測電離圖與實測電離圖在F1層臨頻附近會有一些差異，另外，未考慮層與層連接處的光滑性問題也是該方法的不足之處.

F1層未充分發(fā)展，即F1層峰高處電子濃度剖面梯度為有限數(shù)值，此時，垂測描跡上表現(xiàn)為F1層到F2層的平滑過渡，而不是F1層充分發(fā)展時出現(xiàn)的尖點.這一情況，在垂測電離圖中經常出現(xiàn)，王世凱等(2010)也給出了F1層未充分發(fā)展時基于QPS(準拋物段)模型的反演方法.本文針對F1層未充分發(fā)展情況，并從電離層剖面的光滑性考慮，在Reinisch和Huang(1983)、Reinisch等(1988)提出的電離層剖面反演方法的基礎上，建立了F1層未充分發(fā)展時電離層模型，詳細描述了反演流程，并基于仿真和實測數(shù)據(jù)對反演算法的有效性進行了驗證.

2F1層未充分發(fā)展的剖面模型

現(xiàn)有常用的電離層模型，無論是QPS模型(Dyson和Bennett，1988)、IRI(國際參考電離層)模型(Bilitza，1990)，還是Reinisch和Huang(1983)、Reinisch等(1988)提出的移位切比雪夫多項式模型，所描述的電子濃度剖面在每層的峰高處都近似為拋物形狀，此處梯度剖面為無窮大，也即較低一層剖面在拋物頂點處達到該層電子濃度的最大值，然后到達較高一層，如圖 1a給出的剖面示意圖中紅色實線所示(紅色虛線給出的是該處的拋物形狀)，當F1層充分發(fā)展時，上述模型描述F1層剖面是恰當?shù)?，此時回波描跡表現(xiàn)為在F1層臨頻處“沖上去”，形成一個尖點，如圖 1b給出的回波描跡示意圖中紅色實線所示，然而對于F1層未充分發(fā)展情況，較低一層剖面未到達拋物頂點處就已經達到該層電子濃度的最大值，然后到達較高一層，剖面梯度為一個有限的數(shù)值，如圖 1a給出的示意圖中藍色實線所示(藍色虛線給出的是該處的拋物形狀)，此時回波描跡表現(xiàn)為F1層到F2層描跡的平滑過渡，如圖 1b給出的回波描跡示意圖中藍色實線所示，對于F1層未充分發(fā)展情況，上述常用模型將不能很好地描述F1層未充分發(fā)展情況剖面的變化特征.

針對上述分析的F1層是否充分發(fā)展時剖面存在的差異，本文引入F1層模型設定臨頻這一參數(shù)，基于Reinisch和Huang(1983)、Reinisch等(1988)提出的移位切比雪夫多項式模型建立了F1層未充分發(fā)展情況的電離層模型，下面給出新建模型的詳細描述.

基于Reinisch和Huang(1983)、Reinisch等(1988)提出的移位切比雪夫多項式模型，將電離層建模為包含E層、谷層、F1層、F2層的四層模型，E層和谷層剖面用拋物模型表示，F(xiàn)1層和F2層剖面用移位切比雪夫多項式模型表示，電離層電子濃度剖面具有式(1)所示形式：

(1)

E層和谷層的連接點位于E層峰高hmE，谷層與F1層的連接點位于高度h2處，并且在高度h2處的等離子體頻率等于E層臨頻fCE，谷層包括兩個部分：與E層的連接部分和與F1層的連接部分，這兩部分的連接點位于高度h1處，F(xiàn)1層與F2層連接點位于F1層峰高hmF1處，式(1)中各符號的具體含義如下.

(1) E層：

fNE表示E層等離子體頻率；fCE表示E層臨頻；hmE表示E層峰高；ymE表示E層半厚；hbE=hmE-ymE表示E層底高；

(2) 谷層：

fNV表示谷層等離子體頻率；fCV表示谷層最小等離子體頻率；hmV表示谷層等離子體頻率為fCV時對應的電離層高度；ymV表示谷層半厚；h2=hmE+W，W定義為谷層寬度；

(3) F1層：

Ti(g)為移位切比雪夫多項式，具有式(2)所示 形式：

圖1　F1層未充分發(fā)展時剖面和回波描跡示意圖(a) 剖面示意圖；(b) 回波描跡示意圖.

(2)

(3)

fNF1表示F1層等離子體頻率；fCF1表示F1層臨頻；Ai(i=0～I+1)為移位切比雪夫多項式系數(shù)，且有

(4)

hmF1表示F1層峰高，且有

(5)

對于F1層未充分發(fā)展的情況，F(xiàn)1層模型設定臨頻不再等于fCF1，而是有一定的偏差ΔfC，也就是說，在實際的電離層中，電子濃度還沒有達到F1層模型的最大值，就已經到了F2層了；

(4) F2層：

移位切比雪夫多項式Ti(l)中的l具有式(6)所示形式：

(6)

fNF2表示F2層等離子體頻率；fCF2表示F2層臨頻；Ci(i=0～N+1)為移位切比雪夫多項式系數(shù)，且

(7)

hmF2表示F2層峰高，且

(8)

(1) 在h=hmF1處，有

(9)

(10)

式(10)為后續(xù)反演F2層系數(shù)Ci(i=0～N)的一個約束條件.

(2) 在h=h2處，有

(11)

(12)

式(12)也是后續(xù)反演F1層系數(shù)Ai(i=0～I)的一個約束條件，則由式(11)得到

(13)

(3) 在h=h1處，有

(14)

由式(14)中的上式可以得到

(15)

(16)

由式(14)中的下式可以給出h1的計算結果為

(17)

3反演各層參數(shù)

3.1各層反射回波虛高的計算

在這一節(jié)中，通過計算各層反射回波虛高，得到計算虛高與實測虛高的誤差量，從而用于實現(xiàn)后續(xù)基于最小化誤差量的各層參數(shù)反演.在計算過程中，出于簡化計算過程，同時又不會引入很大的誤差的考量，仍然采用Reinisch和Huang(1983)使用的方法，即在計算電波在E層和谷層傳播的群距離時沒有考慮地磁場的影響，在計算電波在F1層和F2層傳播的群距離時假設地磁場為一定值，即任意位置處的地磁場與垂測站上空300 km高度處的地磁場一致.

(1) E層回波虛高的計算

對于頻率小于等于fCE的電波將在E層反射，回波虛高計算公式為

(18)

其中，f為電波頻率，hr為電波反射點處高度，μ′為群折射指數(shù)，不考慮地磁場時，具有如下形式：

(19)

式中fN表示對應位置處等離子體頻率.基于建立的E層電離層模型，則式(18)可以進一步計算得到

(20)

(2) F1層回波虛高的計算

(21)

(22)

(23)

(24)

(25)

式中，

(26)

(27)

(3) F2層回波虛高的計算

對于頻率大于fCF1、小于等于fCF2的電波將在F2層反射，回波虛高計算公式為

(28)

(29)

式中，

(30)

(31)

式中，

Pi(f)=

(32)

(33)

3.2各層參數(shù)的反演

在本節(jié)中，基于對所選數(shù)據(jù)點的計算得到的虛高和實測虛高的差的平方和最小化準則，實現(xiàn)各層參數(shù)的反演，即實現(xiàn)最小二乘意義上的參數(shù)最優(yōu)化.

(1) E層參數(shù)的反演

由式(1)中上式可知，決定E層剖面的三個參數(shù)主要是fCE、hmE(或hbE)、ymE，本文采用一種區(qū)域搜索的方法實現(xiàn)E層參數(shù)的反演，雖然fCE可由垂測電離圖智能判讀軟件(柳文等，2009)自動給出，誤差小于0.2 MHz，但考慮到進一步提高反演精度，本文對fCE也在一定范圍進行搜索，確定最優(yōu)值.

(34)

使ε達到最小的那組參數(shù)即確定為E層參數(shù).

(2) 設定ΔfC時谷層參數(shù)的反演

由式(13)、式(15)和式(16)可知，得到E層參數(shù)后，谷層只要確定B和W兩個參數(shù)，谷層剖面就可以確定，因此，谷層參數(shù)反演即確定B和W兩個參數(shù).本文基于搜索、迭代方法實現(xiàn)谷層B和W兩個參數(shù)的反演.

在谷層參數(shù)的反演過程中，選擇F1層描跡中大于E層臨頻和F1層描跡最小虛高(記為h″minF1)對應的頻率之間的數(shù)據(jù)用于確定谷層參數(shù)，因為這部分描跡點對谷層參數(shù)比較敏感，假設共有K個點，其對應的工作頻率和虛高分別為fk和h″(fk).

谷層參數(shù)反演的基本步驟如下.

① 設定W=0；

② 設定B=0，I=7；

③ 基于最小二乘法計算F1層剖面系數(shù)Ai(i=0～I+1)；

④ 檢查計算的系數(shù)Ai(i=0～I)是否滿足F1層剖面單調遞增特性，(a)如果不滿足，則I=I-1，如果I<0，則執(zhí)行⑤，否則，執(zhí)行③；(b)如果滿足，對AI+1和前一次迭代記錄的值進行比較，如果二者之差小于某一較小值(例如0.5 km)，則計算K個點實測虛高和計算虛高的誤差平方和，記錄當前B、W以及計算的誤差平方和的值，否則，在沒有超過限定的最大迭代次數(shù)情況下，根據(jù)Ai(i=0～I)自動調整I，按照式(12)更新B的值，執(zhí)行③，若超過限定的最大迭代次數(shù)，則執(zhí)行⑤；

⑤ W=W+1(單位為km)，如果W小于設定的搜索范圍(例如0.7 h″minF1-hmE)，則執(zhí)行②，否則，執(zhí)行⑥；

⑥ 找出記錄的誤差平方和的最小值，該最小值對應的B、W即確定為谷層參數(shù)，如果沒有記錄下有效的B和W，則W=0，B=0.

其中，上述步驟③的具體方法為：

根據(jù)式(21)給出K個點實測虛高h″(fk)和計算虛高h′(fk)的誤差平方和：

(35)

式中，Δh′(fk)).要使ε達到最小，即求解滿足式(36)的系數(shù)Ai(i=0～I)：

(j=0～I)，

(36)

求解上述方程組即可得到系數(shù)Ai(i=0～I)，然后按照式(4)可計算AI+1.

(3) 設定ΔfC時F1層參數(shù)的反演

選取F1層描跡中F1層描跡最小虛高對應的頻率到fCF1之間的數(shù)據(jù)(或者是全部F1層描跡數(shù)據(jù))用于確定F1層參數(shù)，假設共有K個數(shù)據(jù)點，其對應的工作頻率和虛高分別為fk和h″(fk).在讀取了fCF1(由垂測電離圖智能判讀軟件自動給出)的情況下，F(xiàn)1層剖面由系數(shù)Ai(i=0～I+1)完全確定，可采用類似反演谷層參數(shù)時用到的方法來計算這些系數(shù)，這里需要注意的是，在谷層參數(shù)確定后，谷層與F1層交點處的剖面梯度也已經確定，那么，根據(jù)當前數(shù)據(jù)計算出的系數(shù)Ai(i=0～I)必須滿足式(12)，因此，F(xiàn)1層參數(shù)的反演實際上是一個約束優(yōu)化問題，即

(37)

可采用拉格朗日法求解上述問題進行F1層參數(shù)的反演，具體的步驟為：

(38)

② 對式(38)中各個自變量求偏導數(shù)，建立方程組：

(39)

③ 求解上述方程組，得到系數(shù)Ai(i=0～I)，檢查計算的系數(shù)Ai(i=0～I)是否滿足F1層剖面單調遞增特性，(a)如果不滿足，則I=I-1，如果I<0，執(zhí)行④，否則，執(zhí)行①；(b)如果滿足，執(zhí)行④；

④ 按照式(4)計算AI+1.

(4)設定ΔfC時F2層參數(shù)的反演

F2層描跡中的數(shù)據(jù)用于確定F2層參數(shù)，假設共有K個數(shù)據(jù)點，其對應的工作頻率和虛高分別為fk和h″(fk).在讀取了fCF2的情況下，F(xiàn)2層剖面由系數(shù)Ci(i=0～N+1)完全確定，在F1層參數(shù)確定后，F(xiàn)1層與F2層交點處的剖面梯度也已經確定，那么，根據(jù)當前數(shù)據(jù)計算出的系數(shù)Ci(i=0～N)必須滿足式(10)，因此，F(xiàn)2層參數(shù)的反演同樣是一個約束優(yōu)化問題，可采用類似反演F1層參數(shù)時用到的方法來計算這些系數(shù)：

(40)

(5) 谷層、F1層、F2層參數(shù)的最終確定

當F1層未充分發(fā)展時，垂測電離圖智能判讀軟件自動給出的fCF1相對于F1層模型設定臨頻的偏差ΔfC是一個未知量，理論上，ΔfC的取值可以介于0～fCF2-fCF1之間，因此，在反演谷層、F1層、F2層參數(shù)時，將ΔfC在0～fCF2-fCF1內遍歷，選取使所有數(shù)據(jù)點計算虛高與實測虛高誤差和最小的ΔfC對應的谷層、F1層、F2層參數(shù)作為最終的谷層、F1層、F2層參數(shù).

圖2給出了整個反演算法的流程框圖.

圖2　反演算法流程框圖

4理論仿真分析

本節(jié)主要是基于仿真實例對建立的電離層模型和反演算法進行仿真分析.仿真分析的思路是：假定電離層剖面參數(shù)，依此形成真實的電子濃度剖面，基于真實的電子濃度剖面仿真得到真實的回波描跡數(shù)據(jù)，然后基于真實的回波描跡數(shù)據(jù)采用不同的反演算法反演電子濃度剖面，將反演得到的剖面和真實電子濃度剖面進行比較來驗證反演算法的有效性.

表 1給出了兩個仿真實例的電離層參數(shù)，假定F1層充分發(fā)展時的臨頻為7MHz，二者的區(qū)別主要是F1層的發(fā)展程度不同.圖 3和圖 4分別給出了兩個仿真實例反演算法的結果比較.

圖3a和圖4a給出的是反演剖面和真實剖面的比較結果，其中紅色實線是真實電子濃度剖面，綠色

圖3　仿真實例1反演算法比較(a) 剖面比較； (b) 剖面誤差； (c) 虛高描跡比較； (d) 虛高誤差.

圖4　仿真實例2反演算法比較(a) 剖面比較； (b) 剖面誤差； (c) 虛高描跡比較； (d) 虛高誤差.

虛線和藍色虛線分別是基于Reinisch和Huang(1983)以及Reinisch等(1988)提出算法(后續(xù)簡記為Reinisch算法)和本文算法反演的電離層電子濃度剖面.由圖中可以看出，本文反演算法反演的電子濃度剖面和真實剖面誤差較小，圖 3b和圖 4b中紅色實線給出了本文算法反演剖面與構造真實剖面等離子體頻率之差；而Reinisch算法反演的電子濃度剖面和真實剖面誤差相對要大一些，如圖 3b和圖 4b中綠色實線所示，其中在F1層臨頻(峰高)附近，與真實剖面誤差(絕對值)約為0.2 MHz，這主要是因為，Reinisch算法給出的F1層模型剖面在F1層臨頻處接近拋物形狀，對應的剖面梯度為無窮大，而這與F1層未充分發(fā)展情況是不一致的.在仿真實例2中，Reinisch算法反演谷層剖面和真實剖面的最大誤差(絕對值)約為0.7 MHz，這主要是由于反演過程存在解的不唯一性導致.而本文反演算法，要基于不同的F1層模型設定臨頻值得到對應的谷層、F1層和F2層參數(shù)，然后從這些參數(shù)中選取使合成F1層、F2層虛高和真實回波描跡虛高誤差和最小的一組作為最終的電離層參數(shù)，降低了選擇錯誤解的概率.

表1　仿真實例電離層參數(shù)

圖3c和圖4c給出的是真實回波描跡虛高和基于反演結果合成虛高的比較結果，圖3d和圖4d給出了兩種算法反演結果合成虛高與真實回波描跡虛高之差，由圖中可見，基于本文算法反演結果合成的虛高描跡與真實回波描跡誤差較小，兩個仿真實例對應的誤差(絕對值)均值分別為0.01 km和0.04 km，而基于Reinisch算法反演結果合成的虛高描跡與真實回波描跡誤差較大，兩個仿真實例對應的誤差(絕對值)均值分別為14.22 km和3.32 km，并且基于Reinisch算法反演結果合成的虛高描跡在F1層臨頻處“沖了上去”，與真實回波描跡明顯不一致，這也是F1層充分發(fā)展和未充分發(fā)展的區(qū)別所在.

本文算法在反演F1層參數(shù)、F2層參數(shù)時均附加了該層與前一層交點處的剖面梯度保持不變這一約束條件，保證了反演剖面的連續(xù)光滑性，而Reinisch算法只保證了剖面的連續(xù)性.表2和表3分別給出了基于Reinisch算法和本文反演算法反演剖面在各層連接點處的等離子體頻率值和剖面梯度值，用以說明反演剖面的連續(xù)性和光滑性.表中“-”表示在“對應高度-0”處的值，“+”表示在“對應高度+0”處的值.由表中數(shù)值可以看出，在連接點以上和以下趨近于連接點的位置處，等離子體頻率都相等，證明了兩種反演算法得到的剖面都是連續(xù)的；Reinisch算法反演剖面在h2和hmF1處剖面梯度出現(xiàn)了跳變，即左右導數(shù)不相等，如表 2中黑色加粗數(shù)據(jù)所示，說明反演的剖面不光滑，而本文算法反演剖面在各個連接點處的左右導數(shù)都相等(稍許誤差是數(shù)值計算誤差所致)，說明了反演的剖面在各個連接點處保持了光滑性.

表2　Reinisch算法反演剖面各層連接點處計算結果

表3　本文算法反演剖面各層連接點處計算結果

5實測數(shù)據(jù)處理分析

基于中國電波傳播研究所記錄的大量垂測數(shù)據(jù)對反演算法進行了驗證，本節(jié)選取了7組數(shù)據(jù)對算法的有效性加以說明.選取的數(shù)據(jù)均為典型三層(F1層未充分發(fā)展)電離層情況.

圖5給出了典型三層(F1層未充分發(fā)展)電離層情況實測電離圖(左列圖)以及垂測反演結果(右列圖)，在右列圖中，黑色“*”描跡給出的是垂測電離圖智能判讀軟件對實測電離圖的判讀結果，藍色虛線和玫紅色虛線分別是基于Reinisch算法和本文算法反演的電離層電子濃度剖面，綠色實線和紅色實線分別是基于上述兩種算法反演的電子濃度剖面合成的垂測虛高描跡.

Reinisch算法沒有考慮F1層未充分發(fā)展情況，基于反演算法得到的電子濃度剖面合成的垂測描跡在F1層臨頻附近均與實測描跡有很大差別，如圖5中黑色橢圓框入的綠色實線所示，而本文算法考慮了F1層未充分發(fā)展情況，基于反演算法反演的剖面合成的回波描跡與實測描跡在F1層臨頻附近都吻合得較好，如圖5中紅色實線所示.

圖5　典型三層電離層情況垂測電離圖和垂測反演結果(a)—(g)分別為2014年08月03日12時06分、2014年12月11日15時00分、2015年08月23日13時48分、2015年08月23日14時03分、2015年08月23日14時33分、2015年08月23日16時33分和2015年08月24日12時18分的結果，左側為實測電離圖、右側為反演剖面.

圖5(c1)—5(g1)中的實測電離圖顯示，此時F1層或者F2層均存在一定程度的擾動，兩種算法對于頻率間變化比較緩慢的擾動均具有一定的魯棒性，本文算法性能更好一些，如圖5(d2)所示.但是對于明顯存在附加電離層的擾動，如圖 5(c1)和5(f1)所示，出現(xiàn)了F0.5層和F3層，此時基于兩種反演算法反演剖面合成的描跡都沒有表現(xiàn)出實測描跡的擾動特征，因此，對于這類情況，還需要在電離層模型上增加擾動項模型，進行進一步的改進.

由于Es層遮蔽或者吸收的影響，F(xiàn)1層低端常常會出現(xiàn)描跡嚴重缺失，如圖5(a2)和5(d2)所示，此種情況下，兩種算法均能自動反演出電子濃度剖面，實現(xiàn)對提取的實測描跡在最小二乘意義上的最佳擬合，但對應區(qū)域反演剖面可能存在一定的差別.而由于實測數(shù)據(jù)的缺失以及沒有真實剖面的比對，無法說明哪種反演算法對這一區(qū)域的反演結果更加準確.因此，對于描跡缺失問題，應該考慮增加具有實際物理意義的缺失數(shù)據(jù)補償算法，這也是進一步提高反演精度的必要條件.

在中緯度地區(qū)、中午前后，F(xiàn)2層臨頻東西向的梯度一般較小(Davies和Rush，1985).因此，對于中緯度地區(qū)東西向1000 km左右的鏈路，中午前后，電子濃度剖面均勻分布的假設近似成立.在此基礎上，選取了符合上述條件的一條斜測鏈路，基于本文算法反演的電子濃度剖面，在電子濃度剖面為均勻分布的假設下，利用自研的基于三維數(shù)字射線追蹤(Jones和Stephenson，1975；柳文等，2008)和離散電子濃度網格的斜測描跡合成軟件，合成對應該斜測鏈路的斜測描跡與實測電離圖進行對比，對反演電子濃度剖面的準確性做進一步的驗證.

圖 6給出了對比結果，其中，玫紅色點跡表示基于本文提出反演算法反演電子濃度剖面合成的斜測描跡.由圖中可以看出，合成描跡E層低角的群距離和實測數(shù)據(jù)吻合得較好，誤差在斜向探測一個距離分辨率之內(7.5 km)；圖 6a中實測F1層描跡大部分缺失，但根據(jù)描跡的趨勢可以判斷合成的描跡與已有描跡是一致的.圖6b中沒有觀察到明顯的F1層描跡，合成F1層描跡群距離與對應頻率上的實測描跡有一些差別，偏差最大達30 km，這主要是由于電子濃度均勻分布的假設和實際情況有一定的出入所致；圖 6a中合成F2層描跡高低角波都與實測結果一致，圖 6b中隨著頻率的增大，合成F2層描跡群距離與實測值逐漸接近，但MUF比實測值略低，這些差別也主要是由于電子濃度分布的不均勻性造成的.

綜上，由實測數(shù)據(jù)和基于反演結果合成的描跡的對比可知，對于F1層未充分發(fā)展情況，本文算法能夠根據(jù)提取的垂測描跡較為準確地給出上空電離層狀態(tài).

圖6　基于垂測反演結果合成的斜測描跡與實測電離圖對比情況(a) 2014年08月03日12時06分； (b) 2014年12月11日15時00分.

6結論

本文在Reinisch和Huang(1983)、Reinisch等(1988)提出的一種基于移位切比雪夫多項式模型反演電離層剖面的方法的基礎上，針對F1層未充分發(fā)展情況，引入F1層模型設定臨頻這一參數(shù)，建立了F1層未充分發(fā)展時的電離層模型，提出了約束優(yōu)化F1層參數(shù)、F2層參數(shù)的垂測電離圖反演方法，該方法重新歸納如下：利用E層描跡數(shù)據(jù)，通過局部尋優(yōu)實現(xiàn)E層參數(shù)的反演；給出一個F1層模型設定臨頻值，利用F1層較低頻率區(qū)域的數(shù)據(jù)，基于搜索、迭代的方法實現(xiàn)谷層參數(shù)的反演，得到谷參數(shù)后，選取F1層較高頻率區(qū)域回波描跡數(shù)據(jù)(或者全部F1層回波描跡數(shù)據(jù))，在保證剖面連續(xù)光滑的約束條件下，計算F1層剖面多項式系數(shù).同樣，在保證剖面連續(xù)光滑的約束條件下，選取F2層回波描跡數(shù)據(jù)，計算F2層剖面多項式系數(shù)；然后，改變F1層模型設定臨頻值，重復上述過程，得到每個F1層模型設定臨頻值對應的各層剖面參數(shù)；最后，基于F1層、F2層描跡點的計算虛高和實測虛高誤差和最小準則，選取對應的F1層模型設定臨頻下得到的剖面參數(shù)，最終確定電離層剖面.

本文提出的反演算法能夠對F1層未充分發(fā)展的垂測電離圖進行有效反演，并且反演的電子濃度剖面具有連續(xù)光滑特性，反演有效性利用仿真和實測數(shù)據(jù)進行了驗證.需要說明的是，本文提出反演算法只是針對F1層未充分發(fā)展情況，因為在反演F1層和F2層參數(shù)時使用了約束優(yōu)化方法，因此，對于F1層充分發(fā)展情況，式(40)中的約束條件為a=∞，將無法結合約束條件計算出F2層參數(shù)，因此，需要反演算法增加預處理單元，判斷此時F1層的發(fā)展情況.

另外，在本文反演模型中加入擾動項，將有助于對存在擾動情況的電離圖反演精度的進一步提高；對于回波描跡缺失問題，考慮增加具有實際物理意義的缺失數(shù)據(jù)補償算法；以及考慮使用O波、X波聯(lián)合反演進一步提高反演算法的穩(wěn)定性；等等.以上這些都是使算法更適合實際工程應用需要改進的地方，這也是今后需要進一步努力的方向.

References

Bilitza D.1990.The international reference ionosphere.NSSDC/WDC-A-R&S Report 90-22.National Space Science Data Center.Greenbelt，Maryland.Dyson P L，Bennett J A.1988.A model of the vertical distribution of the electron concentration in the ionosphere and its application to oblique propagation studies.Journal of Atmospheric and Terrestrial Physics，50(3)：251-262.

Davies K，Rush C M.1985.High-frequency ray paths in ionospheric layers with horizontal gradients.Radio Science，20(1)：95-110.

Ding Z H，Ning B Q，Wan W X.2007.Real-time automatic scaling and analysis of ionospheric ionogram parameters.Chinese Journal of Geophysics，50(4)：969-978.

Felgate D G，Golley M G.1971.Ionospheric irregularities and movements observed with a large aerial array.Journal of Atmospheric and Terrestrial Physics，33(9)：1353-1369.

Galkin I A，Reinisch B W，Ososkov G A，et al.,1996.Feedback neural networks for ARTIST ionogram proceedings.Radio Science，31(5)：1119-1128.

Gilbert J D，Smith R W.1988.A comparison between the automatic ionogram scaling system ARTIST and the standard manual method.Radio Science，23(6)：968-974.

Huang X Q，Reinisch B W.1982.Automatic calculation of electron density profiles from digital ionograms 2.True height inversion of topside ionograms with the profile-fitting method.Radio Science，17(4)：837-844.

Hunsucker R D.1991.Radio Techniques for Probing the Terrestrial Ionosphere.Springer-Verlag.

Jiang C H，Yang G B，Zhao Z Y，et al.,2014.A method for the automatic calculation of electron density profiles from vertical incidence ionograms.Journal of Atmospheric and Solar-Terrestrial Physics，76：20-29.Jones R M，Stephenson J J.1975.A versatile three-dimensional ray tracing computer program for radio waves in the ionosphere.OT Report 75-76.US Department of Commerce，Office of Telecommunications，US Government Printing Office.Washington，USA,185 pp.Liu W，Jiao P N，Wang S K，et al.,2008.Short wave ray tracing in the ionosphere and its application.Chinese Journal of Radio Science (in Chinese)，23(1)：41-48.

Liu W，Kong Q Y，Chen Y，et al.,2009.Method on ionogram autoscaling based on IRI model.Chinese Journal of Radio Science (in Chinese)，24(2)：218-223.

Reinisch B W，Huang X Q.1983.Automatic calculation of electron density profiles from digital ionograms 3.Processing of bottomside ionograms.Radio Science，18(3)：477-492.

Reinisch B W，Gamache R R，Huang X Q，et al.,1988.Real time electron density profiles from ionograms.Advances in Space Research，8(4)：63-72.

Scotto C，Pezzopane M.2004.Software for the automatic scaling of critical frequency foF2 and MUF(3000)F2 from ionograms applied at the ionospheric observatory of gibilmanna.Annales Geophysicae，47(6)：1783-1790.

Scotto C.2009.Electron density profile calculation technique for autoscala ionogram analysis.Advances in Space Research，44(6)：756-766.

Titheridge J E.1959.The calculation of real and virtual heights of reflection in the ionosphere.Journal of Atmospheric and Terrestrial Physics，17：96-109.Titheridge J E.1961.A new method for the analysis of ionospheric records.Journal of Atmospheric and Terrestrial Physics，21：1-12.

Titheridge J E.1967.The overlapping-polynomial analysis of ionograms.Radio Science，2(10)：1169-1175.

Titheridge J E.1985.Ionogram analysis with the generalized program POLAN.Report UAG-93.World Data Center A for Solar-Terrestrial Physics.Boulder，Colorado.

Wang S K，Liu W，Li S L，et al.,2010.Inversion method of F1 layer on vertical sounding ionogram.Chinese Journal of Radio Science (in Chinese)，25(1)：172-177.

Zheng C Q.1992.The inversion of the vertical sounding ionogram.Chinese Journal of Radio Science (in Chinese)，7(4)：64-74.

附中文參考文獻

丁宗華，寧百齊，萬衛(wèi)星.2007.電離層頻高圖參數(shù)的實時自動度量與分析.地球物理學報，50(4)：969-978.

柳文，孔慶顏，陳躍等.2009.基于IRI模型的垂測電離圖自動判讀算法.電波科學學報，24(2)：218-223.

柳文，焦培南，王世凱等.2008.電離層短波三維射線追蹤及其應用研究.電波科學學報，23(1)：41-48.

王世凱，柳文，李書苓.2010.垂測電離圖F1層反演方法研究.電波科學學報，25(1)：172-177.鄭傳青.1992.垂直探測頻高圖的反演.電波科學學報，7(4)：64-74.

(本文編輯胡素芳)

基金項目中國電科技術創(chuàng)新基金項目(JJ-QN-2013-28)，國防技術基礎科研項目(JSHS2014210A002)，國家自然科學基金(61331012)共同資助.

作者簡介蔚娜，女，1981年生，高級工程師，博士，2008年畢業(yè)于哈爾濱工程大學.目前主要從事雷達信號處理和電波傳播研究工作. E-mail: lenghanbingxue@sina.com

doi:10.6038/cjg20160302 中圖分類號P352

收稿日期2015-10-09，2015-12-23收修定稿

The electron density profile inversion for incompletely developed case of F1layer

WEI Na,LIU Wen,LU Zhuan-Xia,FENG Jing,YANG Long-Quan,GUO Wen-Ling

ChinaResearchInstituteofRadiowavePropagation,Qingdao266107,China

AbstractThe electron density profile inversion from vertical incidence ionograms,which is essential for research in ionospheric structures and movements,wave propagation and space weather applications,has gathered very wide attention.The echo trace of incompletely developed F1 layer is common in vertical incidence ionograms,and it is usually expressed as smooth transition from F1 layer to F2 layer,not a cusp appeared at the critical frequency of F1 layer.However,the existing ionospheric models and inversion algorithms are generally intended for the completely developed F1 layer with the assumptions of a parabolic profile and an infinite slope at the peak of F1 layer,which are not suitable for the profile of incompletely developed F1 layer which achieves the maximum electron density of F1 layer and enters F2 layer at the peak of F1 layer and has a finite slope.Consequently,a F1 layer electron density profile model based on the shifted Chebyshev polynomial for incompletely developed case of F1 layer is introduced with a parameter named model setting critical frenquency.Taking into account the profile smoothness,an electron density profile inversion algorithm with constrained optimization F1 and F2 layer parameters based on the model mentioned above is proposed.The validity of the model and the inversion algorithm is analyzed through the simulation,and the effectivity of the proposed algorithm is further verified by the comparison between the synthesized vertical sounding & oblique sounding traces and the measured data.

KeywordsF1 layer; Incompletely developed; Inversion; Ionosphere profile; Vertical sounding

蔚娜,柳文,魯轉俠等.2016.F1層未充分發(fā)展時的電離層剖面反演.地球物理學報，59(3):778-790,doi:10.6038/cjg20160302.

Wei N,Liu W,Lu Z X,et al.2016.The electron density profile inversion for incompletely developed case of F1layer.Chinese J.Geophys.(in Chinese),59(3):778-790,doi:10.6038/cjg20160302.