劉志學(xué)，宮志華，冷雪冰

(中國白城兵器試驗中心, 吉林 白城 137001)

基于相位差分的彈丸進(jìn)動周期提取

劉志學(xué)，宮志華，冷雪冰

(中國白城兵器試驗中心, 吉林 白城 137001)

進(jìn)動周期是表征炮彈、火箭彈等常規(guī)武器無控彈丸彈道特性的重要指標(biāo)。為有效獲取彈丸進(jìn)動周期，提出了基于相位差分的彈丸進(jìn)動周期提取方法。彈丸進(jìn)動會對連續(xù)波雷達(dá)回波產(chǎn)生微多普勒調(diào)制，建模分析表明，進(jìn)動彈丸的微動速度遵循正弦規(guī)律變化，其相位變化率與進(jìn)動周期呈反比關(guān)系。利用Hilbert變換可以提取微動速度的瞬時相位，由相位差分即可獲得彈丸進(jìn)動周期。仿真實(shí)驗驗證了該方法對信噪比大于7 dB的微動速度提取進(jìn)動周期具有更好的性能。利用該方法對實(shí)際測量數(shù)據(jù)進(jìn)行處理，獲取了彈丸準(zhǔn)確、平滑的進(jìn)動周期數(shù)據(jù)。

兵器科學(xué)與技術(shù)；連續(xù)波雷達(dá)；進(jìn)動周期；微多普勒；Hilbert變換

0 引言

炮彈、火箭彈等常規(guī)武器無控彈丸(以下簡稱彈丸)發(fā)射時，擾動因素使彈軸偏離預(yù)定速度矢量方向，影響射擊精度和飛行穩(wěn)定。高速自旋運(yùn)動形成一個使彈丸圍繞質(zhì)心轉(zhuǎn)動的力矩，使彈軸圍繞速度矢量方向旋轉(zhuǎn)，彈道學(xué)上稱之為進(jìn)動[1]。進(jìn)動特性同彈丸的彈道特性、電磁散射特性一樣, 是彈丸的固有屬性, 可以作為彈丸類目標(biāo)識別的重要特征量，也可以用來驗證彈丸的設(shè)計以及分析落點(diǎn)散布。

相對于彈丸的飛行來講，進(jìn)動屬于微動[2-3]，對雷達(dá)回波產(chǎn)生幅度和相位調(diào)制。其中幅度調(diào)制使目標(biāo)雷達(dá)散射截面(RCS)序列呈現(xiàn)周期性變化，基于RCS序列的目標(biāo)進(jìn)動周期提取方法研究開展的較早，美國開展的Firefly飛行試驗專門研究進(jìn)動特性[4]，但技術(shù)細(xì)節(jié)未披露?，F(xiàn)有關(guān)于提取RCS序列周期的方法包括周期圖法[5-6]、自相關(guān)函數(shù)法[7]、平均幅度差函數(shù)法[8]、循環(huán)幅度差函數(shù)法、循環(huán)自相關(guān)函數(shù)法，以及循環(huán)幅度差結(jié)合自相關(guān)函數(shù)法[9-10]等，文獻(xiàn)[11]利用三角函數(shù)來擬合RCS序列，再求得使擬合誤差最小的進(jìn)動頻率, 文獻(xiàn)[12]對經(jīng)過平動和微動調(diào)制后的彈丸RCS數(shù)據(jù)進(jìn)行分析后,利用希爾伯特- 黃變換(HHT)時頻分析算法提取出了進(jìn)動周期。總的來講，基于RCS序列的目標(biāo)進(jìn)動周期提取利用的是彈丸進(jìn)動對雷達(dá)回波的幅度調(diào)制信息，但該類方法對雷達(dá)回波信噪比、RCS序列循環(huán)平穩(wěn)性和觀察時間要求較高。

進(jìn)動對雷達(dá)回波的相位調(diào)制產(chǎn)生微多普勒信號，并疊加到多普勒信號上[13]。微多普勒調(diào)制信號反映了微動的瞬時相位特性，其頻率表征了目標(biāo)瞬時微動速度。微多普勒目前主要用于彈道導(dǎo)彈、旋翼直升飛機(jī)的分析識別，以及人體手、腳的擺動和橋梁的振動等探測，而關(guān)于彈道目標(biāo)微動特性研究的技術(shù)文獻(xiàn)較少[13-14]。文獻(xiàn)[15]研究了基于經(jīng)驗?zāi)B(tài)分解算法的多分量正弦調(diào)頻信號分離方法,并使用短時傅里葉變換得到了每個本征模態(tài)函數(shù)的瞬時頻率；文獻(xiàn)[16]中提出了利用回波相關(guān)方法和頻譜展寬估計進(jìn)動參數(shù)的方法；文獻(xiàn)[17]利用時域滑窗自相關(guān)處理方法提取了彈丸進(jìn)動周期。上述方法本質(zhì)上是利用微多普勒信號周期間的相關(guān)性來估計周期，短時傅里葉變換對數(shù)據(jù)量有較高要求，自相關(guān)函數(shù)類方法的誤差較大。

本文著重分析彈丸進(jìn)動對連續(xù)波雷達(dá)回波產(chǎn)生的微多普勒調(diào)制模型，針對進(jìn)動微動速度的特點(diǎn),提出了通過Hilbert變換求取微動速度相位，利用相位差分估計彈丸進(jìn)動周期的數(shù)據(jù)處理方法,可以有效提取彈丸進(jìn)動周期。

1 進(jìn)動調(diào)制模型

連續(xù)波雷達(dá)在對彈丸采用尾追式跟蹤測量時，彈丸進(jìn)動示意圖如圖1所示。定義進(jìn)動坐標(biāo)系[1]：以彈丸質(zhì)心為坐標(biāo)系原點(diǎn)o；以彈丸質(zhì)心速度方向線為ox軸，沿質(zhì)心運(yùn)動方向為正；oy軸垂直于ox軸，向上為正；oz軸按右手法則確定，垂直于oxy平面，向右為正。這樣定義的進(jìn)動坐標(biāo)系在外彈道學(xué)中亦稱為彈道坐標(biāo)系。

圖1 進(jìn)動示意圖Fig.1 Schematic diagram of precession

設(shè)彈丸自旋角速度為ω，彈軸進(jìn)動角速度為Ω，進(jìn)動時彈軸與速度方向ox軸夾角為θ. 雷達(dá)發(fā)射電磁波照射彈丸，彈上主要散射點(diǎn)P位于彈丸底部，P點(diǎn)與彈丸質(zhì)心距離為l，則散射點(diǎn)P在彈道坐標(biāo)系中的坐標(biāo)如(1)式所示：

(1)

再建立雷達(dá)測量坐標(biāo)系，如圖2所示。OX軸沿水平線平行射擊面，射擊方向為正；OY軸垂直于OX軸，向上為正；OZ軸按右手法則確定，垂直于OXY平面，向右為正。

圖2 雷達(dá)測量坐標(biāo)系與彈道坐標(biāo)系幾何關(guān)系Fig.2 Geometric relationship of radar coordinate system and trajectory coordinate system

將雷達(dá)測量坐標(biāo)系平移至彈丸質(zhì)心o形成基準(zhǔn)坐標(biāo)系oXNYNZN.ox軸與oXN軸夾角β為速度高低角，oz軸與oZN軸夾角ψ為速度方向角。

彈上主要散射點(diǎn)P在基準(zhǔn)坐標(biāo)系中的坐標(biāo)值可由彈道坐標(biāo)系經(jīng)兩次旋轉(zhuǎn)而成，如(2)式所示：

(2)

通常情況下彈丸的速度方向角ψ≈0，由此可得散射點(diǎn)P在任意時刻到雷達(dá)的瞬時斜距rt，如(3)式所示：

(3)

式中：a=-lcosθ(cosβcosEcosA+sinβsinE)；b=lsinθcosEsinA；c=lsinθ(sinβcosEcosA-cosβsinE)；Rt為雷達(dá)到彈丸質(zhì)心的斜距；A為雷達(dá)測量彈丸方位角；E為雷達(dá)測量彈丸俯仰角。

在一定的觀察時間內(nèi)，方位角A和俯仰角E近似為常數(shù)，因此(3)式對時間求導(dǎo)數(shù)得到彈丸主要散射點(diǎn)的徑向速度為

(4)

式中：

vd(t)為彈丸質(zhì)心的徑向速度；vm(t)為彈丸進(jìn)動引起的微動速度(簡稱微動速度)；vr(t)為彈丸的徑向速度。當(dāng)目標(biāo)遠(yuǎn)離雷達(dá)時Rt?l，這時g≈1,q≈1,則彈丸的徑向速度可以化簡為

vr(t)≈vd(t)+vm(t),

(5)

而微動速度可以表示為

vm(t)=-bΩsin (Ωt)+cΩcos (Ωt)=
dΩcos (Ωt+φ),

(6)

通過上述分析可知，彈丸的徑向速度是質(zhì)心徑向速度與微動速度之和，且微動速度遵循正弦規(guī)律變化,其相位與進(jìn)動周期的關(guān)系為

(7)

式中：Ф(t) 為微動速度相位;Tm為進(jìn)動周期;fm為進(jìn)動頻率。

2 進(jìn)動周期提取方法

彈丸的進(jìn)動頻率一般很低，直接采用頻域處理的方法提取進(jìn)動頻率，其分辨力難以滿足要求。采用時域滑窗自相關(guān)處理方法提取的進(jìn)動周期數(shù)據(jù)率低，受噪聲影響較大。從現(xiàn)實(shí)物理意義上來說，進(jìn)動頻率與微動速度的相位在時域上存在積分匹配的關(guān)系，即通過對微動速度的相位的微分可以得到微動頻率。

彈丸的進(jìn)動周期為非線性變化的，所以微動速度為非平穩(wěn)信號，其頻率和周期的定義方式是以Hilbert變換為基礎(chǔ)的。對微動速度進(jìn)行Hilbert變換,求出解析信號再對其相位微分,可以得到頻率，而周期與頻率互為倒數(shù)。由(6)式可知，微動速度遵循正弦規(guī)律變化，具有窄帶特性，可以應(yīng)用Hilbert變換提取瞬時參數(shù)。為此本文提出了基于Hilbert變換相位差分的進(jìn)動周期提取方法，提取流程如圖3所示。

圖3 進(jìn)動周期提取流程Fig.3 Flow chart of precession period extraction

微動速度vm(t)的Hilbert變換為

(8)

式中：vmi(k)為vm(k)的Hilbert變換；τ是卷積中間變量；vmc(k)為vm(k)的復(fù)數(shù)解析形式，vm(k)是vmc(k)的實(shí)部，vmi(k) 是vmc(k)的虛部。vmc(k)的相位為

(9)

由Φ1(k)估計進(jìn)動周期，需要對其進(jìn)行數(shù)值微分,利用向前差商法可以實(shí)現(xiàn)

Tm(k)=2πTs/[Φ1(k+1)-Φ1(k)],

(10)

式中：Ts為徑向速度的采樣間隔。

從雷達(dá)測量數(shù)據(jù)提取飛行彈丸進(jìn)動周期，由于彈丸質(zhì)心的徑向速度vr(t)和微動速度vm(t)存在固有的線性耦合現(xiàn)象，對雷達(dá)錄取的目標(biāo)徑向速度進(jìn)行自由節(jié)點(diǎn)樣條函數(shù)擬合[18]獲得彈丸質(zhì)心的徑向速度，提取微動速度數(shù)據(jù)vm(t)；對微動速度數(shù)據(jù)vm(t)進(jìn)行Hilbert變換獲得微動速度的解析形式數(shù)據(jù)vmc(t)，求取解析微動速度vmc(t)的相位Ф(t)；對Ф(t)進(jìn)行自由節(jié)點(diǎn)樣條函數(shù)擬合消除測量噪聲的影響；利用(10)式求取進(jìn)動頻率和進(jìn)動周期。

3 仿真及試驗分析

連續(xù)波雷達(dá)測量目標(biāo)徑向速度時不可避免地 引入測量誤差，此時雷達(dá)輸出的彈丸徑向速度由(5)式擴(kuò)展為

vr(t)=vd(t)+vm(t)+n(t),

(11)

式中：n(t)為連續(xù)波雷達(dá)速度測量誤差，服從～[0,σ2]正態(tài)分布，σ2為連續(xù)波雷達(dá)速度測量誤差功率。在考核進(jìn)動周期提取方法的性能時必須考慮速度測量誤差的影響，為此(6)式和(11)式可得微動速度功率與速度測量誤差功率之比，即微動速度信噪比SNR，

(12)

3.1 仿真實(shí)驗

為了驗證文中提出的進(jìn)動周期提取方法的性能，作如下仿真：時間t=0～25 s，采樣間隔Ts=0.01 s，速度測量誤差服從～[0, 0.1]正態(tài)分布，質(zhì)心徑向速度vd=(0.4t2-15t+660) m/s，微動速度vm=Amcos(2πt/Tm)，其中，Am為微動速度振幅，進(jìn)動周期Tm=1/(0.5+0.02t) s. 信噪比SNR從3 dB到14 dB變化，作200次Monte Carlo仿真實(shí)驗，得到的仿真結(jié)果如圖4所示。從仿真結(jié)果可以看出：當(dāng)信噪比大于7 dB時，相位差分法的性能優(yōu)于自相關(guān)法；當(dāng)信噪比小于7 dB時，由于基于Hilbert變換的相位計算對噪聲比較敏感，利用相位差分法提取的進(jìn)動周期會發(fā)生錯誤，進(jìn)動周期估計誤差顯著提升，相位差分法不再可用。

圖4 相位差分法性能分析Fig.4 Analysis of the capability of phase differential method

在仿真實(shí)驗中，信噪比為8 dB時的進(jìn)動周期處理結(jié)果如圖5所示。從仿真結(jié)果可以看出，相位差分法得到的進(jìn)動周期更平滑和準(zhǔn)確。由于Hilbert變換存在端點(diǎn)效應(yīng)，進(jìn)動周期數(shù)據(jù)的起始和結(jié)束部分存在顯著異常，在實(shí)際應(yīng)用中應(yīng)該舍棄異常部分或進(jìn)行進(jìn)一步的處理。

圖5 進(jìn)動周期與估計誤差Fig.5 Precession period and estimated error

圖6 實(shí)測速度Fig.6 Real radial velocity

3.2 應(yīng)用分析

受雷達(dá)布站、彈道特性等影響，只有在彈道下降段，連續(xù)波雷達(dá)才能錄取到進(jìn)動微多普勒信號。圖6為連續(xù)波雷達(dá)測量得到的某型火箭彈下降段徑向速度數(shù)據(jù)，從數(shù)據(jù)曲線中可以明顯看到由彈丸進(jìn)動引起的準(zhǔn)周期變化。

對速度數(shù)據(jù)進(jìn)行自由節(jié)點(diǎn)樣條函數(shù)擬合，求取質(zhì)心徑向速度和微動速度，如圖7所示，其中微動速度包含雷達(dá)測量噪聲。

圖7 質(zhì)心徑向速度和微動速度Fig.7 Radial velocity of center of mass and radial precession velocity

圖8 微動速度相位Fig.8 Phase of radial precession velocity

對微動速度做Hilbert變換得到微動速度的解析形式，計算解析微動速度的相位，如圖8中所示。微動速度的相位是連續(xù)變化的，但是直接求得的相位在±π rad之間，當(dāng)相位跨越±π rad時存在從+π rad到-π rad的跳點(diǎn)，即存在纏繞問題。

對微動速度相位做解纏繞處理，即對從+π rad到-π rad的跳點(diǎn)加上±2π rad的倍數(shù)來解除纏繞現(xiàn)象，使微動速度的相位連續(xù)變化。由于測量誤差對相位的影響較大，需要對微動速度的相位進(jìn)行自由節(jié)點(diǎn)樣條函數(shù)擬合來消除。解纏繞后的微動速度的相位及其擬合結(jié)果如圖9所示。

圖9 解纏繞后的微動速度相位Fig.9 Unwraped radial precession velocity phase

根據(jù)(8)式和(9)式對解纏繞和擬合后的微動速度相位處理，求取進(jìn)動頻率和進(jìn)動周期。采用本文中的Hilbert變換方法處理的某型彈丸的進(jìn)動周期如圖10中紅色曲線所示，圖中藍(lán)色曲線為采用文獻(xiàn)[16]中自相關(guān)法處理得到的進(jìn)動周期。兩種處理方法得到進(jìn)動周期的趨勢非常吻合，證明了該方法的正確性和有效性。

圖10 進(jìn)動周期Fig.10 Precession period

4 結(jié)論

彈丸在運(yùn)動過程中的進(jìn)動現(xiàn)象對雷達(dá)回波產(chǎn)生了微多普勒頻率調(diào)制。連續(xù)波雷達(dá)測量得到的徑向速度數(shù)據(jù)中包含質(zhì)心徑向速度和由進(jìn)動引起的微動速度。對徑向速度時間序列進(jìn)行數(shù)據(jù)處理，利用自由節(jié)點(diǎn)樣條函數(shù)擬合可以得到較為準(zhǔn)確的質(zhì)心徑向速度分量，并由此獲得微動速度，對微動速度進(jìn)行Hilbert變換處理可以得到微動速度的相位信息，對相位做差分可以提取出彈丸的進(jìn)動周期。仿真實(shí)驗結(jié)果表明，在微動速度與速度測量誤差的信噪比大于7 dB時，相位差分法更適于實(shí)際工程應(yīng)用。實(shí)際試驗結(jié)果表明，采用相位差分法提取的進(jìn)動周期數(shù)據(jù)準(zhǔn)確、平滑。

References)

[1] 韓子鵬. 彈箭外彈道學(xué)[M]. 北京: 北京理工大學(xué)出版社, 2008: 153-162. HAN Zi-peng. Projectile and rocket exterior ballistics [M]. Beijing: Beijing University of Aeronautics and Astronauts Press, 2008: 153-162. (in Chinese)

[2] 陳行勇. 微動目標(biāo)雷達(dá)特征提取技術(shù)研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2005. CHEN Hang-yong. Research on radar signature extraction from target with micro-motions[D]. Changsha: National University of Defense Technology, 2005. (in Chinese)

[3] 陳行勇, 黎湘, 郭桂蓉, 等. 微進(jìn)動彈道導(dǎo)彈目標(biāo)雷達(dá)特征提取[J]. 電子與信息學(xué)報，2006，28(4): 643-646. CHEN Hang-yong, LI Xiang, GUO Gui-rong, et al. Radar feature extraction of micro-precession ballistic missile warhead[J]. Journal of Electronics and Information Technology, 2006, 28(4): 643-646. (in Chinese)

[4] Schultz K, Davidson S, Stein A, et al. Range Doppler laser radar for midcourse discrimination: the firefly experiments[C]∥Proceedings of the 2nd Annual AIAA SDIO Interceptor Technology Conference. Albuquerque, New Mexico: AIAA, 1993: 2653-2664.

[5] Liu L H, Wang Z , Hu W D. Procession period extraction of ballistic missile based on radar measurement[C]∥CIE International Conference on Radar. Shanghai, China: Chinese Institute of Electronics,2006.

[6] 張淵, 張愛成. 利用FFT方法分析自旋目標(biāo)RCS變化周期[J].飛行器測控學(xué)報, 2012, 31(1): 26-29. ZHANG Yuan, ZHANG Ai-cheng. Analysis of the RCS variation cycle of spin-stabilized objects using FFT method[J]. Journal of Spacecraft TT & C Technology, 2012, 31(1): 26-29. (in Chinese)

[7] 陳翱. 基于RCS 序列的彈道中段目標(biāo)微動提取技術(shù)[J]. 現(xiàn)代雷達(dá), 2012, 34(6): 78-81. CHEN Ao. A micromoton extraction technique for target in intermediate trajectory based on RCS sequence[J]. Modern Radar, 2012, 34(6): 78-81. (in Chinese)

[8] Ross M, Shafer H. Average magnitude difference function pitch extractor[J]. IEEE Transactions on Acoustic, Speech and Signal Processing, 1974, 22(5): 355-362.

[9] 馮德軍, 丹梅, 馬梁. 一種魯棒的彈道目標(biāo)RCS周期估計方法[J]. 航天電子對抗, 2008, 24(2): 5-8. FENG De-jun, DAN Mei, MA Liang. A robust RCS periodicity estimation algorithm for ballistic target[J]. Aerospace Electronic Warfare, 2008, 24(2): 5-8. (in Chinese)

[10] 康猛, 王春花, 郝明, 等. 彈道目標(biāo)進(jìn)動周期特征提取研究[J]. 現(xiàn)代雷達(dá), 2010, 32(11): 29-32. KANG Meng, WANG Chun-hua, HAO Ming, et al. A study of procession-period extraction method of ballistic target[J]. Modern Radar, 2010, 32(11): 29-32. (in Chinese)

[11] 張仕元. 基于三角函數(shù)擬合和的RCS序列進(jìn)動周期估計[J].電子與信息學(xué)報, 2014, 36(6):1389-1393. ZHANG Shi-yuan. Procession period estimation of RCS sequences based on trigonometric function fitting [J]. Journal of Electronics & Information Technology,2014, 36(6):1389-1393. (in Chinese)

[12] 孟路穩(wěn), 周沫, 察豪,等. 彈道導(dǎo)彈目標(biāo)回波模擬與微動特征提取[J].火力與指揮控制, 2015, 39(10): 1712-1715. MENG Lu-wen, ZHOU Mo, CHA Hao, et al. Echo simulation and micro-motion feature extraction of ballistic missile targets[J]. Fire Control & Command Control, 2015, 39(10): 1712-1715. (in Chinese)

[13] Chen V C. Micro-Doppler effect of micro-motion dynamics: a review[J]. Proceedings of SPIE on Independent Component Analyses, Wavelets, and Neural Networks, 2003, 5102:240-249.

[14] Chen V C, Li F, Ho S, et al. Micro-Doppler effect in radar-phenomenon, model and simulation study [J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(1):2-21.

[15] 牛杰,劉永祥,秦玉亮,等. 一種基于經(jīng)驗?zāi)B(tài)分解的錐體目標(biāo)雷達(dá)微動特征提取新方法[J].電子學(xué)報, 2011, 39(7): 1712-1715. NIU Jie, LIU Yong-xiang, QIN Yu-liang, et al. A new method of radar micro-motion feature extraction of cone target based on empirical mode decomposition[J]. Acta Electronica Sinica, 2011, 39(7): 1712-1715. (in Chinese)

[16] 李康樂, 姜衛(wèi)東, 黎湘. 彈道目標(biāo)微動特征分析與提取方法[J]. 系統(tǒng)工程與電子技術(shù), 2010, 32(1): 115-118. LI Kang-le, JIANG Wei-dong, LI Xiang. Micro-motion feature analysis and extraction methods for ballistic target[J].Journal of Systems Engineering and Electronics, 2010, 32(1): 115-118. (in Chinese)

[17] 李益民,蘇東林, 冷雪冰. 高速旋轉(zhuǎn)彈丸進(jìn)動周期提取[J]. 北京航空航天大學(xué)學(xué)報, 2010, 36(11): 1335-1338. LI Yi-min，SU Dong-lin, LENG Xue-bing．Precession period extraction of high-rate rotating projectile[J]．Journal of Beijing University of Aeronautics and Astronautics，2010，36(11):1335-1338. (in Chinese)

[18] 郭文勝，宮志華，董立濤，等．雷達(dá)測距擬合微分求速方法研究[J]．現(xiàn)代雷達(dá)，2010，32(8)：33-38． GUO Wen-sheng，GONG Zhi-hua，DONG Li-tao，et al．A study on getting radar radial velocity by range fitting and differential method[J]. Modern Radar, 2010, 32(8): 33-38. (in Chinese)

Extraction of Precession Period of Projectile Based on Phase Differential

LIU Zhi-xue, GONG Zhi-hua, LENG Xue-bing

(Baicheng Ordnance Test Center of China, Baicheng 137001, Jilin, China)

For uncontrolled projectile, such as cannonballs, rockets, etc, precession period is an important parameter of the trajectory characteristic. To obtain the precession period of projectile effectively, a precession period extraction method based on phase differential is introduced. Precession of projectile brings micro-Doppler modulation on echoes of continuous-wave radar, and its mathematic model indicates that the micro-velocity produced by precession changes with the law of sine, and the change rate of its phase is inversely proportional to precession period. The phase of the micro-velocity can be obtained by Hilbert transform, and then the precession period can be derived by the phase differential of the micro-velocity. The experimental results show that, for the micro-velocity signal which S/N is greater than 7 dB, precession period can be effectively extracted with the proposed method, and this method is more precise and smooth for actual measuring data.

ordnance science and technology; continuous wave radar; precession period; micro-Doppler; Hilbert transform

2016-04-18

劉志學(xué)(1972—),男,高級工程師。E-mail：liuzx09@qq.com

TJ012.3+6

A

1000-1093(2017)01-0177-07

10.3969/j.issn.1000-1093.2017.01.023