趙國(guó)彥，鄧青林，李夕兵，董隴軍，陳光輝，張楚旋

?

基于EMD和形態(tài)分形維數(shù)的微震波形識(shí)別

趙國(guó)彥，鄧青林，李夕兵，董隴軍，陳光輝，張楚旋

(中南大學(xué)資源與安全工程學(xué)院，湖南長(zhǎng)沙，410083)

針對(duì)現(xiàn)有礦山微震監(jiān)測(cè)系統(tǒng)信號(hào)自動(dòng)識(shí)別難的問(wèn)題，提出基于經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和形態(tài)分形維數(shù)的識(shí)別方法。首先，采用EMD將原始信號(hào)分解為若干個(gè)本征模態(tài)分量(IMF)，選擇前5個(gè)分量進(jìn)行重構(gòu)得到新的信號(hào)。其次，求出處理后信號(hào)的形態(tài)學(xué)分形維數(shù)，利用微震波形和爆破波形分形維數(shù)的差異進(jìn)行信號(hào)識(shí)別。對(duì)50組微震波形和50組爆破波形進(jìn)行試驗(yàn)研究。對(duì)比未經(jīng)EMD處理的形態(tài)學(xué)分形維數(shù)以及經(jīng)EMD處理的盒維數(shù)識(shí)別結(jié)果。研究結(jié)果表明：50組微震波形和50組爆破波形在形態(tài)學(xué)分形維數(shù)為1.4時(shí)具有較高的識(shí)別率；微震波形維數(shù)主要在1.4以下，爆破波形維數(shù)則基本高于1.4；EMD結(jié)合形態(tài)學(xué)分形維數(shù)的識(shí)別效果最好，為微震監(jiān)測(cè)波形識(shí)別提供了新途徑。

微震；EMD；形態(tài)學(xué)；分形維數(shù)；波形識(shí)別

微震監(jiān)測(cè)技術(shù)通過(guò)對(duì)采集的地震波信息進(jìn)行分析處理確定破裂點(diǎn)的空間位置和震級(jí)，根據(jù)累積微震事件的數(shù)量和能量釋放的關(guān)系，估計(jì)因微震活動(dòng)引起的非線性應(yīng)變區(qū)域的發(fā)展，據(jù)此判斷巖體的穩(wěn)定性。該技術(shù)對(duì)于保障礦山人員安全及設(shè)備完整至關(guān)重要，目前已在國(guó)內(nèi)外采礦工程中推廣應(yīng)用[1?5]。然而，傳感器采集的原始數(shù)據(jù)往往摻雜著大量干擾的信息，需要后期人為剔除。為減輕人工勞動(dòng)量，近年來(lái)許多學(xué)者對(duì)波形的自動(dòng)識(shí)別進(jìn)行研究，極大地推動(dòng)了微震監(jiān)測(cè)技術(shù)的發(fā)展。歸結(jié)起來(lái)，微震波形的自動(dòng)識(shí)別特征主要包括頻譜特征、能量分布特征和參數(shù)化特征，分析方法主要有快速傅里葉變換、小波及小波包變換、經(jīng)驗(yàn)?zāi)B(tài)分解EMD等[6?11]。然而，一些分析方法本身存在缺陷，可能會(huì)帶來(lái)不準(zhǔn)確的分析結(jié)果，如傅里葉變換存在時(shí)域和頻域的局部化矛盾，小波及小波包變換則受小波基選擇以及信號(hào)分解層次的限制。根據(jù)信號(hào)本身的局部特征時(shí)間尺度，將1個(gè)復(fù)雜振動(dòng)信號(hào)分解為若干分量之和的經(jīng)驗(yàn)?zāi)B(tài)分解方法具有較大優(yōu)勢(shì)，該方法不需要對(duì)信號(hào)具備先驗(yàn)知識(shí)，是一種自適應(yīng)的信號(hào)處理方式。此外，地震波在巖體內(nèi)傳播時(shí)，由于巖層的散射、反射等效應(yīng)以及巖體對(duì)于地震波具有吸收作用，其能量與頻率是不斷衰減的[12]，因此，無(wú)論采用頻率識(shí)別方式還是能量識(shí)別方式，均存在較大的局限性和隨機(jī)性。但巖體微震與爆破振動(dòng)在觸發(fā)機(jī)制上存在本質(zhì)區(qū)別，直觀表達(dá)出來(lái)的形式就是它們的波形差異。分形維數(shù)則能定性、定量地刻畫這種非線性狀態(tài)，可用于識(shí)別這種差異，但要獲得準(zhǔn)確的分形維數(shù)則需要對(duì)信號(hào)進(jìn)行濾波處理。為此，本文作者在借助經(jīng)驗(yàn)?zāi)B(tài)分解構(gòu)造自適應(yīng)濾波器功能的前提下，提取微震及爆破波形的主要信息，并且選擇精確度較高的形態(tài)學(xué)分形算法對(duì)微震和爆破波形的分形維數(shù)進(jìn)行估計(jì)。

1 EMD方法

對(duì)于任意的復(fù)雜信號(hào)()，HUANG等[13]指出通過(guò)其本身特征時(shí)間尺度來(lái)獲得一組本征模態(tài)函數(shù)IMF，稱為經(jīng)驗(yàn)?zāi)B(tài)分解篩選的方法，其步驟如下。

1) 確定原始信號(hào)()的所有的局部極大值點(diǎn)、極小值點(diǎn)，并分別用3次樣條曲線擬合形成上、下包 絡(luò)線。

2) 求出原始的數(shù)據(jù)序列與上下包絡(luò)線均值()之差：

3) 檢驗(yàn)式(1)得到的()是否滿足IMF的條件。若不滿足，則將其視為新的()，重復(fù)以上步驟直到滿足IMF的條件為止，此時(shí)()為第1階IMF，記為1()，通常包含信號(hào)的最高頻成分。

4) 求出原始信號(hào)()與1()之差得到剩余序列1()：

5) 將1()作為新的原始信號(hào)，重復(fù)以上步驟，直到第階的剩余信號(hào)R()成為單調(diào)函數(shù)或者比預(yù)期值小，不能再篩分為止，即可得到一系列IMF分量1()，2()，…，n()以及剩余項(xiàng)n()，這樣初始信號(hào)()就分解成若干IMF分量與1個(gè)殘余項(xiàng)之和，即

(3)

式中：余項(xiàng)n()為信號(hào)的趨勢(shì)項(xiàng)，各IMF分量則表示信號(hào)從高到低不同頻率段的成分。

EMD將任意復(fù)雜的信號(hào)自適應(yīng)地分解為一組本征模態(tài)分量IMF，需要滿足以下2個(gè)條件：

1) 在整個(gè)時(shí)間數(shù)據(jù)序列之中，極值點(diǎn)的數(shù)量(包括極大值與極小值)與過(guò)零點(diǎn)的數(shù)量至多相差1個(gè)。

2) 在任意1個(gè)時(shí)間點(diǎn)上，由局部極值點(diǎn)構(gòu)成的包絡(luò)線平均值為0(局部極大值點(diǎn)確定上包絡(luò)線，局部極小值點(diǎn)確定下包絡(luò)線)。

在EMD分解過(guò)程中，因?qū)嶋H的篩分很難保證局部均值為0，而過(guò)多的重復(fù)篩選會(huì)導(dǎo)致IMF失去實(shí)際的物理意義，可通過(guò)限制標(biāo)準(zhǔn)偏差來(lái)實(shí)現(xiàn)篩分停止：

式中：為信號(hào)的時(shí)間長(zhǎng)度；I?1()和I()為2個(gè)連續(xù)處理結(jié)果的時(shí)間序列；D通常取0.2~0.3。

2 形態(tài)學(xué)分形維數(shù)及其計(jì)算

2.1 分形維數(shù)

分形是自然界復(fù)雜事物存在的普遍特征。非平穩(wěn)隨機(jī)振動(dòng)信號(hào)的不規(guī)則程度能夠以分形維數(shù)來(lái)描述，計(jì)算分形維數(shù)的常見(jiàn)方法有盒維數(shù)、網(wǎng)格維數(shù)、關(guān)聯(lián)維數(shù)等，其中基于“覆蓋”思想的算法簡(jiǎn)單清晰。該方法假設(shè)1個(gè)實(shí)平面集合2，和分別為集合中1個(gè)非空子集和1個(gè)尺度為的閉集，則的分形維數(shù)定義為[14]

式中：()為覆蓋所需的個(gè)數(shù)。

2.2 形態(tài)學(xué)分形算法

估計(jì)分形維數(shù)的關(guān)鍵是在不同尺度下對(duì)信號(hào)進(jìn)行度量，利用多尺度形態(tài)學(xué)可以實(shí)現(xiàn)[15]。在集合覆蓋信號(hào)過(guò)程中，以其上確界函數(shù)即結(jié)構(gòu)元素()進(jìn)行等效變換，分析尺度范圍=1, 2, …,，其算法包含膨脹及腐蝕2種基本算子。對(duì)于一維離散時(shí)間信號(hào)()(=0, 1, 2, …)，每個(gè)尺度的膨脹和腐蝕結(jié)果分別為：

(7)

(9)

式中：M為信號(hào)的Minkowski-Bouligand維數(shù)；為常數(shù)。對(duì)上式和進(jìn)行最小二乘擬合，得到的直線斜率為M，即為最終所要求得的()分形維數(shù)。

2.3 參數(shù)選擇

一般地，選擇單位結(jié)構(gòu)元素()為{0,0,0}，因該結(jié)構(gòu)在保證維數(shù)估計(jì)不受信號(hào)幅值范圍影響的同時(shí)，還能使算法的運(yùn)算量減少。最大分析尺度原則上為小于/2的正整數(shù)(為離散信號(hào)采樣點(diǎn)數(shù))，當(dāng)數(shù)據(jù)長(zhǎng)度比較大時(shí)，適當(dāng)減小可減少計(jì)算量[16]，本文取256。

3 基于EMD和形態(tài)分形維數(shù)的微震波形識(shí)別

3.1 原理及步驟

從本質(zhì)上來(lái)講，礦山微震以及爆破產(chǎn)生的波形均屬于非平穩(wěn)隨機(jī)振動(dòng)(如圖1和圖2所示)，因而其本身已經(jīng)具備可估計(jì)分形維數(shù)的特點(diǎn)。然而，由于礦井下復(fù)雜作業(yè)環(huán)境的影響，系統(tǒng)采集到的波形往往受到各種噪聲的污染，因此，直接對(duì)原始波形進(jìn)行分形維數(shù)估計(jì)會(huì)得到不準(zhǔn)確的結(jié)果，必須事先對(duì)波形進(jìn)行去噪預(yù)處理。

圖1 典型的巖體微震波形

圖2 典型的爆破振動(dòng)波形

經(jīng)驗(yàn)?zāi)B(tài)分解EMD是一種無(wú)需先驗(yàn)知識(shí)的非平穩(wěn)振動(dòng)信號(hào)處理方法，可將信號(hào)自適應(yīng)地分解為有限個(gè)IMF，IMF按信息重要程度由高到低階排列，且由于EMD算法本身還存在一定缺陷，容易產(chǎn)生虛假分量[17?18]，因此，本文在結(jié)合大量試驗(yàn)結(jié)果的前提下選擇前5階IMF提取信號(hào)的主要信息，作為去噪后的最終信號(hào)。綜上所述，基于EMD和形態(tài)分形維數(shù)的微震波形識(shí)別步驟為：

1) 利用EMD將原始信號(hào)()分解為若干個(gè)IMF；

由此計(jì)算出的典型微震與爆破波形形態(tài)學(xué)分形維數(shù)分別為1.292 9和1.543 2，可見(jiàn)兩者差異較大。

3.2 實(shí)例應(yīng)用

為尋找微震波形與爆破波形的分形特征差異，確定兩者形態(tài)學(xué)維數(shù)的定量差異，分別從原始波形數(shù)據(jù)庫(kù)隨機(jī)選取50組微震事件波形、50組爆破事件波形，利用matlab編程，計(jì)算各個(gè)波形的形態(tài)學(xué)分形維數(shù)，見(jiàn)表1。

表1 微震波形與爆破波形形態(tài)學(xué)分形維數(shù)

從表1可以看出：微震事件波形的形態(tài)學(xué)分形維數(shù)基本上低于爆破事件波形的形態(tài)學(xué)分形維數(shù)。因此，可以用形態(tài)學(xué)分形維數(shù)來(lái)辨識(shí)這2類波形。圖3所示為微震及爆破波形的形態(tài)學(xué)分形維數(shù)分布圖。

圖3 微震與爆破波形的形態(tài)學(xué)分形維數(shù)分布

從圖3可知：微震波形與爆破波形的形態(tài)學(xué)分形維數(shù)為1.4時(shí)呈現(xiàn)明顯分界，分形維數(shù)在1.4以上以爆破事件為主，1.4以下以微震事件為主。若以形態(tài)學(xué)分形維數(shù)1.4作為辨識(shí)微震與爆破事件的判別因子，那么，在抽取的樣本中，分形維數(shù)高于1.4的微震事件有4個(gè)，識(shí)別正確率為92%，而分形維數(shù)低于1.4的爆破事件則有3個(gè)，識(shí)別正確率為94%；總體事件識(shí)別的正確率達(dá)到93%。從波形復(fù)雜度的角度看，微震事件產(chǎn)生的波形大多數(shù)較單一，而井下生產(chǎn)爆破多有多段微差延時(shí)，波形會(huì)相對(duì)復(fù)雜。因此，以形態(tài)學(xué)分形維數(shù)1.4作為該礦山微震系統(tǒng)波形識(shí)別因子是可行的。

3.3 結(jié)果驗(yàn)證

為驗(yàn)證形態(tài)學(xué)分形維數(shù)辨識(shí)的效果，一方面，以未經(jīng)EMD處理的原始波形形態(tài)學(xué)分形維數(shù)比較，另一方面，選擇分形常用的盒維數(shù)計(jì)算方法進(jìn)行對(duì)比。圖4和圖5所示分別為未處理的原始波形形態(tài)學(xué)分形維數(shù)分布圖以及用盒維數(shù)代替形態(tài)學(xué)分形維數(shù)的微震與爆破波形維數(shù)分布圖。

從圖4和圖5可見(jiàn)：微震與爆破原始波形的形態(tài)學(xué)維數(shù)以及經(jīng)EMD的分形盒維數(shù)在總體上也呈爆破高、微震低的趨勢(shì)，但無(wú)論是以形態(tài)學(xué)分形維數(shù)1.4~1.5還是以盒維數(shù)1.30~1.32中任意1個(gè)維數(shù)作為兩者的識(shí)別因子，其總體正確識(shí)別率均低于90%，說(shuō)明不經(jīng)過(guò)預(yù)處理的波形以及采用本身估計(jì)精度較低的盒維數(shù)計(jì)算方法，均可能影響最終結(jié)果，使得高維數(shù)微震波形和低維數(shù)爆破波形容易產(chǎn)生混淆。以上分析結(jié)果表明了基于EMD的形態(tài)學(xué)分形維數(shù)識(shí)別微震和爆破波形的有效性及合理性。

圖4 未經(jīng)EMD處理的微震與爆破波形的形態(tài)學(xué)分形維數(shù)分布

圖5 微震與爆破波形的分形盒維數(shù)分布

4 結(jié)論

1) EMD具有自適應(yīng)濾波器的功能，能夠提取振動(dòng)信號(hào)的主要特征與信息，從而提高了原始信號(hào)的信噪比。

2) 分形維數(shù)具有刻畫不同波形特征的能力?；谛螒B(tài)學(xué)分形維數(shù)的算法估計(jì)精度高，對(duì)經(jīng)EMD處理的100組樣本進(jìn)行試驗(yàn)，求出微震波形與爆破波形的識(shí)別維數(shù)為1.4，且識(shí)別正確率高達(dá)93%，相比未經(jīng)EMD處理的原始波形形態(tài)學(xué)分形維數(shù)以及經(jīng)EMD處理的盒維數(shù)識(shí)別，更加合理有效，可為后續(xù)大量樣本的模式自動(dòng)識(shí)別提供一種新的方法。

[1] SALVONI M, DIGHT P M. Rock damage assessment in a large unstable slope from microseismic monitoring-MMG century mine (Queensland, Australia) case study[J]. Engineering Geology, 2016, 210: 45?56.

[2] URBANCIC T I, TRIFU C I. Recent advances in seismic monitoring technology at Canadian mines[J]. Journal of Applied Geophysics, 2000, 45(4): 225?237.

[3] POTVIN Y. Strategies and tactics to control seismic risks in mines[J]. The Journal of the Southern African Institute of Mining and Metallurgy, 2009, 109(3): 177?186.

[4] 楊承祥, 羅周全, 唐禮忠. 基于微震監(jiān)測(cè)技術(shù)的深井開(kāi)采地壓活動(dòng)規(guī)律研究[J].巖石力學(xué)與工程學(xué)報(bào), 2007, 26(4): 818?824.YANG Chengxiang, LUO Zhouquan, TANG Lizhong. Study on rule of geostatic activity based on microseismic monitoring technique in deep mining[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4): 818?824.

[5] 董隴軍, 李夕兵, 唐禮忠, 等. 無(wú)需預(yù)先測(cè)速的微震震源定位的數(shù)學(xué)形式及震源參數(shù)確定[J].巖石力學(xué)與工程學(xué)報(bào), 2011, 30(10): 2057?2067. DONG Longjun, LI Xibing, TANG Lizhong, et al. Mathematical functions and parameters for microseismic source location without pre-measuring speed[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 2057?2067.

[6] ZHAO Guoyan, MA Ju, DONG Longjun, et al. Classification of mine blasts and microseismic events using starting-up features in seismograms[J]. Transactions of Nonferrous Metals Society of China, 2015, 25(10): 3410?3420.

[7] 趙國(guó)彥, 鄧青林, 馬舉. 基于FSWT時(shí)頻分析的礦山微震信號(hào)分析與識(shí)別[J]. 巖土工程學(xué)報(bào), 2015, 37(2): 306?312. ZHAO Guoyan, DENG Qinglin, MA Ju. Recognition of mine microseismic signals based on FSWT time-frequency analysis[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(2): 306?312.

[8] 陸菜平, 竇林名, 吳興榮, 等. 巖體微震監(jiān)測(cè)的頻譜分析與信號(hào)識(shí)別[J]. 巖土工程學(xué)報(bào), 2005, 27(7): 772?775. LU Caiping, DOU Linming, WU Xingrong, et al. Frequency spectrum analysis on microseismic monitoring and signal differentiation of rock material[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(7): 772?775.

[9] 朱權(quán)潔, 姜福興, 于正興, 等. 爆破震動(dòng)與巖石破裂微震信號(hào)能量分布特征研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2012, 31(4): 723?730. ZHU Quanjie, JIANG Fuxing, YU Zhengxing, et al. study on energy distribution characters about blasting vibration and rock fracture microseismic signal[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(4): 723?730.

[10] 朱權(quán)潔, 姜福興, 尹永明, 等. 基于小波分形特征與模式識(shí)別的礦山微震波形識(shí)別研究[J]. 巖土工程學(xué)報(bào), 2012, 34(11): 2036?2042. ZHU Quanjie, JIANG Fuxing, YIN Yongming, et al. Classification of mine microseismic events based on wavelet-fractal method and pattern recognition[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(11): 2036?2042.

[11] 吳響, 錢建生, 王海燕, 等. 微震信號(hào)多尺度非線性特征提取與辨識(shí)研究[J]. 儀器儀表學(xué)報(bào), 2014, 35(5): 969?975. WU Xiang, QIAN Jiansheng, WANG Haiyan, et al. Study on multi-scale nonlinear feature extraction and signal identification for microseismic signal[J]. Chinese Journal of Scientific Instrument, 2014, 35(5): 969?975.

[12] 李振春, 王清振. 地震波衰減機(jī)理及能量補(bǔ)償研究綜述[J]. 地球物理學(xué)進(jìn)展, 2007, 22(4): 1147?1152. LI Zhenchun, WANG Qingzhen. A review of research on mechanism of seismic attenuation and energy compensation[J]. Progress in Geophysics, 2007, 22(4): 1147?1152.

[13] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[C]// Proceedings of the Royal Society of A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1998: 903?995.

[14] 李兵, 張培林, 任國(guó)全, 等. 基于數(shù)學(xué)形態(tài)學(xué)的分形維數(shù)計(jì)算及在軸承故障診斷中的應(yīng)用[J].振動(dòng)與沖擊, 2010, 29(5): 191?194. LI Bing, ZHANG Peilin, REN Guoquan, et al. Mathematic morphology based fractal dimension calculation and its application in fault diagnosis of roller bearings[J]. Journal of Vibration and Shock, 2010, 29(5): 191?194.

[15] 張亢, 程軍圣, 楊宇. 基于局部均值分解與形態(tài)學(xué)分形維數(shù)的滾動(dòng)軸承故障診斷方法[J].振動(dòng)與沖擊, 2013, 32(9): 90?94. ZHANG Kang, CHENG Junsheng, YANG Yu. Roller bearing fault diagnosis based on local mean decomposition and morphological fractal dimension[J]. Journal of Vibration and Shock, 2013, 32(9): 90?94.

[16] 王鳳利, 段樹(shù)林, 于洪亮, 等. 基于EEMD和形態(tài)學(xué)分形維數(shù)的柴油機(jī)故障診斷[J].內(nèi)燃機(jī)學(xué)報(bào), 2012, 30(6): 557?562. WANG Fengli, DUAN Shulin, YU Hongliang, et al. Fault diagnosis of diesel engine based on EEMD and morphology fractal dimension[J]. Transactions of CSICE, 2012, 30(6): 557?562.

[17] 韓東穎, 李庚, 時(shí)培明. 基于EMD和分形盒維數(shù)的旋轉(zhuǎn)機(jī)械耦合故障診斷方法研究[J].振動(dòng)與沖擊, 2013, 32(15): 209?214. HAN Dongying, LI Geng, SHI Peiming. Diagnosis method for coupling faults of rotary machinery based on EMD and fractal box dimension[J]. Journal of Vibration and Shock, 2013, 32(15): 209?214.

[18] 程軍圣, 于德介, 楊宇. 基于EMD和分形維數(shù)的轉(zhuǎn)子系統(tǒng)故障診斷[J].中國(guó)機(jī)械工程, 2005, 16(12): 1088?1091. CHENG Junsheng, YU Dejie, YANG Yu. Fault diagnosis for rotor system based on EMD and fractal dimension[J]. China Mechanical Engineering, 2005, 16(12): 1088?1091.

(編輯 陳燦華)

Recognition of microseismic waveforms based on EMD and morphological fractal dimension

ZHAO Guoyan, DENG Qinglin, LI Xibing, DONG Longjun, CHEN Guanghui, ZHANG Chuxuan

(School of Resources and Safety Engineering, Central South University, Changsha 410083, China)

Considering that it is difficult to recognize signals automatically in microseismic monitoring system of mine, the new recognition method was proposed based on empirical mode decomposition (EMD) and morphological fractal dimension. Firstly, origined-signal was decomposed to some intrinsic mode functions (IMF) by EMD, and new signal was obtained by adding former five components. Then, morphological fractal dimension of the new signal was calculated to recognize microseismic waveform and blasting waveform. The 50 microseimic waveforms and 50 blasting waveforms were tested. Finally, the dimensions were calculated with two different methods, i.e. morphological fractal dimension without EMD and fractal box dimension with EMD. The results show that 50 microseimic waveforms and 50 blasting waveforms show high recognition ratio in morphological fractal dimension of 1.4. Dimensions of microseismic waveform and blasting waveform are mainly below and over 1.4, respectively. The new recognition method has the best recognition result, and provides a new way for waveform recognition.

microseimic; EMD; morphology; fractal dimension; waveform recognition

10.11817/j.issn.1672-7207.2017.01.023

TD76；TD235

A

1672?7207(2017)01?0162?06

2016?01?10；

2016?03?24

國(guó)家自然科學(xué)基金資助項(xiàng)目(51374244) (Project(51374244) supported by the National Natural Science Foundation of China)

趙國(guó)彥，教授，博士生導(dǎo)師，從事采礦工程、礦山安全和巖石力學(xué)與工程等研究；E-mail: 810703752@qq.com