基于不同散射機制特征的海雜波時變多普勒譜模型?

張金鵬 張玉石 李清亮 吳家驥

1)(中國電波傳播研究所,電波環境特性及模化技術重點實驗室,青島 266107)

2)(西安電子科技大學電子工程學院,西安 710071)

1 引 言

雷達海雜波多普勒譜本質上是一種用于表征海雜波在各頻帶上能量分布的功率譜,由于海面風力及重力、張力的作用通常導致海面散射體處于運動狀態,而這種雷達與散射體之間的相對運動又會引起雷達相干脈沖間的相位變化,使其產生多普勒頻移,因此也稱為多普勒譜.海雜波的多普勒譜建模研究對采用多普勒處理技術的雷達進行有效的海雜波抑制和目標檢測具有重要的意義[1?3].

海雜波多普勒譜的建模需基于相應的物理機理和特性研究.早在20世紀60年代,Pidgeon[4]就利用C波段小擦地角海雜波數據分析了多普勒譜的特性,認為海面散射主要由海面小尺度毛細波產生的Bragg散射起主導作用,并由此解釋海雜波的多普勒譜特征.到90年代,多個獨立的海雜波觀測實驗發現,單純通過Bragg散射機制并不能完全解釋海雜波多普勒譜的特性.因此,Lee等[5]和Walker[6]開展了專門的海雜波譜特性研究實驗,通過結合同步錄取的海浪視頻信息,研究了海面演化過程中各海面結構對海雜波譜的貢獻,分析了海面不同散射機制與海雜波譜特性的關聯.近年來,以Rosenberg等[7?9]為代表的澳大利亞國防科技部研究人員對海雜波的幅度、相關和譜等特性進行了研究,其中包括不同雷達參數和不同海況下的海雜波多普勒譜特性分析.Corretja等[10]對小擦地角下的海雜波多普勒譜特性進行了統計分析,并對比了經驗模型與實驗數據的符合程度.

根據海雜波多普勒譜產生的物理機理,在海雜波觀測時間不同的情況下,其多普勒譜表現出不同的形狀與非平穩特性[11,12].因此,海雜波多普勒譜的建模主要從兩方面開展研究.一是針對平均多普勒譜形狀的建模,即較長時間(通常大于重力波周期,秒級)海雜波譜的平均特性建模.Lee等[5]首先采用多組實測數據分析的方式,研究了平均多普勒譜形狀的特點和規律,提出了譜線形狀分解為三種具有物理意義基函數的建模方法,每種基函數用于表征不同的散射機制.Walker[6,13]在上述建模方法的基礎上,利用造浪池海雜波數據深入分析了海浪從產生到破碎全過程的多普勒譜變化規律,建立了一種簡化的三分量海雜波平均多普勒譜模型.該模型利用三個高斯函數分別表征Bragg,白冠和破碎三種散射機制的譜分量,在大部分情況下可以成功對實測海雜波多普勒譜形狀進行建模.

海雜波多普勒譜建模的第二個方面是短時動態多普勒譜的建模,即對較短觀測時間內(通常小于重力波周期,大于白冠和破碎散射的去相關時間)譜形狀及其變化特性的建模[14].Miller和Dawber[15]在短時譜建模中做了開創性工作,他們分析了多波段雷達的海雜波數據,建立了高斯函數形式的譜結構,并由兩個服從伽馬分布的變量對其進行調制的短時多普勒譜模型.考慮到實際的短時譜結構通常是非高斯的,Ward等[16]提出了采用復合K分布中的伽馬函數對短時譜強度進行調制,利用兩個高斯函數對譜形狀進行描述的短時多普勒譜建模方法.該方法在后來的短時譜建模中得到了成功應用,但由于該模型中假定海面的Bragg散射譜分量具有零多普勒頻移,且白冠散射和破碎波散射共享相同的多普勒頻移和展寬,因此該模型在部分實測海雜波短時譜建模中精度不高.

本文考慮上述傳統模型的缺點,提出了一種新的動態譜模型,稱為時變多普勒譜模型.該模型能夠同時適用于平均多普勒譜與短時多普勒譜的建模,其性能通過黃海海域實測的P波段和S波段岸基海雜波數據進行了驗證.

2 海雜波多普勒譜建模方法

海雜波的多普勒譜特性與海表面的運動和擾動狀態密切相關,雷達與海面散射體之間的相對運動使得電磁波產生多普勒頻移,而海面散射體運動的隨機性使得多普勒譜具有一定的展寬.因此,海雜波多普勒譜的建模方法需要基于多普勒譜的物理機理分析.

在假設雷達分辨單元內海浪相對雷達的運動速度基本不變的情況下,海雜波的多普勒頻移可以表示為

多普勒展寬可以表示為

式中,λ為電磁波波長;V為海浪相對雷達的運動速度;φ為海浪速度方向與雷達視向的夾角;θ1,θ2分別為對應分辨單元前后兩邊緣的擦地角.

對實際的動態海面而言,它是在近似周期的風浪和涌浪上疊加著小尺度的波紋、泡沫和浪花,每種散射體的運動速度和方向不盡相同.對海雜波多普勒譜而言,其譜特性直接與動態海面的多種散射機制有關,相應地表現為多普勒譜包含多種譜分量[17,18].通常情況下,譜形狀的建模需根據具體情形考慮如下三種散射機制譜分量.

1)Bragg散射譜分量

海面Bragg散射指的是當海浪波長在入射電磁波方向的投影等于電磁波半波長的整數倍時產生的諧振(相干)散射.海浪波長與電磁波波長的關系為

式中,Λ為諧振海浪波長,λ為電磁波波長,θ為擦地角.

根據(3)式給出的產生諧振散射的條件可知,Bragg散射通常發生在兩種尺度的海浪結構上.當雷達頻率較低時,布拉格散射主要來自于海面的重力波結構,而當雷達頻率較高時,Bragg散射主要來自于海面的張力波結構.對存在兩種海浪尺度的復合海面模型而言,Bragg散射回波的多普勒頻移為

式中,vB代表Bragg諧振散射波速度;vD代表重力波的漂移和軌道速度,用于對張力波進行調制.

對于P及以下波段而言,由于Bragg散射主要來自海面重力波結構,而重力波可沿雷達視向接近或遠離觀測點運動,因此海雜波多普勒譜中經常出現兩個關于零頻對稱的譜峰[19].

2)破碎散射譜分量

海浪的破碎散射指的是強風驅使下的海浪破碎波對電磁波的準鏡面反射,該散射機理使得在較高海況下,L及以上波段的小擦地角海雜波多普勒譜中經常含有獨特的破碎散射譜分量[20].

由于破碎散射來源于準鏡面反射,且海浪破碎波的運動速度明顯高于Bragg散射波,因此海雜波的破碎散射譜分量具有以下特點[16]:

3)白冠散射譜分量

白冠散射指的是雷達波照射在海浪破碎后形成的泡沫浪花上的后向散射.海浪的破碎散射和白冠散射是順序發生的,白冠散射之前必然已經產生破碎散射,而破碎散射之后卻不總有白冠散射發生.由于海面白冠層呈現了體散射特點,極化不敏感性導致了HH和VV極化下該散射幾乎具有相同的強度.白冠散射持續時間長(秒級)但相關時間極短(幾個毫秒).體現在海雜波的多普勒譜上,白冠散射的多普勒譜頻移取決于重力波的相速度,遠大于Bragg散射的譜頻移,而由于白冠散射來源于近隨機的泡沫浪花的后向體散射,因此白冠散射的多普勒譜具有類似于噪聲的寬多普勒譜特點.通常情況下,一般的實測海雜波多普勒譜可以結合上述的三種譜分量對其進行描述和建模,但由于實際的動態海面還可能受一些不確定因素的擾動,因此特殊的海雜波多普勒譜需要根據實際情況進行物理機理分析與建模.

3 海雜波多普勒譜模型

3.1 傳統多普勒譜模型

在長時平均多普勒譜和短時動態多普勒譜建模中,實用性較高且被廣泛接受的模型分別為引言中介紹的Walker模型[13]和Ward模型[16].

Walker海雜波多普勒譜模型是一種簡化的Lee模型,適用于對觀測時間大于重力波周期(通常為幾秒)的海雜波多普勒譜進行建模,屬于一種長時平均譜模型.該模型假定三種不同散射機制的譜分量都使用高斯譜線函數來描述,對于HH極化,譜模型形式為

對于VV極化,譜模型形式為

式中,f為頻率;fB和fP分別為對應Bragg諧振波速度和重力波相速度的多普勒頻率;B,W和S分別代表Bragg、白冠和破碎散射分量的強度,wB,wW和wS分別代表三種譜分量的多普勒展寬.

Ward模型適用于對觀測時間小于重力波周期但大于海浪破碎散射和白冠散射去相關時間的海雜波多普勒譜進行建模,屬于一種短時多普勒譜模型.該模型假定多普勒譜形狀被一個符合Γ分布的隨機功率變量所調制,在某一短觀測時間內,該變量為常數,多普勒譜形狀分解成諧振波高斯譜分量和非諧振波高斯譜分量兩部分.其模型形式如下:

其中

(7)式分子中第一項對應于Bargg散射(諧振波散射)的多普勒譜分量,多普勒頻移為零且多普勒展寬為wBragg;第二項對應于海浪破碎波散射和白冠散射(非諧振波散射)的多普勒譜分量,具有多普勒頻移fd和多普勒展寬wBreaking;x為隨觀測時間區間變化且服從Γ分布的功率調制因子,A為譜幅度.

3.2 時變多普勒譜模型

上述兩種傳統的多普勒譜模型在實際的譜建模中都存在一定的局限性且適用范圍互不重疊,本文在此基礎上提出了一種時變多普勒譜模型,克服了傳統模型的譜形狀建模的局限性,且從物理意義上其適用范圍包含長時平均譜和短時動態譜建模.該模型的建立主要考慮以下三個方面.

1)Ward短時譜模型認為Bragg散射具有零多普勒頻移.然而,從散射機理上考慮,由于海浪存在大尺度重力波的運動,實際的Bragg散射是存在多普勒頻移的,且當雷達波長較長、趨近于重力波長時,多普勒頻移更為明顯.因此,本文提出的時變多普勒譜模型中引入了Bragg散射的多普勒頻移.

2)Ward模型假定海雜波破碎散射和白冠散射機制共同產生一個高斯譜分量,具有相同的多普勒頻移和展寬.然而,根據多普勒譜物理機理可知,兩種散射的產生機理和相關時間并不相同,譜分量具有不同的頻移和展寬特性.因此,在本文提出的時變多普勒譜模型中,認為海雜波多普勒譜由來自Bragg,白冠和破碎三種散射機制的三個譜分量組成,且頻移和展寬不同.

3)Walker譜模型中白冠散射和破碎散射共享相同的多普勒頻移,這與實際情況是不相符的.在實際的動態海面中,海面白冠通常是附著在重力波上以其相速度運動的,而海面的破碎波(這里指未完全破碎之前的卷浪結構)通常在重力波相速度的基礎上還附加了一個由瞬時風和垂直重力加速度引起的附加速度,因此本文在提出的時變多普勒譜模型中,將破碎散射譜分量的譜頻移認為由重力波相速度引起的頻移和附加速度引起的頻移兩部分組成.

鑒于以上三點考慮,并將三個譜分量的譜強度假設為受海雜波觀測時間(譜估計時間)區間影響的隨機變量,提出海雜波時變多普勒譜模型如下:

式中

這里,IB,IW和IS分別為當雷達觀測時間為Δt情況下隨時間動態變化的Bragg,白冠和破碎散射譜分量的強度;ΨB,ΨW和ΨS為三種散射譜分量的譜線形狀基函數;fB和fG分別為對應Bragg諧振波和重力波相速度的多普勒頻率;Δfs表示由瞬時風和垂直重力加速度引起的在重力波相速度基礎上的附加速度頻移量,當重力波相速度頻移量fG為正時,Δfs前的符號取正,反之相反.該符號在物理意義上表征附加速度是沿著重力波相速度方向的.

值得注意的是,該時變多普勒譜模型中,三個譜分量的強度之和IB+IW+IS符合K分布中的調制分量Γ分布.在對(9)—(13)式給出的時變多普勒譜進行參數優化時,需附加約束條件:

此條件可以保證白冠散射和破碎散射的多普勒頻移(主要來自于重力波相速度)大于Bragg散射的譜頻移,這與海面散射機理是相符的.

4 模型性能的實驗數據驗證

本文基于黃海海域實測的岸基P波段和S波段雷達海雜波數據,利用提出的時變多普勒譜模型對實測數據的短時譜和平均譜進行了建模,并分別與傳統的Ward短時譜模型和Walker平均譜模型建模結果進行了對比,驗證了模型性能.

4.1 短時譜建模性能驗證

圖1給出了一組岸基P波段雷達海雜波數據在40 s時間內的短時多普勒譜,通過非參數化的Welch法估計得到.測量雷達為全相參脈沖體制,脈沖重復頻率為1 kHz.相應的測量條件為:2級海況、順浪向,擦地角為5.2°.譜估計過程中,先后將海雜波數據分成400個和40個觀測時間區間,每個區間的長度分別為100 ms和1000 ms,相應的短時譜估計中每個觀測時間區間的脈沖數分別為100和1000.

從圖1可以看出,在觀測時間為100 ms的情況下,多普勒譜形狀中關于零頻對稱的兩個Bragg峰非常明顯,這是由于觀測時間較短,可以獲得海雜波動態譜精細結構的原因.然而,當觀測時間較長(1000 ms)時,兩個Bragg峰已經分辨不開,表現為一個零頻附近的主峰.若觀測時間繼續變長,當大于重力波周期時,海雜波譜將會隨時間基本不變,變為平均多普勒譜.

圖1 不同觀測時間情況下的P波段海雜波短時多普勒譜(a)100 ms;(b)1000 msFig.1.Short-time doppler spectra at P band for the measured time of(a)100 ms and(b)1000 ms.

分別采用傳統的Ward短時多普勒譜模型和本文提出的時變多普勒譜模型對圖1給出的實測海雜波短時譜進行建模,結果如圖2和圖3所示.譜建模中我們引入粒子群(PSO)算法[21]用于實測譜形狀的優化,PSO初始參數設置如表1所列,各模型參數的搜索范圍根據實測譜確定.將圖2和圖3分別與圖1進行對比,可以得到如下結論.

1)兩種短時譜模型都能夠實現對隨時間動態變化的海雜波短時多普勒譜建模,建模結果基本能夠描述短時譜的時間非平穩性.

2)時變譜模型建模效果優于Ward短時譜模型.從譜的形狀上來看,時變譜模型能夠較好地描述短時譜的各個譜分量,例如圖3(a)中時變譜的建模結果能夠顯示出多普勒譜關于零頻對稱的兩個Bragg峰,與圖1(a)的實測譜一致.而Ward短時譜的建模結果中,由于假設了Bragg散射譜分量的頻移為零,兩個Bragg峰被建模為一個移動到零頻上的峰,譜形狀的精細結構消失.對譜形狀的建模結果差異是由譜模型的函數形式決定的,Ward短時譜模型由兩個高斯函數組成,因此只能描述含兩個譜峰的多普勒譜,而時變譜模型在Ward模型基礎上引入了海浪破碎波的譜分量,因此可以描述含三個譜峰的多普勒譜.

圖2 Ward短時譜模型對P波段海雜波短時多普勒譜的建模效果(a)觀測時間100 ms;(b)觀測時間1000 msFig.2.The short-time doppler spectra at P band modeled by the Ward model for the measured time of(a)100 ms and(b)1000 ms.

圖3 時變多普勒譜模型對P波段海雜波短時多普勒譜的建模效果(a)觀測時間100 ms;(b)觀測時間1000 msFig.3.The short-time doppler spectra at P band modeled by the introduced time-varying model for the measured time of(a)100 ms and(b)1000 ms.

表1定量地給出了兩種譜模型在上述P波段短時譜建模中的時間與誤差.對比建模誤差可以看出,時變譜模型的建模誤差總體上小于Ward模型.當短時譜估計中觀測時間較長(1000 ms)時,時變譜模型相對Ward模型建模精度提高的更大,這是因為當觀測時間較長時,海雜波多普勒譜的精細結構不顯著,更易于建模的原因.如圖1(b)中實測短時譜的精細結構不如圖1(a)明顯,譜建模難度降低.在建模優化時間方面,時變譜模型情況下的單譜(即一個觀測時長的譜)的優化時間約為Ward模型的兩倍,這是由于前者的模型參數較多,PSO優化時種群規模為后者情況的兩倍.但單譜優化時間小于1 s,保證了實際應用中的可行性.

表1 P波段短時譜建模時間與建模誤差Table 1.The times and RMS errors in modeling the short-time doppler spectra at P band.

圖4給出了一組岸基S波段雷達海雜波數據在100 s時間內的短時多普勒譜,通過非參數化的Welch法估計得到.海雜波數據的測量模式為岸基波束駐留模式,譜估計過程中根據雷達的相干脈沖數限制,每個觀測時間區間對應的脈沖數為100,時間長度約為230 ms.海雜波數據的測量條件為:2級海況、順浪向,擦地角為5.5°.

從圖4可以看出,該組數據多普勒譜的形狀隨時間起伏變化,主要體現為負頻移上的多普勒主峰與正頻移上較弱的多普勒次峰,主峰的頻移與展寬隨時間變化較慢,而次峰隨時間變化很快,該現象恰恰體現了實際海面速度分量的不斷變化,表現出譜的短時動態特性.

同時采用兩種譜模型對圖4給出的短時譜進行建模,優化算法配置同P波段一致,結果如圖5和表2所列.對比圖4和圖5可以看出:

圖5 岸基S波段海雜波短時多普勒譜建模結果(a)Ward模型;(b)時變譜模型Fig.5.The short-time doppler spectra at S band modeled by(a)the Ward model and(b)the introduced time-varying model.

表2 S波段短時譜建模時間與建模誤差Table 2.The times and RMS errors in modeling the short-time doppler spectra at S band.

1)時變譜模型能夠實現對隨時間動態變化的海雜波短時多普勒譜建模,建模結果基本能夠描述短時譜的時間非平穩性;

2)時變譜模型建模精度和穩定性優于Ward譜模型,時變譜模型對于負頻移上的多普勒主峰以及正頻移上的多普勒次峰刻畫基本準確,而Ward模型建模結果基本體現不出多普勒次峰,造成多普勒譜中的譜分量丟失.

表2定量地列出了兩種譜模型在上述S波段短時譜建模中的時間與誤差.可以看出,本文時變譜模型的建模誤差明顯小于Ward模型,且在單譜優化時間上,兩種模型情況下都小于1 s,符合實際應用可行性要求.

4.2 平均譜建模性能驗證

圖6給出了兩種譜模型對岸基P波段海雜波數據平均多普勒譜的建模效果.兩組典型數據對應的海況分別為3級和4級.測量雷達參數與圖1所示數據的雷達參數相同,擦地角為3.6°.實測平均譜估計所使用的雜波序列長度為60000,對應的雷達觀測時間為1 min,PSO算法的參數設置與表1相同.從圖6(a)可以看出,在3級海況下,P波段海雜波平均譜具有兩個明顯的譜峰,且在?50 Hz左右具有一個較弱的類似譜峰的突起.本文時變譜模型和Walker譜模型都能夠實現對前兩個明顯譜峰的準確描述,然而對于后一個較弱的譜突起結構,時變譜模型的建模結果體現出了該細微結構,而Walker譜模型忽略了該結構.從圖6(b)可以看出,在4級海況下,P波段海雜波平均譜表現出三個明顯的譜峰,說明了較高海況下海浪的速度分量更為復雜.在該情況下,時變譜模型可以實現整個譜形狀(三個譜峰)的準確建模,而Walker模型忽略了正頻移上的譜峰,且對中心譜峰強度的描述也存在明顯偏差.

圖7給出了兩組典型岸基S波段海雜波數據的平均多普勒譜.實測譜對應的海況分別為3級和4級,擦地角為4.1°.譜估計所使用的雜波序列長度為60000,對應的雷達觀測時間為1 min.從圖7可以看出,兩種海況下S波段的平均譜形狀相對P波段較為簡單,主要體現為單峰形式.相對3級海況,4級海況下的多普勒譜展寬更大,尤其在頻率?300—?50 Hz仍具有一定幅度的功率密度,這與海面白冠散射能夠引起多普勒譜的展寬效應是對應的.說明在較高海況下,海面出現白冠散射,且該散射機制在海雜波中占有相當的比重.

圖7還給出了兩種譜模型對于上述數據建模的對比結果.從圖7可以看出,由于譜形狀表現為單峰結構,Walker譜模型和時變譜模型的建模結果基本一致,均能實現對譜線結構的較準確的描述,僅在一些細微結構處時變譜模型顯得更為精確.在該波段下,時變譜模型的優勢一般.

圖6 岸基P波段海雜波平均多普勒譜建模效果(a)3級海況;(b)4級海況Fig.6.The modeled long-time averaged doppler spectra at P band for(a)sea state 3 and(b)sea state 4.

圖7 岸基S波段海雜波平均多普勒譜建模效果(a)3級海況;(b)4級海況Fig.7.The modeled long-time averaged doppler spectra at S band for(a)sea state 3 and(b)sea state 4.

表3列出了與圖6和圖7對應的兩種譜模型在平均多普勒譜建模中的模型參數與定量的建模誤差.根據模型參數可以分析各個譜分量對總多普勒譜的貢獻,例如,對于P波段3級海況下時變多普勒譜的建模結果而言,IBh＞IW＞IS說明海雜波數據中布拉格散射的貢獻依次大于白冠散射和破碎波散射的貢獻,wW＞wB＞wS說明白冠散射產生的譜展寬依次大于布拉格散射和破碎波散射產生的譜展寬.對比建模誤差可以看出,時變譜模型的建模誤差總體上小于Walker模型,其中在P波段優勢更為顯著,在S波段兩者建模效果相當,這是由兩個波段平均譜的特征決定的,與圖6和圖7結論一致.

表3 兩種譜模型在平均譜建模中的建模參數與建模誤差Table 3.The parameters and RMS errors of the two spectrum models in modeling the long-time averaged doppler spectra of sea clutter.

5 結 論

本文在海雜波傳統多普勒譜模型的基礎上,提出了一種時變多普勒譜模型.通過對兩個波段海雜波短時多普勒譜的建模驗證,表明時變譜模型能夠更清晰地描述短時譜的各個譜分量,建模精度優于Ward短時譜模型,尤其體現在短時譜的觀測時間較長的情況下.在平均多普勒譜建模中,當實測譜的結構體現為多峰,形狀較為復雜時,Walker平均譜模型忽略掉了部分譜形狀的細微結構,而時變譜模型可以對整個譜形狀較為準確地建模.當實測平均譜的形狀較為簡單時,兩種譜模型性能基本一致.

[1]Ward K D,Watts S 2010IET Radar,Sonar Navig.4 146

[2]Li D C,Shui P L 2016IET Radar,Sonar Navig.10 699

[3]Rosenberg L,Watts S 2017 IET Radar,Sonar Navig.11 1340

[4]Pidgeon V W 1968 J.Geophys.Res.73 1333

[5]Lee P H Y,Barter J D,Beach K L 1995 IEE Proc.Radar,Sonar Navig.142 252

[6]Walker D 2000 IEE Proc.Radar,Sonar Navig.147 114

[7]Rosenberg L 2014 IEEE T.Aero.Elec.Sys.50 406

[8]Watts S,Rosenberg L,Bocquet S 2016 IET Radar,Sonar Navig.10 24

[9]Watts S,Rosenberg L,Bocquet S 2016 IET Radar,Sonar Navig.10 32

[10]Corretja V,Meslot V,Kemkemian S,Montigny R 2017 Proceedings of the IEEE Radar Conference Seattle,WA,USA,May 8–12,2017 p1046

[11]Watts S,McDonald M,Rosenberg L 2015 Proceedings of the IEEE Radar Conference Johannesburg,South Africa,October 27–30,2015 p91

[12]Jiang W Z,Yuan Y L,Wang Y H,Zhang Y M 2012 Acta Phys.Sin.61 124213(in Chinese)[姜文正,袁業立,王運華,張彥敏2012物理學報61 124213]

[13]Walker D 2001 IEE Proc.Radar,Sonar Navig.148 73

[14]Watts S 2012 IEEE T.Aero.Elec.Sys.48 3303

[15]Miller R J,Dawber W N 2002 Proceedings of the International Radar Conference Edinburgh,UK,October 15–17,2002 p444

[16]Ward K D,Tough R J A,Watts S 2013 Sea Clutter:Scattering,the K Distribution and Radar Performance(2nd Ed.)(London:The Institution of Engineering and Technology)

[17]Raynal A M,Doerry A W 2010 Doppler Characteristics of Sea Clutter(New Mexico:Sandia National Laboratories)

[18]Duncan J R,Keller W C,Wright J W 1974 Radio Sci.9 809

[19]Plant W J,Keller W C 1990 J.Geophys.Res.95 16299

[20]Zhang Y M,Wang Y H,Zhao C F 2010 Chin.Phys.B 19 084103

[21]Zhang Y S,Zhang J P,Li X,Wu Z S 2014 Chin.Phys.B 23 108402