熊雄,楊日杰,韓建輝,郭新奇
(1.海軍航空工程學(xué)院 電子信息工程系,山東煙臺264001;2.海軍航空工程學(xué)院指揮系,山東 煙臺264001)
在潛艇聲隱身性能越來越好的條件下,航空磁探儀由于幾乎不受傳播介質(zhì)特性的影響,成為一種有效的反潛探測設(shè)備,在航空反潛中得到廣泛應(yīng)用。雖然航空磁探潛受海洋環(huán)境復(fù)雜傳播介質(zhì)的影響較小,但是隨著航空磁探儀靈敏度越來高,潛艇磁異常信號也越來越容易受到各種背景擾動的干擾。在反潛機(jī)飛行高度較低的情況下,海洋中的風(fēng)浪等海水運(yùn)動切割地球磁場,激發(fā)感應(yīng)電磁場,對航空磁探儀的工作性能產(chǎn)生重要的影響。大量的實(shí)驗(yàn)表明海浪產(chǎn)生的感應(yīng)電磁場噪聲與要檢測的目標(biāo)磁場量級、頻帶基本接近,這些噪聲都是不可忽略的干擾源[1-2]。然而海浪磁噪聲很難直接測量,基于線性波浪理論對海浪磁場進(jìn)行建模和仿真是研究海浪感應(yīng)磁場噪聲特點(diǎn)及能量分布特性等因素的一種重要方法[3]。
Weaver在海浪波高為常數(shù)的假設(shè)下,給出了理想海域條件下單頻重力波產(chǎn)生的感應(yīng)磁場的理論表達(dá)式,并形成了經(jīng)典 Weaver海浪磁場模型[4]。Ochadlick對Weaver海浪磁場模型進(jìn)行驗(yàn)證,驗(yàn)證結(jié)果表明該模型具有較強(qiáng)的可信性[5]。近年來,有學(xué)者分別基于Weaver海浪磁場模型給出了有限深海域磁場模型和海浪磁場矢量理論模型以及計(jì)算方法[6-7]。但是以上模型和計(jì)算方法都是基于單頻海浪重力波,而實(shí)際中海浪是一種復(fù)雜的海水運(yùn)動,其浪高隨頻率變化,實(shí)際應(yīng)用中描述海浪特征的方法是海浪譜分析,完整的海浪譜由頻率譜和方向函數(shù)組成[8]。文獻(xiàn)[9]基于海浪譜推導(dǎo)了航空磁探儀接收到的海浪磁噪聲功率譜的理論表達(dá)式,文獻(xiàn)[10-11]給出了基于海浪頻率譜等分法的海浪磁場數(shù)值仿真方法,但是該方法仿真速度慢,而且沒有考慮方向函數(shù)的影響。文獻(xiàn)[8]在頻率等分法的基礎(chǔ)上采用傅里葉反變換來模擬海浪,該方法可以減少仿真次數(shù),提高仿真速度,但是該方法無法應(yīng)用到海浪磁場仿真。
針對以上問題,本文基于Weaver單頻波海浪磁場模型,建立地理坐標(biāo)系下任意方向傳播實(shí)際海浪磁場數(shù)值模型,給出快速數(shù)值仿真算法,實(shí)現(xiàn)海浪磁場的快速仿真,并對仿真結(jié)果進(jìn)行校驗(yàn)和分析。
根據(jù)Weaver海浪磁場理論,在地磁場中海洋重力波感應(yīng)電場E,磁感應(yīng)強(qiáng)度B滿足麥克斯韋電磁理論:
式中:μ為海水磁導(dǎo)率,ε為海水介電常數(shù)。海水電流傳導(dǎo)密度為J=σ(E+V×BE),σ為海水電導(dǎo)率,BE表示地磁場矢量,V為海洋重力波速度矢量。
通常認(rèn)為海水傳導(dǎo)電流密度遠(yuǎn)大于式(2)中等號右邊第2項(xiàng)的位移電流,因此可以忽略位移電流。則海浪感應(yīng)磁場B可以表示為[11]
為了求解B,Weaver引入速度勢φ,且定義V=?φ,并求解得到沿X軸方向傳播單頻波海浪磁場強(qiáng)度。本文將建立地理坐標(biāo)系下任意方向傳播海浪磁場強(qiáng)度數(shù)學(xué)模型。
建立地理直角坐標(biāo)系OXYZ,OXY位于平均海平面,OZ垂直向上。OW為海浪傳播方向,OW與OX軸的夾角為θ,航空磁探儀在海平面上方高度Zm沿著OM直線飛行,飛行路徑OM與OX軸的夾角為β,如圖1所示。Z>0為空氣介質(zhì),Z<0為海水介質(zhì)。ON為磁北方向,地磁場矢量BE表示為BE=,I表示磁傾角,γ表示磁北方向與X軸的夾角,如圖2所示。
圖1 地理坐標(biāo)系Fig.1 Geographic coordinate system
圖2 地磁場矢量示意圖Fig.2 Schematic of the geomagnetic vector
假設(shè)海水是不可壓縮無旋流體,根據(jù)文獻(xiàn)[10],波浪沿θ方向傳播,以簡諧運(yùn)動描述單頻波海浪流體運(yùn)動,則流體擾動速度勢可以表示為
式中:Ω=xcos θ+ysin θ,a、ω 、k分別表示單頻波幅度、頻率、波數(shù),g為重力加速度,k和ω的散布關(guān)系可以表示為
將式(4)代入式(3)求解得到坐標(biāo)點(diǎn)(x,y,z)處t時刻單頻波海浪磁場B(x,y,z,t)為
式中:hB(z,θ)為磁場幅度矢量。
海浪磁場傳播經(jīng)過海水和空氣兩層介質(zhì)。根據(jù)邊界條件z=0處海浪感應(yīng)磁場的垂直分量連續(xù),并且結(jié)合式(7)可以得到磁場幅度標(biāo)量
其中:
則t時刻海平面上方坐標(biāo)點(diǎn)(x,y,z)處標(biāo)量磁探儀探測到單頻重力波標(biāo)量海浪磁場為
式中:ε為海浪初始相位,在(0,2π)上均勻分布。
海浪是一種復(fù)雜的隨機(jī)運(yùn)動過程,在海洋學(xué)研究中利用隨機(jī)過程來描述海浪是進(jìn)行海浪研究的主要途徑之一。為了模擬實(shí)際的海洋環(huán)境,根據(jù)Longuet-higgins線性波浪理論,t時刻位于坐標(biāo)點(diǎn)(x,y)處波面 ζ(x,y,t)可以表示為無限個隨機(jī)相位正弦波的疊加:
式中:εn為第n個組成波的初始相位,an、ωn、kn第n個組成波的幅度、頻率、波數(shù),θ為波浪主傳播方向。
根據(jù)線性理論,海浪產(chǎn)生的磁場也可以表示為無限個單頻重力波海浪感應(yīng)磁場的疊加。由式(9)、(10)可以得到磁探儀靜止條件下探測海浪磁場為
其中:
在航空磁探測過程中,磁探儀是隨著飛機(jī)運(yùn)動的,因此其接收到的海浪磁噪聲不僅隨時間變化而且隨觀測位置變化。根據(jù)式(11)可以得到以速度v按照航向β飛行的航空磁探儀探測海浪磁場為
其中:
式(11)、(12)給出了實(shí)際海浪磁場的數(shù)值計(jì)算模型,但是模擬實(shí)際的海浪磁場需要根據(jù)數(shù)值模型給出有效的數(shù)值仿真算法。由式(11)、(12)可知,海浪磁場與波浪頻率以及相應(yīng)頻率的波高有關(guān),而且海浪具有三維不規(guī)則性,不僅波高不同、頻率不同,而且會從各個方向傳到某一點(diǎn),除沿主風(fēng)向產(chǎn)生的主浪以外,在主浪向兩側(cè)±π/2角度范圍內(nèi)都有諧波的擴(kuò)散。描述海浪三維不規(guī)則性常用的方法是海浪譜[12-13]。
海浪譜定義為單位頻率間隔和方向間隔內(nèi)的海浪平均能量密度,二維海浪譜又叫方向譜,方向譜S(ω,θ)是由頻率和角度相關(guān)的2個函數(shù)組成,可表示為
式中:S(ω)為海浪頻譜,G(ω,θ)為海浪方向分布函數(shù)。
海浪頻譜比較容易觀測,國外根據(jù)大量海浪觀測資料,提出了許多的海浪頻譜模型。Person-Moscowitz譜模型簡稱P-M譜,能較好地描述風(fēng)速為0~20 m/s之間的海浪譜。P-M譜模型表示為
式中:S(ω)為能量頻譜,ω為頻率,α=0.008 1,β=0.74,g為重力加速度,U為海面上19.5 m處風(fēng)速,譜峰頻率為ωn=8.565/U。
根據(jù)ITTC的觀測資料方向分布函數(shù)可以簡化表示為
根據(jù)方向譜可以得到在 ωi-Δωi/2~ωi+Δωi/2頻段和θi-Δθi/2 ~ θi+Δθi/2角度內(nèi)海浪波高ai,j,可以表示為[12-13]
海浪磁場仿真算法具體步驟如下:
1)海浪頻段的選擇。為了提高仿真速度和仿真時間,需要對海浪頻段進(jìn)行估計(jì)。設(shè)定風(fēng)速,根據(jù)海浪譜表達(dá)式(13)對海浪譜的頻段進(jìn)行估計(jì),選擇有限的頻段ω1~ωn來數(shù)值計(jì)算。
2)進(jìn)行頻段和方向的離散化采樣。根據(jù)海浪譜密度函數(shù),對頻譜和方向進(jìn)行離散化采樣,頻率采樣間隔為Δω,方向的采樣間隔為Δθ。
3)計(jì)算每個離散網(wǎng)格上海浪波高。根據(jù)式(16)可以得到ωi和θj對應(yīng)網(wǎng)格下的海浪的波高ai,j。
4)產(chǎn)生隨機(jī)相位εn。利用隨機(jī)數(shù)生成原理產(chǎn)生[0,2π)之間均勻分布的隨機(jī)數(shù)。
5)單頻波海浪磁場的計(jì)算。設(shè)定坐標(biāo)位置點(diǎn)(x,y,z),根據(jù)式(9)計(jì)算 ωi和 θj對應(yīng)網(wǎng)格下的單頻重力波產(chǎn)生的磁場信號模值。
6)海浪磁場信號的合成。根據(jù)式(11)或式(12)進(jìn)行單頻波海浪磁場信號的線性疊加,計(jì)算磁探儀靜止或者運(yùn)動條件下接收到的磁噪聲信號。
一種簡單的采樣方法就是區(qū)間等分法:將頻率區(qū)間和方向區(qū)間分別進(jìn)行M、N等分,取固定大小的采樣子區(qū)域Δω×Δθ,使,將每個采樣子區(qū)域中心對應(yīng)的頻率和方向角作為單元波的頻率和方向,按照3.1節(jié)仿真算法將不同振幅和頻率的單頻波合成得到海浪磁場。區(qū)間等分法算法簡單、容易實(shí)現(xiàn),仿真過程中為了盡可能的精確,采樣數(shù)相對要大,時空消耗大,不適合在線計(jì)算。
四叉樹分解合成算法具有節(jié)省存儲空間,提高運(yùn)算速度等優(yōu)點(diǎn),適合于快速計(jì)算,廣泛應(yīng)用于地形學(xué)圖形繪制和圖像處理[14-17]。相對于其他多叉樹算法,四叉樹具有結(jié)構(gòu)簡單,檢索效率高的優(yōu)點(diǎn)。四叉樹分解的基本思想是將二維平面按4個象限進(jìn)行遞歸分割,直到子象限的數(shù)值符合設(shè)定的條件,從而得到一棵四分叉的倒向樹。四叉樹分解的示意圖如圖3所示。在四叉樹分解中,一個根結(jié)點(diǎn)有4個子結(jié)點(diǎn),這4個子節(jié)點(diǎn)按順序標(biāo)為東北(NE)、西北(NW)、西南(SW)、東南(SE)4個子區(qū)域,4個子區(qū)域?qū)⒃瓐D形區(qū)域四等分。依此判斷4個子區(qū)域是否滿足進(jìn)一步分解的條件,如果不滿足分解條件則,子圖形成為葉子節(jié)點(diǎn)并存儲該節(jié)點(diǎn),如果滿足分解條件,則子圖形成為根節(jié)點(diǎn)進(jìn)一步分解為4個節(jié)點(diǎn),依此遞歸循環(huán)直至分解結(jié)束。
按照四叉樹分解的思想,提出基于四叉樹分解的海浪磁場快速優(yōu)化仿真算法。其優(yōu)化方法是在3.1節(jié)海浪磁場數(shù)值仿真算法步驟2)中采用等能量四叉樹分解方法,其步驟2)可以分解為以下步驟:
1)設(shè)定每個網(wǎng)格最小采樣能量與總能量的比例PE。
2)根據(jù)步驟1)中選擇的頻段和角度范圍,將S(ω,θ)進(jìn)行四叉樹遞歸分解。
3)判斷每個子節(jié)點(diǎn)是否滿足該網(wǎng)格的能量小于或者等于設(shè)定比例,若不滿足條件則繼續(xù)分解該網(wǎng)格,若滿足條件,則結(jié)束分解并記錄葉子節(jié)點(diǎn)網(wǎng)格信息。
4)將每個網(wǎng)格葉子節(jié)點(diǎn)按照樹形鏈表結(jié)構(gòu)記錄,在后續(xù)的仿真過程中采用樹形鏈表遍歷的方法快速遍歷每個葉子節(jié)點(diǎn),得到步驟3)中所需的每個分解單元的信息。
圖3 四叉樹分解示意圖Fig.3 Schematic of quadtree division
基于海浪磁場快速仿真算法,對不同條件下海浪磁場進(jìn)行數(shù)值仿真計(jì)算,并分析仿真性能。基本的仿真條件為:地磁場BE=50 000 nT,地磁傾角為60°,地磁偏角為 10°,重力加速度為 9.8 m/s2,海水的磁導(dǎo)率為4π×10-7H/m,海水電導(dǎo)率為5 S/m,采樣頻率為10 Hz。
設(shè)定海況等級為3級,對應(yīng)標(biāo)準(zhǔn)風(fēng)速為8.23 m/s,得到海浪波譜圖如圖4所示。從圖4可以看出,海浪波譜能量主要分布在0.5~1.5 Hz頻帶范圍和-π/2~π/2角度范圍內(nèi)。
圖4 海況等級為3級海浪方向譜Fig.4 Ocean wave directional spectrum under sea state level 3
對圖4中3級海況下海浪方向譜進(jìn)行等能量四叉樹離散化分解,圖5為網(wǎng)格能量為總能量的0.5%時方向和頻段的離散化結(jié)果。從圖5離散化結(jié)果可以看出,四叉樹等能量分解法在能量密度低的區(qū)域采樣稀疏,在能量密度高的區(qū)域采樣密集。
圖5 基于四叉樹方向譜等能量離散化Fig.5 Equal energy division of directional spectrum based on quadtree
為了客觀比較區(qū)間等分法和四叉樹分解法的仿真速度,對不同的能量百分比條件下區(qū)間等分法和四叉樹分解法的計(jì)算次數(shù)進(jìn)行比較,結(jié)果如表1??梢钥闯瞿芰康确值谋壤叫?,四叉樹分解算法與等間隔分解算法所需要的計(jì)算次數(shù)比例越小。
表1 仿真速度比較Table 1 Simulation speed comparison
波浪傳播主方向?yàn)?5°,磁探儀高度50 m,基于四叉樹分解的優(yōu)化仿真算法仿真不同海況等級下磁探儀靜止條件下采樣磁異常時間歷程信號如圖6所示,并利用Welch功率譜計(jì)算方法進(jìn)行譜分析得到功率譜如圖7所示。從圖7譜分析可知,海浪磁場能量隨海況的增長而迅速的增加。隨著海況的增長,中心頻率是逐漸向低頻方向移動的。這與海浪譜的分布特征是吻合的。
圖6 不同海況下靜止磁探儀采樣海浪磁場信號仿真Fig.6 Simulation on ocean wave generated magnetic field signals sampled by a staying magnetometer under different sea state levels
圖7 不同海況等級靜止磁探儀采樣信號功率譜Fig.7 Power spectrum of ocean wave generated magnetic field signals sampled by a staying magnetometer under different sea state levels
根據(jù)線性海浪理論,磁探儀靜止條件下采樣仿真結(jié)果應(yīng)該能反映線性隨機(jī)海浪磁場外觀上和統(tǒng)計(jì)上的特征,整體統(tǒng)計(jì)特征表現(xiàn)為上下對稱,均值為零,其正態(tài)性偏度和峰度應(yīng)為0和3。對圖6仿真數(shù)據(jù)進(jìn)行統(tǒng)計(jì),得到均值、偏度和峰度如表2所示。根據(jù)表2統(tǒng)計(jì)數(shù)據(jù)可見,海浪磁場仿真結(jié)果與線性隨機(jī)海浪外觀上和統(tǒng)計(jì)上的理論特征吻合。
表2 靜止磁探儀采樣時域統(tǒng)計(jì)檢驗(yàn)Table 2 Time-domain statistical analysis of samples by stationary magnetometer
海況等級為4級條件,海浪傳播主方向60°。磁探儀飛行方向?yàn)?5°,飛行高度50 m,基于四叉樹分解的優(yōu)化仿真算法仿真不同速度下航空磁探儀探測的海浪磁場時間歷程信號如圖8,并利用Welch法得到頻譜如圖9所示。
根據(jù)圖7、9可以看出,隨著飛行速度的增加,磁探儀接收到的海浪磁場信號存在明顯的多普勒效應(yīng)。磁探儀靜止條件下采樣的信號能量主要集中在0~0.4 Hz,而航空磁探儀運(yùn)動速度50 m/s時接收信號的主要頻段集中在0.2~1.6 Hz,存在明顯的頻帶擴(kuò)展和頻率移動。隨著航空磁探儀運(yùn)動速度的增加多普勒效應(yīng)越明顯,這與文獻(xiàn)[9]中實(shí)驗(yàn)分析的結(jié)果是一致的。
圖8 不同飛行速度下運(yùn)動磁探儀采樣海浪磁場信號仿真Fig.8 Simulation on ocean wave generated magnetic field signals sampled by a magnetometer moving with different speeds
圖9 運(yùn)動磁探儀不同飛行速度采樣海浪磁場信號功率譜Fig.9 Power spectrum of ocean wave generated magnetic field signals sampled by a magnetometer moving with different speeds
圖10 運(yùn)動磁探儀不同飛行角度采樣海浪磁場信號功率譜Fig.10 Power spectrum of ocean wave generated magnetic field signals sampled by a magnetometer moving with different angles
海況等級為4級,海浪傳播主方向60°,航空磁探儀運(yùn)動速度80 m/s,基于四叉樹分解的優(yōu)化仿真算法仿真運(yùn)動磁探儀不同運(yùn)動角度采樣信號的功率譜分析結(jié)果如圖10所示。從圖10分析可知,當(dāng)飛行方向與海浪傳播方向接近時,海浪磁場信號頻率更加集中,頻帶范圍更窄,隨著飛向方向與海浪傳播方向夾角的增大,海浪磁場信號帶寬變窄,并且明顯向低頻方向擴(kuò)展。
基于Weaver海浪磁場模型推導(dǎo)了地理坐標(biāo)系下沿任意方向傳播海浪單頻重力波感應(yīng)磁場數(shù)學(xué)模型,在此基礎(chǔ)上基于線性理論推導(dǎo)了實(shí)際海浪磁場數(shù)值模型?;谟^測海浪譜給出了實(shí)際海浪磁場數(shù)值仿真算法,在海浪頻率和海浪方向上更加真實(shí)地描述實(shí)際海浪,并且采用四叉樹理論對仿真算法進(jìn)行了仿真速度優(yōu)化。時頻域驗(yàn)證和分析結(jié)果表明,該仿真算法仿真速度快,仿真結(jié)果與理論和實(shí)際情況吻合。該仿真算法能夠仿真遠(yuǎn)離海岸任意精度的連續(xù)海洋波浪產(chǎn)生的磁場,可用于航空磁探儀海浪磁噪聲背景消除研究,可為進(jìn)一步的海浪磁場實(shí)驗(yàn)研究提供指導(dǎo)。
[1]NELSON J B.Aeromagnetic noise during low-altitude flights over the Scotian shelf[R].Dartmouth:Defence Research &Development,2002:15-18.
[2]LILLEY F E,HITCHMAN A P,MILLIGAN P R,et al.Sea-surface observations of the magnetic signals of ocean swells[J].Geophysical Journal International,2004,159(2):565-572.
[3]李洪平,張海濱.海洋背景磁場模擬計(jì)算及東中國海表層磁場分布[J].海洋技術(shù),2008,27(3):70-74.LI Hongping,ZHANG Haibin.Simulation of ocean background magnetic field and its distribution in the East China Sea[J].Ocean Technology,2008,27(3):70-74.
[4]WEAVER J T.Magnetic variations associated with ocean waves and swell[J].Journal of Geophysical Research,1965,70(8):1921-1929.
[5]OCHADLICK A R.Experimental demonstration of weaver’s model of magnetic fields from ocean waves[R].Washington,DC:Naval Air Development Center,1980:10-12.
[6]閆曉偉,閆輝,肖昌漢.海浪感應(yīng)磁場矢量的模型研究[J].海洋測繪,2011,31(6):8-11.YAN Xiaowei,YAN Hui,XIAO Changhan.Research on model of induced magnetic vector of ocean waves[J].Hydrographic Surveying and Charting,2011,31(6):8-11.
[7]張揚(yáng),康崇,呂金庫.有限深海域的海浪感應(yīng)磁場[J].中國慣性技術(shù)學(xué)報(bào),2012,20(5):594-595.ZHANG Yang,KANG Chong,LYU Jinku.Inductive magnetic field with deep ocean waves[J].Journal of Chinese Inertial Technology,2012,20(5):594-595.
[8]JOCELYN F.Realistic simulation of ocean surface using wave spectra[C]//Proceedings of the First International Conference on Computer Graphics Theory and Applications,Portugal:CCSD Press,2006:76-83.
[9]唐勁飛,龔沈光,林春生.經(jīng)典海浪譜下運(yùn)動飛行器接收到的海浪磁場的噪聲[J].海軍工程大學(xué)學(xué)報(bào),2002,14(6):54-58.TANG Jinfei,GONG Shenguang,LIN CHunsheng.Wave generated magnetic noise received by moving airborne magnetometers under classical spectrum assumptions[J].Journal of Naval University of Engineering,2002,14(6):54-58.
[10]YAAKOBI O,ZILMAN J,MILOH T.Detection of the electromagnetic filed included by the wake of a ship moving in a moderate sea state of finite depth[J].Journal of Engineering Mathematics,2011,70(3):17-27.
[11]SEMKIN S V,SMAGIN V P.The effect of self-induction on magnetic field generated by sea surface waves[J].Atmospheric and Oceanic Physics,2012,48(2):207-213.
[12]STEELE K E,EARLE M D.Directional ocean wave spectra using buoy azimuth,pitch and roll derived from magnetic field components[J].IEEE Journal of Oceanic Engineering,1991,16(4):427-433.
[13]STEELE K E,TENG C C,WANG D W.Wave direction measurements using pitch and roll buoys[J].IEEE Journal of Oceanic Engineering,1992,19(4):349-375.
[14]郭利進(jìn),師五喜,李穎,等.基于四叉樹的自適應(yīng)柵格地圖創(chuàng)建算法[J].控制與決策,2011,26(11):1690-1693.GUO Lijin,SHI Wuxi,LI Ying,et al.Mapping algorithm using adaptive size of occupancy grids based on quad-tree[J].Control and Decision,2011,26(11):1690-1693.
[15]劉揚(yáng),宮阿都,李京.基于數(shù)據(jù)分層分塊的海量三維地形四叉樹簡化模型[J].測繪學(xué)報(bào),2010,39(4):410-413.LIU Yang,GONG Adu,LI Jing.A model for massive 3D terrain simplification based on data block partition and quad-tree[J].Acta Geodaetica et Cartographica Sinica,2010,39(4):410-413.
[16]曾凱,楊華,翟月,等.光電成像干擾圖像質(zhì)量評估[J].電子與信息學(xué)報(bào),2011,33(9):2164-2166.ZENG Kai,YANG Hua,ZHAI Yue,et al.Quality assessment of photoelectric image interference[J].Journal of E-lectronics& Information Technology,2011,33(9):2164-2166.
[17]尹長林,詹慶明,許文強(qiáng),等.大規(guī)模三維地形實(shí)時繪制的簡化技術(shù)研究[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2012,37(5):556-558.YIN Changlin,ZHAN Qingming,XU Wenqiang,et al.Simplification technology for real-time rendering of largescale 3D Terrain[J].Geomatics and Information Science of Wuhan University,2012,37(5):556-558.