基于顯著性特征的大冰雹識別模型*

王萍 潘躍

1 引言

冰雹災害是由強對流天氣系統引起的一種劇烈的氣象災害,它出現的范圍較小,時間短促,但來勢猛烈,并常常伴隨著狂風、強降水、急劇降溫等陣發性災害性天氣過程.新一代天氣雷達是監測甚至預測包括冰雹在內的強對流天氣的重要平臺,其中最具代表性的是美國于1996年規模化應用于全美的配有較豐富應用軟件的多普勒天氣雷達氣象業務系統,該系統所提供的冰雹指數與所替代的冰雹指數相比,使美國強冰雹預警的擊中率維持在70%的同時,TS評分從26%增至42%[1-3].目前,新一代天氣雷達及其應用軟件算法與美國基本相同,該系統每隔6 min提供一組多個探測仰角下的回波圖像,特別是其反射率因子圖和徑向速度圖能夠直觀展現對流云團的結構、形態以及對流場的分布,從而促進了人們對短時強降水、強冰雹、災害性大風等不同類型的強對流天氣特點及其演變規律的認識.就強冰雹而言,人們認識到：

1)凡發生冰雹的云團,至少有一部分越過了環境融化層,即零度層,越過的部分越多、越過的程度越高、零度層的高度(與地表溫度有關)越低,降強冰雹的可能性越大,反映在雷達回波的反射率因子圖上,高仰角的冰雹回波強度會持續很高;

2)在強冰雹云團回波的反射因子圖上,通常會呈現一種“懸垂”結構,即位于中高層的較強回波區域會“探”到其低層較強回波區域之“外”,換言之,即在底層弱回波區域之上“懸垂”著一部分較強的回波體[4];

3)強冰雹云團的低層回波出現明顯的梯度不均等結構,即從回波高強度中心遞減到低強度邊界所經距離的差異極大,導致強回波中心偏倚于回波圖像的一側[5],上述懸垂現象就多發生于該低層回波圖像的極高或較高梯度的方位;

4)與短時強降水回波相比,強冰雹的高強度回波范圍所占比例更大[6].

另外,還認識到一些評定強冰雹的充分條件,即回波的反射率因子圖像一旦出現“三體散射”或“鉤”、徑向速度圖像一旦出現“強中氣旋”,則肯定出現強冰雹[7].

目前,氣象上普遍使用的業務系統提供的強冰雹指數算法是考慮上述特點1)設計的,因為實際上許多短時強降水回波也會有一部分甚至大部分擴展到零度層高度以上,因此,基于這一強冰雹指數的強冰雹誤報(空報)率會很高,據美國給出的統計,空報率高達49%[8].在國內,由于這一過高(經常會高于49%)的空報率,許多氣象臺只能將其作為強冰雹的起報條件,有的甚至直接用其輔助預報“短時強降水”[9].為此,人們展開了一些積極的研究和探索,其中,所提出的“指標加權法”[10]增加了液態水含量、中氣旋、回波形態、三體散射等多個預報因子,通過對各預報因子分段賦值再求和的方法,獲得對冰雹的TS評分的測試結果達到75%;“相似演變聚類法”[11]雖然還是主要選擇回波高度及強度類指標作因子,但通過與已發生的冰雹過程貼近度的計算,使對冰雹的TS評分達到65%等.

本文在考慮冰雹指數的同時,新增了反映上述特點2)—4)的特征提取算法,并以短時強降水回波為反例,在對諸項特征顯著性分析的基礎上,聯合使用這些特征訓練出強冰雹識別模型.推出一種新型冰雹指數,最后對新的冰雹指數與傳統的冰雹指數進行了對比測試.

2 強冰雹回波的特征提取

根據氣象理論和大量觀測經驗可知,發生強對流天氣時的多個仰角的多普勒雷達回波圖中會出現顯著的強“單體”,即在多個仰角的反射率圖像中出現以大于等于50 dBZ的區域為核,從“核”向外反射率逐漸由高(50 dBZ)到低(25 dBZ)過渡的區域.本文首先以多仰角的反射率回波單體為對象,根據上小節指出的冰雹回波特點,構建用于鑒別冰雹回波的系列特征.

2.1 多仰角單體回波的匹配

在強對流可能出現的時段,多普勒天氣雷達多采用9仰角的掃描模式,通過9張回波強度圖像捕獲對流體的立體結構信息,因此,在多單體共存的通常情況下,需要確定各單體在最多9張仰角圖中的匹配關系.

2.1.1 摒棄高度信息的重疊率規則

圖1是雷達通過兩仰角方式探測垂直線型對象A的示意圖.在探測仰角α1下的錐面探測圖像I1中,探測A為點P1(r1,θ1,h1),同理,在探測仰角α2下,探測A為點P2(r2,θ2,h2),將點P1和點P2向水平面投影,得到P1′(r1′,θ1′)和P2′(r2′,θ2′),其中容易推知：θ1=θ2,r′1=r′2,即點p′1和點p′2完全重合.推而廣之,若對象A為任意水平截面形態完全相同的垂直物體,則在其錐面圖像I1和I2的探測結果雖不相同(s1/=s2),但其水平投影依然完全重合 (s′1=s′2). 一個強對流云體幾乎是垂直的,但不同高度的形態、面積大小會有差異,強冰雹的“懸垂”結構還會使這種差異變大,因此,在對同一單體的多仰角探測圖像進行匹配時,將重合條件放寬至小于1的一個比率ρ,即來自相鄰兩仰角錐面中對象的水平投影區域重疊率超過ρ,則認為兩者屬于同一單體.

圖1 垂直對象水平投影示意圖

2.1.2 單體的縱向匹配算法

1)對投影后的九個仰角反射率因子圖I1,I2,···,I9進行二值化處理,保留30 dBZ以上區域并進行形態學運算,除去孤立塊、填補孔洞、去除毛刺;

2)取每個連通的區域最小外包矩形Mpip,p=1,···,9,ip=1,···,np;

3)將高層Mpip逐層逐個與低層Mqjq進行比較,若對任意的q＞0,q＜p＜10都存在一個Mqj′

q使得

其中,Mα為常數,這里取0.2,記Mk=M′k={M1j′1,M2j′2,···,Mpj′p},k=1,···,n(n為單體數);并將滿足(2)式的Mqj′q作標記;

4)用未進行比較或未被標記的Mpip重復3),直到所有Mqjq都被標記;

5)若最低仰角反射率因子圖I1中有未被標記的M1i1,直接將M1i1記為Mk;

6)若低層的同一連通區域的最小外包矩形被多次標記,如圖2所示,低層連通區域的最小外包矩形a與其上層兩個連通區域的最小外包矩形b和c都匹配,此時a被標記兩次,則取消相關標記,提高分割閾值,對多次標記的矩形框內區域重新進行閾值分割,返回4),直至各層各單體都被一一匹配.

圖2 低層同一區域被上層的多區域標記示例

2.2 傳統冰雹指數簡介

首先,根據單體在某一仰角反射率圖像上的強度值,計算與地面冰雹災害潛勢密切相關的“冰雹動能”E[12]

其中,Z為單體在該仰角圖像上的最大反射率值,W(Z)為權重

該權重使最大反射率強度與冰雹動能正相關,并使最大反射率值低于ZL的回波區域不對“冰雹動能”做出貢獻.

再求環境融化層(零度層)H0以上的冰雹動能加權通量：

其中權重函數

表明,單體躍出零度層高度的程度(H-H0)與總冰雹動能正相關,同時認為,低于H0的冰雹動能不會對降雹做出貢獻.該冰雹動能加權通量就是傳統的強冰雹指數,在用其進行強冰雹監測和預警時,通常使用根據美國Oklahoma州和Florida州風暴實況的相關分析得出的預警閾值函數：

即若SHI＞WT,則進行強冰雹預警[13].

可見,傳統強冰雹指數屬于反映強冰雹回波特點(1)的一種特征.

2.3 單體的懸垂度特征

根據氣象預報經驗,當對流風暴的低仰角反射率因子圖具有寬闊的弱回波區或有界弱回波區,且在它們上方存在強反射率因子核時,最有利于強冰雹的發生.圖3(a)—(c)給出了單體非強風暴、單體強風暴和超級單體風暴反射率因子垂直結構示意圖,實線為單體低層回波強度的等值線,虛線為中層回波強度大于20 dBZ的輪廓線,實心點為高層最大回波強度處.可見,具有冰雹潛勢的單體強風暴(圖3(b))和超級單體風暴(圖3(c)),其反射率因子結構在低層高梯度側呈現強回波“懸垂”.

圖3 風暴單體反射率因子垂直結構示意圖 (a)單體非強風暴;(b)單體強風暴;(c)超級單體風暴

2.3.1區域平均梯度

為定位出懸垂關注區,需要尋找強風暴在低層回波圖像中反射率快速及較快速遞減的區域,也就是高梯度區.

圖4 從RU過渡到RD的局部區域示意圖

圖4 是某單體低層回波從高反射率RU等值線過渡到低反射率RD等值線的局部情況,在關注區域 A 中,從點 pi到 p′i的梯度

因為(RU-RD)為常數,G(pi,p′i)與 ‖pi-p′i‖ 成反比,為降低運算時間開銷,特別定義區域A從RU過渡到RD的平均梯度為

即

其中,SA為區域A的面積.可見區域A的平均梯度與其平均寬度成反比,若設區域A的長度為常數,則區域A的平均梯度與其面積成反比,因此,尋找最大梯度區,等價于尋求平均寬度最小區域或最小面積區域.

2.3.2 懸垂區域關注系數

考慮到出現“懸垂”的高梯度區具有一定的長度,且高梯度區一般不會突變到最低或次低梯度區,首先對單體進行處理,保留RU至RD的區域 (缺省值：RU=45 dBZ,RD=30 dBZ),如圖 5(a)—(c)所示,然后以單體的強反射率核心為中心作射線將單體分為 8份 A1,A2,···,A8(見圖 5(d)),分別求出它們的面積 S1,S2,···,S8.一般情況下,區域面積越小,其平均梯度越大,將其按照從小到大的順序排序,考慮到較大梯度區域不一定剛好落入一個分區內以及懸垂一定不會出現在最低及次低梯度區之上,設懸垂區域關注系數β={0.3,0.3,0.2,0.1,0.05,0.05,0,0}.

2.3.3 懸垂度提取算法

1)一個匹配好的風暴云團單體k,經多仰角匹配后得到外包矩形集合Mk={Mk1,Mk2,···,MkN},其中k=1,···,K,K為一個時刻回波圖像中出現的單體數,對每一個k執行以下過程(為簡單起見,略去k不記).

圖5 八射線區域分割示意圖 (a)單體低仰角反射率圖;(b)保留反射率值大于等于RD的區域;(c)減去(b)區域中反射率值大于等于RU的區域;(d)用八條射線將(c)區域等分

2)對M1,求取梯度區

將ΔA用射線等分成8份,并求取相應的面積ΔSi,對八等份區域根據ΔSi的大小按升序編號n,n=1,···,8.

3)對 Mj,j=1,···,N,做與 2)相同的八等份劃分,得到Mji中的有效回波面積Sji,并賦予編號,即Sji(n).

4)計算第m層回波對m-1層回波的局部懸垂度

5)計算總懸垂度

2.4 高回波比

雹暴通常與大片的強雷達回波相聯系,高回波比用于反映此項特點,定義為50 dBZ及以上區域占單體30 dBZ及以上區域的比重.該特征不涉及高度信息,因此可借助單體匹配成功的外包矩形集合 M={M1,M2,···,MN},N ≤ 9 求得.

設矩形區域Mi中,取值大于等于50 dBZ的像素數為Ni-50,取值大于等于30 dBZ的像素數為Ni-30,則定義高回波比特征如下：

2.5 峰度特征

圖6(a),(b)是某個短時強降水單體和某個強冰雹單體的低仰角(0.5°)反射率因子圖,圖7是它們的強度分布直方圖.直觀地看,兩者有著高強度值成分略少、低強度比例相差不大、單側尾部較厚的分布規律.基于分布直方圖的“峰度”是一個突出“厚重尾部”的4階統計量,定義式如下：

其中,ˉY,s是其均值和標準差.因為強冰雹的最高反射率普遍高于短時強降水,或者說強冰雹的較高反射率所占份額普遍高于短時強降水,因此強冰雹的均值普遍高于短時強降水的均值,從而導致強冰雹的“峰度”普遍高于短時強降水的“峰度”.

圖6 強冰雹與短時強降水反射率圖對比 (a)短時強降水;(b)強冰雹

圖7 強冰雹單體和短時強降水單體反射率強度分布直方圖

另外,本文還設計了有效厚度特征ΔH,即大于45 dBZ的單體核的最大高度與零度層高度H0之差：

2.6 特征的顯著性分析

為檢驗所提特征的顯著性,隨機選擇來自14過程的強冰雹單體60個和來自12個過程的短時強降水單體52個進行測試,形成基于單個特征的直方圖分布,如圖8所示.

圖8 關于強冰雹樣本和短時強降水樣本的各特征統計直方圖 (a)傳統冰雹指數統計直方圖;(b)懸垂度統計直方圖;(c)高回波比統計直方圖;(d)峰度統計直方圖;(e)有效厚度統計直方圖

設用于測試特征有效性的強冰雹樣本和短時強降水樣本分別來自兩個來自正態分布的總體,且方差相同.下面根據測試數據(圖8)推斷由相關特征描述的這兩個總體有無顯著性差異,并由此推斷所提特征的有效性,即假設強冰雹均值等于短時強降水均值,使用(17)式的服從t分布的統計量展開置信水平(1-α)的檢驗.計算結果如表1所示.

式中,ˉx,ˉy分別為強冰雹和強降水樣本的特征均值,S21,S22分別為對應的方差.

取顯著性水平α=0.01,則tα/2(n1+n2-2)=t0.005(112),查表得t0.005(112)＜2.660.由于5個特征計算結果均使 ti＞ 2.66 ＞ t0.005(112),i=1,···,5,因此推翻原假設,認為本文中所選各項特征在短時強降水和強冰雹之間均具有顯著性差異.

表1 5個特征的單因子顯著性差異檢驗

3 強冰雹識別模型

上小節所提5項特征均是站在強冰雹回波特點的角度上設計和實現的.圖8和表1所展現的實測結果表明,它們在短時強降水和強冰雹兩類風暴之間的統計性差異還是很明顯的.另外,根據氣象理論,零度層高度也是左右冰雹強弱的重要特征.因此,本文特別將短時強降水這種極易與強冰雹混淆的單體作為反例,訓練基于這5項顯著性特征及零度層高度的分類模型,并選用支持向量機擔當此任.

3.1 支持向量機分類模型

支持向量機(SVM)是根據結構風險最小原理,在特征空間或特征變換空間尋找一個最優分類超平面,使得兩類樣本中離最優分類超平面最近的樣本點到超平面的距離最大,當這個超平面位于特征空間時,該支持向量機是線性的,位于特征變換空間時則為非線性的.例如圖9最優超平面為H,直觀地看,離開超平面H越遠(即di越大)的樣本,其特征與另一類相差越大,就越屬于這個類型.

圖9 線性可分SVM示意圖

設線性可分的樣本集 (xi,yi),i=1,2,···,n;其中xi的類別編號yi={1,-1},超平面方程為

使得

在面對非線性分類問題時,只需將(19)式中的點積(xi·x)置換成核函數K(xi,x),即可將其轉化為某個隱性的高維空間中的線性問題,再在該高維空間中得到最優分類超平面[14].在這時的約束條件上添加一個松弛變量ξi≥0,使優化問題變為

其中懲罰因子C＞0.

非線性支持向量機的最優判斷函數為

文獻[15]研究表明,利用基于升維思想的支持向量機方法對來自實際氣候系統的非平穩過程存在穩定的預測能力.文獻[16,17]研究表明,結合聚類技術、模糊邏輯的支持向量機在解決混沌時間序列的預測問題時顯現優勢.

3.2 基于支持向量機的強冰雹識別模型

選用冰雹指數SHI,懸垂度G,高回波比ρH,峰度μ4,有效厚度ΔH和零度層高度H0六項特征,從天津地區2002—2008年間26次風暴過程的多普勒雷達數據中選出120組學習樣本(強冰雹(d≥20 mm)60組,短時強降水(降水量≥20 mm/h)60組),利用libsvm庫文件編寫代碼,訓練強冰雹識別模型,過程如下.

1)歸一化處理：設特征向量

對Xi,i=1,2,···,5,使用對各個特征測試過程中獲得的Ximin和Ximax;對X6,根據零度層高度一般處于2—5 km之間的經驗,令X6min=2,X6max=5,利用 (24)式,對樣本 x=(x1,x2,···,x6)T的 6 個分量歸一化：

2)確定核函數：采用最常用的徑向基函數作為支持向量機的核函數,表達式為

3)訓練樣本以獲得模型參數：通過對樣本的訓練,得到最優懲罰參數C=32,最優核參數γ=0.5.

4)獲得樣本點到超平面距離d.

3.3 構建新冰雹指數

通過訓練樣本得到的由為數不多的支持向量(樣本)獲得的超平面,將隱性的高維空間一分為二,強冰雹樣本均落入該超平面的一側,直觀地看,離開該超平面越遠的強冰雹樣本與其反例(短時強降水)樣本越不像.換句話說,遠離超平面的樣本較超平面附近的樣本其強冰雹屬性更強.因此,文本借助(26)式將此距離歸一化,并稱之為新冰雹指數：

(26)式中,參數C為決定μ變化到±1快慢的整數,根據2＜dmax＜3的實驗結果,令

解得C＞1.2,因此取C=2.

4 實驗測試

4.1 有效性測試

從京津冀地區2002—2010年期間17個強冰雹過程中隨機選取未參加訓練的單體樣本33個、15個短時強降水過程中選取未參加訓練的單體樣本31個、7個小冰雹過程中選取單體樣本13個,組織對本文提出的新冰雹指數與傳統的冰雹指數的對比性測試.圖10給出了所有測試樣本的傳統冰雹指數與零度層高度的散點圖,其中直線WT=57.5H0-121是基于零度層高度的閾值函數;圖11是這些樣本關于新冰雹指數的分布直方圖,其中橫坐標是對新冰雹指數按照四舍五入的方式保留1位小數的整定結果.表2是兩種方法關于強冰雹擊中率POD、誤報率FAR以及臨界成功指數CSI(即TS評分)的統計結果,其中擊中率、誤報率及臨界成功指數定義如下：

4.2 時效性測試

選擇發生于河南、廣東的4個強冰雹過程的5組體掃數據,組織對本文提出的新冰雹指數與傳統冰雹指數的時效性測試,測試樣本的冰雹發生時間等詳細信息及測試結果見表3.容易看出,對于5組體掃數據所展現的4個強冰雹過程,本文方法全部識別正確,且識別時間較強冰雹發生時刻至少提前了30 min;而對于相同的測試過程,PUP所提供的傳統冰雹指數僅對河南的強冰雹過程的識別結果與本文相似,其中對于河南濮陽雷達的探測過程,PUP的預警提前量比本文算法提前了2個體掃12 min,至于廣東的3個過程,PUP的一個預警提前量大幅度減小、一個僅擊中前11個體掃且于強冰雹尚未發生時就出現了錯誤識別、第三個過程的提前識別量從本文方法的30 min減小到6 min.

表2 -1 以短時強降水單體做反例的強冰雹識別

表2 -2 以短時強降水和小冰雹單體做反例的強冰雹識別

圖10 所有測試樣本關于傳統冰雹指數與零度層高度的散點圖

圖11 所有測試樣本關于新冰雹指數的分布直方圖

表3 時效性測試結果

5 分析和小結

1)本文設計了4個新的基于風暴回波強度(反射率因子)圖像的特征,連同傳統的冰雹指數,力圖從不同的側面去貼合強冰雹的物理概念,體現強冰雹單體回波的結構特點,圖8給出的測試結果表明,就單個特征而言,它們對強冰雹和強降水單體的區分能力均較高,但因在兩類間有一定的重疊又不盡善盡美.就傳統冰雹指數而言,許多短時強降水由于其強回波頂也會較高(圖8(e),約7成的短時強降水的強回波頂越過零度層高度4 km以上),因而獲得的冰雹指數值也較大,這就是為什么雖然傳統強冰雹指數在具有較高的強冰雹擊中率(93.93%)的同時,會伴隨著很高的誤報率(32.61%)的原因.

2)本文在關注強冰雹單體自身特點的同時,注意到它們與非強冰雹單體量上的差異,令從不同角度抽取的多個顯著性特征在模式識別的概念下協同工作,加上支持向量機這種適應性極高的非線性分類器,使得本文歸納出的新冰雹指數呈現出優良的工作效果(表2-1：POD=96.97%,FAR=3.13%,CSI=91.18%).

3)因為本文提出的新冰雹指數實質上是從由強冰雹和短時強降水這兩類對象的非線性分類器的判斷結果轉化而來,這個分類器應該同時具有識別短時強降水的能力,站在識別強降水的角度得到統計結果如表4所示.可見此模型識別降水的能力依然很強(POD=96.77%,FAR=3.23%,CSI=93.75%).

表4 短時強降水識別

4)因為決定本文新冰雹指數的“懸垂度”、“有效厚度”、“高回波比”等5個關鍵特征于強冰雹的形成期間就逐漸顯著起來,因此新冰雹指數具有識別出正在形成的(即尚未發生的)強冰雹的能力,表3給出的測試結果表明,這一提前識別的時效性優于傳統的冰雹指數.

5)按照模式識別的概念,判別模型只適用于對與訓練樣本相同類型的未知樣本進行識別,而強冰雹、短時強降水與弱冰雹乃至弱降水同屬對流天氣,在線識別時,弱冰雹等是不可回避的識別對象,如果不顧及它們的存在,定會造成在線識別的誤報率升高(表2：FAR從3.13%升至17.95%;表4：FAR從3.23%升至21.05%).

分析圖11容易發現,就預報指數而言,有46.15%的小冰雹落入了強冰雹區,余下的53.85%落入了短時強降水區,這應該是提取特征時沒有考慮小冰雹單體特點的緣故.另外,我們注意到小冰雹與短時強降水的指數幾乎不發生重疊,對此,本文提出如下解決思路：

1)利用關鍵特征(如強回波比)的閾值將弱降水剔除;

2)增加強降水的特有特征,如液態水含量、時間累積等,將新冰雹指數為負的小冰雹乃至誤報為短時強降水的強冰雹從短時強降水中剝離出來;

3)對弱冰雹和強冰雹展開新特性的設計及其差異性研究,如中氣旋及其強度等,尋求從強冰雹中剝離出弱冰雹的有效方法;

4)添加算法,通過“三體散射”、“強中氣旋”等的識別,首先將滿足強冰雹充分條件的單體挑選出來;

5)建立分層識別模型,第一層剔除弱降水,第二層識別滿足充分條件的強冰雹,第三層使用本文算法將待識別單體暫時歸為強冰雹或短時強降水,第四層進行是否為弱冰雹的鑒別.

以上思路也正是本文工作的進一步的延續,同時本文作者將進一步關注和致力解決冰雹尺寸的預測問題.

[1]Makitov V 2007 Atmo.Res.83 380

[2]Seraf i n R J,Wilson J W 2000 Bull.Amer.Meteorol.Soc.81 501

[3]Holleman I,Wessels H R A,Onvlee J R A,Barlag S J M 2000 Phys.Chem.Earth B-Hydrol.Oceans Atmos.25 1293

[4]Lopez L,Sanchez J L 2009 Atmos.Res.93 358

[5]Baeck M L,Smith J A 1998 Weather Forecast.13 416

[6]Edwards R,Thompson R L 1998 Weather Forecast.13 277

[7]Lemon L R 1998 Weather Forecast.13 327

[8]Doviak R J,Bringi V,Ryzhkov A,Zahrai A,Zrnic D 2000 J.Atmos.Ocean.Technol.17 257

[9]Liang Q Q,Lin L X 2008 Meteorol.Sci.Technol.2 150(in Chinese)[梁巧倩,林良勛2008氣象科技2 150]

[10]Zhang G C 2011 Strong Convection Weather Analysis and Forecasting(1st Ed.)(Beijing：China Meteorological Press)p186(in Chinese)[章國材2011強對流天氣分析與預報(第一版)(北京：氣象出版社)第186頁]

[11]Zhou Y F,Zhu Y J 2007 J.Anhui Agricult.Sci.35 9637(in Chinese)[周葉芳,朱擁軍2007安徽農業科學35 9637]

[12]Witt A,Eilts M D,Stumpf G J,Johnson J T,Mitchell E D,Thomas K W 1998 Weather Forecast.13 286

[13]Lenning E,Fuelberg H E,Watson W I 1998 Weather Forecast.13 1029

[14]Zhang J S,Dang J L,Li H C 2007 Acta Phys.Sin.56 67(in Chinese)[張家樹,黨建亮,李恒超2007物理學報56 67]

[15]Wang F F,Zhang Y R 2012 Acta Phys.Sin.61 084101(in Chinese)[王芳芳,張業榮2012物理學報61 084101]

[16]Cai J W,Hu S S,Tao H F 2007 Acta Phys.Sin.56 6820(in Chinese)[蔡俊偉,胡壽松,陶洪峰2007物理學報56 6820]

[17]Cui W Z,Zhu C C,Bao W X,Liu J H 2005 Acta Phys.Sin.54 3009(in Chinese)[崔萬照,朱長純,保文星,劉君華2005物理學報54 3009]