基于時(shí)域縮放ERA方法的工作模態(tài)辨識(shí)技術(shù)

王 亮，張 妍， 蔡毅鵬，南宮自軍，陳國平

(1.中國運(yùn)載火箭技術(shù)研究院，北京 100076； 2.南京航空航天大學(xué) 振動(dòng)工程研究所，南京 210016)

【信息科學(xué)與控制工程】

基于時(shí)域縮放ERA方法的工作模態(tài)辨識(shí)技術(shù)

王 亮1，張 妍1， 蔡毅鵬1，南宮自軍1，陳國平2

(1.中國運(yùn)載火箭技術(shù)研究院，北京 100076； 2.南京航空航天大學(xué) 振動(dòng)工程研究所，南京 210016)

研究了基于時(shí)域縮放選帶，使用振動(dòng)模態(tài)特征系統(tǒng)實(shí)現(xiàn)算法(eigensystem realization algorithm，ERA)的工作模態(tài)辨識(shí)技術(shù)。分析了ERA環(huán)境激勵(lì)模態(tài)辨識(shí)方法的理論，給出了時(shí)域縮放選帶技術(shù)和真實(shí)模態(tài)挑選技術(shù)的方法以及工作流程。通過算例研究了采用時(shí)域縮放技術(shù)的ERA工作模態(tài)辨識(shí)方法。研究結(jié)果表明：使用該方法可以有效地辨識(shí)出結(jié)構(gòu)的模態(tài)，且可以采用滑動(dòng)時(shí)間窗的方法獲得模態(tài)隨時(shí)間的變化結(jié)果；用于辨識(shí)的數(shù)據(jù)長(zhǎng)度需要包含所需識(shí)別模態(tài)的6～8周期；阻尼辨識(shí)結(jié)果的偏差隨阻尼比理論值變大而變小，對(duì)于理論值從1%～5%，辨識(shí)偏差從175%減小至10%；一階和二階阻尼比辨識(shí)結(jié)果偏差基本一致；振型辨識(shí)結(jié)果MAC值均在0.99以上。

模態(tài)辨識(shí)；時(shí)域縮放；ERA；工作模態(tài)

由于飛行器在飛行過程中的激勵(lì)無法精確測(cè)量，另外受到氣動(dòng)加熱、外部氣動(dòng)力、飛行過載、飛行的剛體運(yùn)動(dòng)和發(fā)動(dòng)機(jī)推力等外部因素的影響，不僅影響飛行器的模態(tài)頻率，還影響飛行器的模態(tài)阻尼特性。設(shè)計(jì)結(jié)果與飛行結(jié)果的差異，使得控制能力挖潛、設(shè)計(jì)裕度的準(zhǔn)確把握和研制難度的降低等方面無法深入開展，難以滿足未來戰(zhàn)術(shù)飛行器研制的發(fā)展需求。因此，基于遙測(cè)數(shù)據(jù)，采用環(huán)境激勵(lì)模態(tài)辨識(shí)技術(shù)對(duì)飛行過程中的結(jié)構(gòu)模態(tài)參數(shù)進(jìn)行辨識(shí)，對(duì)設(shè)計(jì)參數(shù)進(jìn)行驗(yàn)證顯得非常重要。

環(huán)境激勵(lì)模態(tài)辨識(shí)技術(shù)是將互相關(guān)函數(shù)同傳統(tǒng)時(shí)域模態(tài)分析法相結(jié)合的方法，將響應(yīng)之間的互相關(guān)函數(shù)代替?zhèn)鹘y(tǒng)時(shí)域模態(tài)分析法中的自由振動(dòng)響應(yīng)或脈沖響應(yīng)函數(shù)。該方法不需要測(cè)量激勵(lì)型號(hào)，而僅依靠各通道的時(shí)域響應(yīng)數(shù)據(jù)進(jìn)行系統(tǒng)的模態(tài)辨識(shí)。環(huán)境激勵(lì)模態(tài)辨識(shí)的理論和思想在20世紀(jì)70年代初期出現(xiàn)，經(jīng)過幾十年的發(fā)展，形成了多種模態(tài)辨識(shí)的方法[1-2]。

工作模態(tài)分析研究最早可以追溯到1968年cole的單階模態(tài)測(cè)試的隨機(jī)減量法[3]。1973年Ibrahim提出了一種參數(shù)識(shí)別的方法，該方法僅利用時(shí)域信號(hào)即可進(jìn)行識(shí)別工作，經(jīng)多年的不斷完善形成了獨(dú)具一格的Ibrahim時(shí)域法(ITD法)。1983年Metgeay提出了單參考點(diǎn)復(fù)指數(shù)法[4]，其核心是最小二乘估計(jì)，后來Leuridan和vold進(jìn)一步發(fā)展了多參考點(diǎn)復(fù)指數(shù)法(PRCE)。1984年P(guān)appa發(fā)表了特征系數(shù)實(shí)現(xiàn)法的專著，該方法利用了線性系統(tǒng)的狀態(tài)方程和系統(tǒng)最小實(shí)現(xiàn)理論，屬于多輸入多輸出的模態(tài)參數(shù)辨識(shí)方法，通過構(gòu)造Hankel矩陣，利用奇異值分解技術(shù)，確定相互描述狀態(tài)方程的系統(tǒng)矩陣和輸入、輸出矩陣，構(gòu)成最小階的系統(tǒng)，通過求解，得到系統(tǒng)的模態(tài)參數(shù)。

90年代以后，隨著測(cè)試技術(shù)、信號(hào)分析技術(shù)和計(jì)算機(jī)技術(shù)的發(fā)展，模態(tài)參數(shù)辨識(shí)理論研究及應(yīng)用獲得很大進(jìn)展。美國SADIA國家實(shí)驗(yàn)室的James和Carne[5]在1995年證明了系統(tǒng)脈沖響應(yīng)與白噪聲激勵(lì)時(shí)兩點(diǎn)之間的響應(yīng)互相關(guān)函數(shù)有相似的表達(dá)式，從而將運(yùn)用脈沖響應(yīng)函數(shù)進(jìn)行參數(shù)辨識(shí)擴(kuò)展到應(yīng)用相關(guān)函數(shù)進(jìn)行參數(shù)識(shí)別，這就是自然激勵(lì)技術(shù)(NExT)方法，并將該方法用于高速汽輪機(jī)葉片在工作狀態(tài)下固有頻率和阻尼比的識(shí)別。

近年來，國內(nèi)外在模態(tài)參數(shù)識(shí)別領(lǐng)域研究在許多理論和應(yīng)用中取得了較多的成果[6-8]。其中練繼建等[9]對(duì)基于熵降噪的水工結(jié)構(gòu)振動(dòng)模態(tài)特征系統(tǒng)實(shí)現(xiàn)算法進(jìn)行了研究。劉興漢等[10]對(duì)改進(jìn)的隨機(jī)子空間法進(jìn)行了研究。于開平等[11]用小波分析方法對(duì)結(jié)構(gòu)系統(tǒng)的脈沖響應(yīng)函數(shù)進(jìn)行小波變換利用小波變換的幅值、相位與阻尼比、頻率的關(guān)系進(jìn)行參數(shù)辨識(shí)。王彤等[12]提出了一種基于頻域空間域分解的工作模態(tài)分析方法。黃琴等[13]提出并實(shí)現(xiàn)了一種基于隨機(jī)減量技術(shù)和復(fù)模態(tài)指示因子函數(shù)法的新型頻域運(yùn)行狀態(tài)模態(tài)參數(shù)識(shí)別方法。

綜上所述，針對(duì)工作模態(tài)辨識(shí)中結(jié)構(gòu)模態(tài)阻尼比辨識(shí)不穩(wěn)定的問題，本文在不改變ERA辨識(shí)算法的核心思想的基礎(chǔ)上，提出了算法的前后處理方法，增加已有算法的阻尼比辨識(shí)穩(wěn)定性。其中以已知特性的懸臂梁為例，提出了基于時(shí)域縮放選帶，使用振動(dòng)模態(tài)特征系統(tǒng)實(shí)現(xiàn)算法的工作模態(tài)辨識(shí)技術(shù)，研究了時(shí)域縮放選帶技術(shù)和真實(shí)模態(tài)挑選技術(shù)的方法以及工作流程。

1 模態(tài)辨識(shí)技術(shù)

ERA法屬于一種多輸入多輸出的時(shí)域整體模態(tài)參數(shù)辨識(shí)方法[14]。對(duì)于n維線性系統(tǒng)，當(dāng)振動(dòng)系統(tǒng)響應(yīng)采用位移或速度傳感器測(cè)量時(shí)，振動(dòng)方程用向量表示為

(1)

其中Α、B、G別為系統(tǒng)矩陣、控制矩陣和觀測(cè)矩陣。假設(shè)激勵(lì)點(diǎn)數(shù)是L，測(cè)量響應(yīng)點(diǎn)數(shù)是M，則B的階數(shù)是(2n×L)階，激勵(lì)列陣F(t)的階數(shù)為L(zhǎng)，輸出向量Z為M階列陣。

設(shè)離散時(shí)間點(diǎn)為k=0,1,2,…,采樣時(shí)間間隔為Δt，則t=t0+kΔt，簡(jiǎn)記kΔt為k，因此，連續(xù)系統(tǒng)以加速度表示輸出的離散形式為

(2)

(3)

Z變換后得到傳遞函數(shù)

H(z)=z-2GA2(I-z-1A1)-1B1+z-1GB2+GB

(4)

構(gòu)造Hankel矩陣

(5)

整理得

(6)

其中

α、β分別為能觀、能控指數(shù)。

在式(6)中，令k=1，對(duì)H(0)做奇異值分解H(0)=U∑VT，可推導(dǎo)出

(7)

設(shè)系統(tǒng)矩陣A的特征值矩陣為Λ，特征矢量矩陣為ψ′，由指數(shù)矩陣的性質(zhì)而知A1的特征矢量與A的相同，A1的特征值矩陣為

Z=eΛΔ t=diag(z1,z2,…,z2n)

(8)

由此可確定各個(gè)模態(tài)振動(dòng)的固有頻率、阻尼比和模態(tài)矩陣。

固有頻率為

阻尼比為

模態(tài)矩陣為

Φ=Gψ

2 時(shí)域縮放選帶技術(shù)

本文針對(duì)工作模態(tài)辨識(shí)中結(jié)構(gòu)模態(tài)阻尼比辨識(shí)不穩(wěn)定的問題，在不改變ERA辨識(shí)算法核心思想的基礎(chǔ)上，提出了算法的前后處理方法，增加已有算法的阻尼比辨識(shí)穩(wěn)定性。

時(shí)域縮放選帶技術(shù)：研究時(shí)域縮放技術(shù)與ERA相結(jié)合，在ERA方法前后增加前后處理方法，前者采用小頻帶濾波方法縮小關(guān)心頻帶，達(dá)到放大有效信息的作用，去除無用頻帶的信息，后者在辨識(shí)結(jié)果中采用統(tǒng)計(jì)平均的方法進(jìn)行模態(tài)頻率和阻尼比數(shù)據(jù)處理，可提高有測(cè)試噪聲影響情況下的結(jié)果魯棒性。

實(shí)施流程圖如圖1所示。具體操作時(shí)，首先通過滑動(dòng)矩形窗滑動(dòng)時(shí)域信號(hào)，接著根據(jù)矩形窗內(nèi)信號(hào)的功率譜密度分析結(jié)果，選擇待辨識(shí)模態(tài)的頻帶范圍，再對(duì)原始各通道信號(hào)進(jìn)行選帶濾波和重采樣，再根據(jù)信號(hào)的數(shù)據(jù)長(zhǎng)度自動(dòng)生成ERA方法的階次，進(jìn)行ERA方法的模態(tài)辨識(shí)，最后將辨識(shí)結(jié)果與功率譜密度曲線諧振峰峰值進(jìn)行對(duì)比，去除阻尼比過大的模態(tài)后，篩選出較為真實(shí)的模態(tài)。

3 算例

3.1 梁模型介紹

算例計(jì)算時(shí)，對(duì)懸臂梁懸臂端進(jìn)行白噪聲激勵(lì)，將各處加速度響應(yīng)作為輸入進(jìn)行工作模態(tài)辨識(shí)。懸臂梁模型如圖2所示，各參數(shù)如表1所示。

表1 懸臂梁參數(shù)

圖1 加入時(shí)域縮放技術(shù)和有效模態(tài)篩選技術(shù)的ERA工作模態(tài)辨識(shí)技術(shù)流程

圖2 懸臂梁模型

采用wilson-θ方法用800 Hz頻帶白噪聲進(jìn)行激勵(lì)懸臂梁端部，獲得各位置加速度響應(yīng)，典型位置如懸臂梁端部加速度時(shí)域響應(yīng)，其功率譜密度曲線如圖3所示。

圖3 懸臂梁端部響應(yīng)的功率譜密度曲線

從圖3可以看出，在400 Hz頻帶內(nèi)，存在兩階模態(tài)，分別位于20 Hz和120 Hz，從200 Hz向上，響應(yīng)信噪比較差，因此以下針對(duì)這兩階模態(tài)進(jìn)行模態(tài)辨識(shí)研究。

3.2 辨識(shí)用數(shù)據(jù)長(zhǎng)度研究

為了研究可有效辨識(shí)出模態(tài)的數(shù)據(jù)長(zhǎng)度的需求，首先對(duì)信號(hào)使用功率譜密度分析，查看譜密度較大的頻率，將該頻率定為估計(jì)模態(tài)位置，按該頻率可計(jì)算出估計(jì)模態(tài)的周期，按其周期的1倍到10倍的數(shù)據(jù)長(zhǎng)度進(jìn)行模態(tài)辨識(shí)，對(duì)比模態(tài)辨識(shí)結(jié)果，分析模態(tài)辨識(shí)用的數(shù)據(jù)長(zhǎng)度要求。

研究發(fā)現(xiàn)，由于需要對(duì)時(shí)域縮放后的信號(hào)進(jìn)行重采樣，當(dāng)數(shù)據(jù)包含關(guān)心模態(tài)頻率5個(gè)周期以下時(shí)，重采樣后的數(shù)據(jù)構(gòu)成不了Hankel矩陣，因此無法進(jìn)行模態(tài)辨識(shí)，而當(dāng)數(shù)據(jù)長(zhǎng)度大于關(guān)心模態(tài)周期5個(gè)周期時(shí)，使用ERA方法可正常辨識(shí)模態(tài)，需要注意的是，數(shù)據(jù)長(zhǎng)度大于8倍關(guān)心模態(tài)頻率時(shí)，數(shù)據(jù)長(zhǎng)度過長(zhǎng)，在構(gòu)建Hankel矩陣時(shí)只使用到前段數(shù)據(jù)，而后段數(shù)據(jù)并未參與構(gòu)件Hankel矩陣，因此在使用時(shí)域縮放技術(shù)后采用ERA方法進(jìn)行模態(tài)辨識(shí)時(shí)，數(shù)據(jù)的長(zhǎng)度最有可能需要包含所需識(shí)別模態(tài)的6～8周期。

3.3 工作模態(tài)辨識(shí)

這里以一階模態(tài)辨識(shí)為例，研究時(shí)域縮放技術(shù)與ERA技術(shù)的結(jié)合對(duì)模態(tài)辨識(shí)的有利作用。

根據(jù)圖3典型測(cè)點(diǎn)功率譜密度分析結(jié)果，選擇頻帶10～40 Hz進(jìn)行時(shí)域縮放濾波，再根據(jù)圖1流程進(jìn)行模態(tài)辨識(shí)，辨識(shí)結(jié)果如表2所示，典型振型辨識(shí)結(jié)果如圖4所示，其中計(jì)算MAC值計(jì)算見式(9)。其中阻尼比理論值從1%～5%共設(shè)置了4個(gè)阻尼比，用于對(duì)比分析阻尼比理論值選擇大小對(duì)辨識(shí)結(jié)果的影響。

表2 一階模態(tài)辨識(shí)結(jié)果

(9)

圖4 模態(tài)振型辨識(shí)與理論計(jì)算結(jié)果對(duì)比

從以上辨識(shí)結(jié)果可以發(fā)現(xiàn)：

1) 采用時(shí)域縮放和ERA相結(jié)合的方法辨識(shí)可辨識(shí)出相應(yīng)的模態(tài)，并且具備一定的精度，該結(jié)合方法有效；

2) 從頻率辨識(shí)的結(jié)果可以看出，根據(jù)功率譜密度諧振峰參考方法以及過大阻尼篩選方法辨識(shí)可篩選出真實(shí)模態(tài)?；瑒?dòng)矩形窗后辨識(shí)出的模態(tài)頻率較為準(zhǔn)確，與理論模態(tài)頻率值的最大偏差為4.3%左右；

3) 從阻尼辨識(shí)的結(jié)果可以看出，阻尼比理論值越大，辨識(shí)結(jié)果與理論值偏差越小，尤其是各段數(shù)據(jù)辨識(shí)結(jié)果的平均值，從理論值1%～5% 4個(gè)工況，平均值偏差從175%減小至10%；

4) 從振型辨識(shí)的結(jié)果可以看出，滑動(dòng)矩形窗后篩選出的模態(tài)振型較為準(zhǔn)確，與理論模態(tài)振型值的MAC值均在0.99以上，最差為0.993左右，證明辨識(shí)模態(tài)與理論模態(tài)振型基本一致。

3.4 辨識(shí)算法對(duì)信號(hào)信噪比的魯棒性

根據(jù)以上分析，在3.3小節(jié)的基礎(chǔ)上，研究增加不同量級(jí)噪聲對(duì)模態(tài)頻率和阻尼比辨識(shí)結(jié)果的影響。這里對(duì)比了噪聲量級(jí)從信號(hào)標(biāo)準(zhǔn)差的0%～50%周期4種情況，辨識(shí)結(jié)果如表3所示。

表3 不同數(shù)據(jù)長(zhǎng)度的模態(tài)辨識(shí)結(jié)果

對(duì)各節(jié)點(diǎn)數(shù)據(jù)同時(shí)增加不同噪聲量級(jí)辨識(shí)模態(tài)時(shí)，從辨識(shí)結(jié)果可以看出：

1) 各階模態(tài)頻率辨識(shí)結(jié)果一致性較好，與理論計(jì)算結(jié)果較為一致；

2) 增加噪聲后，無論是低頻還是高頻模態(tài)，對(duì)模態(tài)阻尼比辨識(shí)結(jié)果取均值，各階模態(tài)阻尼比與理論值有一定的偏差，但偏差最大20%左右，隨著噪聲量級(jí)，阻尼比辨識(shí)結(jié)果比較穩(wěn)定，沒有出現(xiàn)辨識(shí)結(jié)果偏差明顯變大的情況，說明本文的方法對(duì)信號(hào)的信噪比的魯棒性較強(qiáng)。

4 結(jié)論

1) 采用時(shí)域縮放和ERA相結(jié)合的方法辨識(shí)可辨識(shí)出相應(yīng)的模態(tài)，并且具備一定的精度和魯棒性；

2) 在使用時(shí)域縮放技術(shù)后采用ERA方法進(jìn)行模態(tài)辨識(shí)時(shí)，數(shù)據(jù)的長(zhǎng)度最有可能需要包含所需識(shí)別模態(tài)的6～8周期；

3) 從頻率辨識(shí)的結(jié)果可以看出，根據(jù)功率譜密度諧振峰參考方法以及過大阻尼篩選方法辨識(shí)可篩選出真實(shí)模態(tài)，滑動(dòng)矩形窗后辨識(shí)出的模態(tài)頻率較為準(zhǔn)確；

4) 從阻尼辨識(shí)的結(jié)果可以看出，辨識(shí)阻尼結(jié)果偏差稍大，認(rèn)為增加測(cè)試噪聲后，辨識(shí)得出阻尼比最大偏差在20%左右；

5) 從振型辨識(shí)的結(jié)果可以看出，滑動(dòng)矩形窗后篩選出的模態(tài)振型較為準(zhǔn)確，與理論模態(tài)振型值的MAC值均在0.99以上，辨識(shí)模態(tài)與理論模態(tài)振型基本一致。

[1] SLAVIE.J，SIMONOVSKIL，BOLTEZARM.Damping Identification using a continuous wavelet transform:application to real data.Journal of Sound and Vibration[J].2003,21(2)：291-307.

[2] 譚冬梅，姚三，翟偉廉.振動(dòng)模態(tài)的參數(shù)識(shí)別綜述[J].華中科技大學(xué)學(xué)報(bào):自然科學(xué)版,2002，19(3):73-79.

[3] JOHNL CHRIDTIAN ASMUSSEN.Modal analysis Based on the Random Decrement Technique [D].Denmark:Aalborg University，1997.

[4] 鄭敏，申凡，陳同鋼.采用互相關(guān)復(fù)指數(shù)法進(jìn)行工作模態(tài)參數(shù)識(shí)別[J].南京理工大學(xué)學(xué)報(bào)，2002， 26(2):113-116.

[5] JAMES G H，CARNE T G.The Natural Excitation Technique(NExT)for modal parameter extraction from operating structures[J].The International Journal of Analytical and Experimental Modal Analysis，1995，10(4):260-277.

[6] 李惠彬.大型工程結(jié)構(gòu)模態(tài)參數(shù)識(shí)別技術(shù)[M].北京:北京理工大學(xué)出版社，2007:30-76.

[7] 孫曉蘭，王太勇.基于相關(guān)函數(shù)的振動(dòng)結(jié)構(gòu)工作模態(tài)參數(shù)識(shí)別方法[J].天津大學(xué)學(xué)報(bào)，2007，40(4):503-506.

[8] KAMMER DC.Sensor Placement for On-Orbit Modal Identification and Correlation of Large Space Structures[J].Journal of Guidance，Control，andDynamies,1991，14(2):251-258.

[9] 練繼建，李火坤，張建偉.基于奇異熵定階降噪的水工結(jié)構(gòu)振動(dòng)模態(tài)ERA識(shí)別方法[J].技術(shù)科學(xué)，2008，38(9):1398-1413.

[10]劉興漢，王躍宇.基于cholesky分解的改進(jìn)的隨機(jī)子空間法研究[J].宇航學(xué)報(bào)，2007，28(3):608-652.

[11]于開平，鄒經(jīng)湘.模態(tài)參數(shù)辨識(shí)的小波變換方法[J].宇航學(xué)報(bào)，1999，20(4):72-78.

[12]王彤,張令彌.工作模態(tài)分析的頻域空間域分解法及其應(yīng)用[J].航空學(xué)報(bào)，2006，27(1):62-66.

[13]黃琴,王彤,張海黎.基于隨機(jī)減量的運(yùn)行模態(tài)頻域分析方法[J].南京航空航天大學(xué)學(xué)報(bào)，2011，46(6):770-773.

[14]伍特輝.聲激勵(lì)下結(jié)構(gòu)模態(tài)參數(shù)識(shí)別的試驗(yàn)研究[D].南京:南京航空航天大學(xué),2009.

(責(zé)任編輯楊繼森)

Study on the Operational Mode Identification Based on ERA with the ZOOM Method

WANG Liang1，ZHANG Yan1，CAI Yi-peng1，NANGONG Zi-jun1, CHEN Guo-ping2

(1.China Academy of Launch Vehicle Technology, Beijing 100076, China; 2.Institute of Vibration Engineering Research, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Based on the ZOOM method, the operational mode identification was investigated under ambient excitation based on ERA method. Firstly, the theory of the ERA method was introduced. Secondly, the strategy for decomposition and implementation was put forward, including the ZOOM method and the modes’ filtration method. At last, an example as studied, where the cantilever beam was built and the white noise exciting was added. It is found that the mode can be realized by using the method, and the variety of the mode with the time can be present by using the sliding time window. A longer time period is most likely needed to capture at least 6～8 cycles of the mode for operational mode identification. The warp of the damping ratio identification becomes smaller as the theoretical value increases, where the warp changes from 175% to 10% as the theoretical value increases from 1% to 5%. The warps of the firstly and secondly order damping ratio identification are basically consistent. And the MAC of the model shape is all above 0.99.

mode identification; ZOOM; ERA; operational mode

2016-08-07；

2016-09-11

國家重點(diǎn)實(shí)驗(yàn)室2015年開放課題(MCMS-0115G01)；國防技術(shù)基礎(chǔ)科研項(xiàng)目(JSZL2015203B002)

王亮(1985—)，男，博士，高級(jí)工程師，主要從事導(dǎo)彈載荷與環(huán)境設(shè)計(jì)研究。

10.11809/scbgxb2017.01.023

王亮，張妍， 蔡毅鵬，等.基于時(shí)域縮放ERA方法的工作模態(tài)辨識(shí)技術(shù)[J].兵器裝備工程學(xué)報(bào)，2017(1):97-101.

format:WANG Liang，ZHANG Yan，CAI Yi-peng, et al.Study on the Operational Mode Identification Based on ERA with the ZOOM Method[J].Journal of Ordnance Equipment Engineering,2017(1):97-101.

TU311.3

A