秦 晅,蔡建超,劉少勇,卞愛飛
(中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,地球內(nèi)部多尺度成像湖北省重點實驗室,湖北武漢430074)
基于經(jīng)驗?zāi)B(tài)分解互信息熵與同步壓縮變換的微地震信號去噪方法研究
秦 晅,蔡建超,劉少勇,卞愛飛
(中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,地球內(nèi)部多尺度成像湖北省重點實驗室,湖北武漢430074)
針對微地震信號具有隨機性、非平穩(wěn)性與時頻耦合的特點以及經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)的模態(tài)混疊問題,提出了基于經(jīng)驗?zāi)B(tài)分解互信息熵與同步壓縮變換(Synchrosqueezing Transform,SST)的微地震信號去噪方法。首先對微地震信號進行經(jīng)驗?zāi)B(tài)分解,獲得從高頻到低頻排列的固有模態(tài)函數(shù)(Intrinsic Mode Function,IMF)分量;然后求取相鄰固有模態(tài)函數(shù)分量之間的互信息熵,從而辨識出高頻與低頻部分的分界;最后利用同步壓縮變換提取高頻部分的有效信號,將其與低頻部分重構(gòu),實現(xiàn)微地震信號的有效去噪。利用不同噪聲強度的理論模型和實際資料,對本文方法與直接舍棄高頻成分的去噪方法進行了對比,結(jié)果表明,本文方法能夠很好地去除微地震信號中的混疊噪聲,并將有效信號從噪聲中提取出來,提高了資料的信噪比。
微地震信號;經(jīng)驗?zāi)B(tài)分解;同步壓縮變換;互信息熵;重構(gòu);噪聲壓制
微地震監(jiān)測技術(shù)作為獲得地層裂縫分布的一種較為有效方法,已經(jīng)在油氣田開發(fā)與煤礦開采等領(lǐng)域得到了廣泛應(yīng)用[1]。由于微地震波能量微弱,檢波器采集到的微地震信號存在大量隨機噪聲,因而后續(xù)的初至拾取、震源定位及裂縫解釋工作會受到很大影響[2-3],因此,需要對微地震信號進行濾波、去噪等處理,將微地震有效信號從噪聲中分離出來。
傳統(tǒng)的傅里葉變換是信號去噪的有效方法,但對于具有非平穩(wěn)信號特點的微地震信號去噪效果不佳[4-5]。連續(xù)小波變換在低頻部分具有較高的頻率分辨率,在高頻部分具有較高的時間分辨率,能夠在時頻域很好地區(qū)分非平穩(wěn)信號的突變部分,但去噪結(jié)果受到Heisenberg不確定性和小波基函數(shù)選取的影響[6-9]。經(jīng)驗?zāi)B(tài)分解(EMD)是一種信號自適應(yīng)分解方法,已經(jīng)在生物醫(yī)學(xué)的信號處理、機械設(shè)備的檢測與故障診斷、地震勘探的信號處理[10-14]等領(lǐng)域得到了廣泛的應(yīng)用,但其穩(wěn)定性會受到模態(tài)混疊現(xiàn)象的影響[15-16]。
針對EMD中出現(xiàn)的模態(tài)混疊現(xiàn)象,WU等[15]在信號中加入功率譜密度均勻分布的白噪聲,使其在不同尺度上具有連續(xù)性的特征,給出了總體平均經(jīng)驗?zāi)B(tài)分解方法;YEH等[16]通過加入正、負成對的輔助噪聲,提出了互補集合經(jīng)驗?zāi)B(tài)分解方法。上述方法能在一定程度上抑制模態(tài)混疊現(xiàn)象,但計算的復(fù)雜程度會顯著上升。
本文根據(jù)微地震信號具有隨機性、非平穩(wěn)性與時頻耦合的特點,兼顧EMD出現(xiàn)的模態(tài)混疊問題,提出了基于EMD互信息熵與同步壓縮變換(SST)的微地震信號去噪方法。模型測試和實際資料處理驗證了方法的有效性與實用性。
在實際微地震信號去噪中,有效信號往往位于低頻部分,所以通常把高頻成分的固有模態(tài)函數(shù)(IMF)分量作為噪聲直接去除,重構(gòu)原信號剩余的IMF分量即可實現(xiàn)壓制噪聲的目的[11,17-18]。但是,簡單地舍棄高頻成分會丟失部分有效信息,造成去噪效果不佳,因此需要發(fā)展一種有效剔除混疊噪聲并保留有效信號的方法。
1.1 IMF互信息熵成分辨識
EMD算法將時間序列信號自適應(yīng)地分解成一組含有不同尺度的IMF分量,能夠較好地處理非平穩(wěn)信號,它不依賴于基函數(shù)的選取。EMD得到的IMF必須滿足以下兩個條件[11-16]:①IMF中極值點的個數(shù)與過零點的個數(shù)相等或最多相差一個;②在任意數(shù)值點處,由局部極大值與極小值確定的包絡(luò)線均值為0。
經(jīng)過EMD后,微地震信號s(t)可以表示為:
(1)
式中:SIMFi(t)表示頻率從高到低排列的IMF分量;rn(t)表示殘差分量,代表信號的平均趨勢。將公式(1) 改寫為:
(2)
互信息熵用于衡量兩個隨機變量之間的統(tǒng)計依存性[18-21],其表達式為:
(3)
式中:p(xi,yj)為聯(lián)合概率分布;p(xi)和p(yj)為邊緣概率分布;X和Y表示不同的IMF分量;r和s分別是X和Y的符號數(shù)量。
假設(shè)高頻成分為噪聲干擾,低頻成分為有效信號,利用IMF分量之間的互信息熵關(guān)系可以辨識出高頻與低頻成分的分界。由EMD得到的IMF分量特征可知:高頻噪聲對每一個IMF分量的依存性逐漸降低,而低頻有效信號對每一個IMF分量的依存性逐漸增強。因此,可以假設(shè)高頻成分與低頻成分是部分相互統(tǒng)計獨立的。由互信息熵的特點可知:兩個相互獨立的隨機變量之間的互信息熵應(yīng)等于0[18,22]。所以,在計算每個相鄰IMF分量之間的互信息熵時會出現(xiàn)局部極小值,只需搜索第一個局部極小值即得到高頻成分與低頻成分的分界。于是可得以下搜索目標函數(shù):
(4)
式中:k為高頻成分與低頻成分分界處的IMF分量序號。
1.2 基于同步壓縮變換的混疊噪聲壓制改進算法
SST算法是在小波變換的基礎(chǔ)上發(fā)展起來的一種時頻域重排算法,但與原始重排算法不同的是,其在提高時頻分辨率的同時,也能夠支持信號的重構(gòu)[23-27]。由于SST能夠獲得較高頻率精度的時頻譜,使得各頻率曲線之間不存在交叉項,因此可以較好地解決頻率的混疊問題。利用SST實現(xiàn)高頻成分混疊噪聲壓制的步驟如下:
1) 對高頻成分H(t)進行連續(xù)小波變換,得到小波變換系數(shù)W(ωl,b),并計算其瞬時頻率ω(a,b)。
(5)
式中:n表示信號采樣點數(shù)。
3) 利用計算得到的瞬時頻率建立(b,a)→(b,ω(a,b))之間的映射關(guān)系,對時間-尺度平面的能量進行重新分配,將其轉(zhuǎn)化為時間-頻率平面,使得中心頻率ωl相鄰區(qū)間[ωl-Δω/2,ωl+Δω/2]的值被壓縮到ωl上[28],獲得SST系數(shù)T(ωl,b),其離散形式為
(6)
4) 利用SST獲得的高分辨率時頻譜確定有效信號的時間與頻率范圍,再對其進行SST重構(gòu),獲得信號x(t),實現(xiàn)對高頻成分混疊噪聲的壓制。
(7)
將高頻成分混疊噪聲壓制后的信號x(t)與剩余的低頻成分有效信號進行重構(gòu),即可實現(xiàn)對微地震信號噪聲的去除。
2.1 模型測試
為了驗證本文去噪方法的可行性,采用如圖1所示的井中微地震觀測系統(tǒng),利用有限差分正演方法模擬了一個井中16級檢波器接收的微地震記錄(圖2),數(shù)據(jù)采樣間隔為0.5ms,采樣點數(shù)為1000。將圖2模擬的微地震記錄加入高斯白噪聲,加噪后數(shù)據(jù)的信噪比為0.3,如圖3 所示,從中抽取一道數(shù)據(jù)進行去噪處理。圖4為抽取的單道記錄加噪前、后波形及其頻譜,對單道加噪信號進行自適應(yīng)EMD,得到8個IMF分量以及對應(yīng)的頻譜,如圖5所示??梢钥闯雒恳粋€IMF分量都有不同的振幅,而且分解的頻率按照從高到低的順序排列。
利用互信息熵公式(3)計算相鄰IMF分量之間的互信息量,結(jié)果如表1所示。按照本文1.1節(jié)所給的高頻部分與低頻部分分界點判定方法,可知IMF34為出現(xiàn)的第一個極小值,所以IMF4即為加噪信號的高頻部分與低頻部分分界點,圖5中IMF4的波形以及對應(yīng)的頻譜也驗證了判定結(jié)果的正確性。圖6為將IMF4之后的分量疊加獲得的低頻部分信號和SST的結(jié)果,從中可以看出低頻部分的噪聲含量較低,可以將低頻部分當作有效信號進行處理。
圖1 井中微地震觀測系統(tǒng)
圖2 微地震記錄
圖3 加噪后微地震記錄
圖4 原信號、含噪信號及其頻譜a 原信號; b 原信號頻譜; c 含噪信號; d 含噪信號頻譜
圖5 含噪信號EMD及其頻譜
IMF分量IMF12IMF23IMF34IMF45IMF56IMF67IMF78系數(shù)0.93290.98640.57790.68481.43800.86171.0580
利用1.2節(jié)給出的SST混疊噪聲壓制方法,對IMF4之前的分量疊加組成的高頻部分進行同步壓縮變換,去除混疊噪聲,如圖7所示。圖7a為IMF4之前的分量疊加組成的高頻部分信號,從中可以看出噪聲基本上分布在高頻部分。圖7b和圖7c分別為高頻部分短時傅里葉變換(SIFT)和SST時頻譜,可以看出,SIFT時頻譜模糊,不能很好地分辨出有效信號,而SST時頻譜具有較高的頻率精度,各時頻曲線間不存在交叉項,提高了有效信號脊線的分辨率,也證明了SST能夠較好地改善頻率混疊現(xiàn)象。圖7d為SST自適應(yīng)閾值去噪后的時頻譜,可以看出高頻部分的噪聲得到了很好的壓制。圖7e為對高頻部分有效信號進行提取后得到的有效信號,將其與剩余的低頻部分進行重構(gòu),即可得到降噪后的信號。
圖8對比了單道加噪信號直接舍棄高頻分量的去噪結(jié)果和本文方法的去噪結(jié)果,圖9對比了圖3所示的微地震加噪記錄直接舍棄高頻分量的去噪結(jié)果和本文方法的去噪結(jié)果??梢钥闯?直接舍棄高頻分量去噪的方法雖然能夠壓制噪聲,但同時使有效信號頻率出現(xiàn)部分缺失,并且還混疊了部分噪聲,造成了有效信號的波形畸變與衰減;而本文方法不僅很好地壓制了噪聲,而且使有效信號的波形得到了較好保持。定義去噪后信號與原信號的差值為含有的噪聲,信號與噪聲能量的比值為信噪比RSN[27,29]:
圖6 低頻部分信號(a)和同步壓縮變換結(jié)果(b)
圖7 同步壓縮變換混疊噪聲壓制效果a 高頻部分信號; b SIFT時頻譜; c SST時頻譜; d SST自適應(yīng)閾值去噪后時頻譜; e SST提取的有效信號
圖8 單道加噪信號應(yīng)用不同方法去噪的效果對比a 直接舍棄高頻分量; b 本文方法
圖9 加噪記錄應(yīng)用不同方法去噪的效果對比a 直接舍棄高頻分量; b 本文方法
(8)
式中:M表示原信號;S為去噪后信號。由公式(8)計算直接舍棄高頻分量方法去噪后的信噪比為2.43,本文方法去噪后的信噪比為6.67,定量說明了本文方法去噪的有效性。
不斷地改變含噪信號的信噪比,觀察不同噪聲水平下兩種方法的去噪效果,并計算去噪后的信噪比,結(jié)果如表2所示??梢钥闯?在不同噪聲水平下,兩種方法都能提高含噪信號的信噪比,但是本文方法具有更好的去噪效果。
2.2 應(yīng)用實例
本文截取某油田壓裂的一段井中微地震有效事件數(shù)據(jù)進行了去噪效果分析,如圖10所示。原始數(shù)據(jù)由16級三分量檢波器接收,采樣間隔為0.25ms,采樣點數(shù)為1000。圖10a中1~16,17~32,33~48道分別為數(shù)據(jù)的z,x,y分量,由于數(shù)據(jù)的信噪比較低,噪聲幾乎淹沒了有用的微地震信號,無法清晰地識別出P波、S波。圖10b為對原始數(shù)據(jù)進行EMD后直接去除高頻成分的結(jié)果,可以看出,雖然高頻噪聲被濾除,但是微地震有效信號同相軸仍然不夠清晰,有效信號波形發(fā)生了畸變。圖10c為利用SST方法[27]提取的結(jié)果,可以看出,SST方法雖然較好地提取了有效信號,但是混疊了部分噪聲,說明噪聲壓制還不夠徹底。圖10d為利用本文方法去噪后的結(jié)果,原始數(shù)據(jù)的噪聲得到了很好壓制,有效信號的波形也得到了很好保持,為后續(xù)的初至拾取、震源定位及裂縫解釋工作奠定了基礎(chǔ)。
表2 不同噪聲水平下不同方法去噪后的信噪比
圖10 實際微地震數(shù)據(jù)去噪效果對比a 原始數(shù)據(jù); b 直接舍棄高頻成分去噪結(jié)果; c SST提取結(jié)果; d 本文方法去噪結(jié)果
本文針對微地震信號具有隨機性、非平穩(wěn)性與時頻耦合的特點,以及EMD存在的模態(tài)混疊問題,利用同步壓縮變換能提高信號時頻分辨率,同時支持信號重構(gòu)的優(yōu)點,提出了基于EMD互信息熵與同步壓縮變換的微地震信號去噪方法。不同噪聲強度的理論模型測試和實際資料分析表明,該方法能夠很好地去除信號中的混疊噪聲,較好地將有效信號從噪聲中提取出來,提高資料信噪比,為后續(xù)的初至拾取、震源定位及裂縫解釋等工作奠定了基礎(chǔ)。
[1] MAXWELL S.Microseismic:growth born from success[J].The Leading Edge,2010,29(3):338-343
[2] 宋維琪,何可,郭全仕,等.微地震資料自適應(yīng)濾波方法研究[J].石油物探,2013,52(3):229-233 SONG W Q,HE K,GUO Q S,et al.Micro-seismic data adaptive filtering method[J].Geophysical Prospecting for Petroleum,2013,52(3):229-233
[3] 楊心超,朱海波,崔樹果,等.P波初動震源機制解在水力壓裂微地震監(jiān)測中的應(yīng)用[J].石油物探,2015,54(1):43-50 YANG X C,ZHU H B,CUI S G,et al.Application of P-wave first-motion focal mechanism solutions in microseismic monitoring for hydraulic fracturing[J].Geophysical Prospecting for Petroleum,2015,54(1):43-50
[4] ALVANITOPOULOS P F,PAPAVASILEIOU M,ANDREADIS I,et al.Seismic intensity feature construction based on the Hilbert-Huang transform[J].IEEE Transactions on Instrumentation and Measurement,2012,61(2):326-337
[5] GACI S.The use of wavelet-based denoising techniques to enhance the first-arrival picking on seismic traces[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(8):4558-4563
[6] SINHA S,ROUTH P S,ANNO P D,et al.Spectral decomposition of seismic data with continuous-wavelet transform[J].Geophysics,2005,70(6):19-25
[7] BEENAMOL M,PRABAVATHY S,MOHANALIN J.Wavelet based seismic signal de-noising using Shannon and Tsallis entropy[J].Computers & Mathematics with Applications,2012,64(11):3580-3593
[8] 盛冠群,李振春,王維波,等.基于小波分解與高階統(tǒng)計量的微地震初至拾取方法研究[J].石油物探,2015,54(4):388-395 SHENG G Q,LI Z C,WANG W B,et al.A new automatic detection method of microseismic events based on wavelet decomposition and high-order statistics[J].Geophysical Prospecting for Petroleum,2015,54(4):388-395
[9] 王小杰,欒錫武.基于小波分頻技術(shù)的地層Q值提取方法研究[J].石油物探,2015,54(3):260-266 WANG X J,LUAN X W.Q value extraction method based on wavelet frequency division technology[J].Geophysical Prospecting for Petroleum,2015,54(3):260-266
[10] BATTISTA B M,KNAPP C,MCGEE T,et al.Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data[J].Geophysics,2007,72(2):H29-H37
[11] 李月,彭蛟龍,馬海濤,等.過渡內(nèi)蘊模態(tài)函數(shù)對經(jīng)驗?zāi)B(tài)分解去噪結(jié)果的影響研究及改進算法[J].地球物理學(xué)報,2013,56(2):626-634 LI Y,PENG J L,MA H T,et al.Study of the influence of transition IMF on EMD de-noising and the improved algorithm[J].Chinese Journal of Geophysics,2013,56(2):626-634
[12] 董烈乾,李振春,楊少春,等.基于經(jīng)驗?zāi)B(tài)分解的f-x域面波壓制方法[J].石油地球物理勘探,2013,48(1):42-48 DONG L Q,LI Z C,YANG S C,et al.Ground roll suppression inf-xdomain based on empirical mode decomposition[J].Oil Geophysical Prospecting,2013,48(1):42-48
[13] HAN J,MIRKO V D B.Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding[J].Geophysics,2015,80(6):KS69-KS80
[14] 薛雅娟,曹俊興.聚合經(jīng)驗?zāi)B(tài)分解和小波變換相結(jié)合的地震信號衰減分析[J].石油地球物理勘探,2016,51(6):1148-1155 XUE Y J,CAO J X.Seismic attenuation analysis using ensemble empirical mode decomposition and wavelet transform[J].Oil Geophysical Prospecting,2016,51(6):1148-1155
[15] WU Z,HUANG N E.A study of the characteristics of white noise using the empirical mode decomposition method[J].Proceedings Mathematical Physical & Engineering Sciences,2004,460(2046):1597-1611
[16] YEH J R,SHIEH J S,HUANG N E.Complementary ensemble empirical mode decomposition:a novel noise enhanced data analysis method[J].Advances in Adaptive Data Analysis,2010,2(2):135-156
[17] 賈瑞生,趙同彬,孫紅梅,等.基于經(jīng)驗?zāi)B(tài)分解及獨立成分分析的微震信號降噪方法[J].地球物理學(xué)報,2015,58(3):1013-1023 JIA R S,ZHAO T B,SUN H M,et al.Micro-seismic signal denoising method based on empirical mode decomposition and independent component analysis[J].Chinese Journal of Geophysics,2015,58(3):1013-1023
[18] 梁喆,彭蘇萍,鄭晶.基于EMD和互信息熵的微震信號自適應(yīng)去噪[J].計算機工程與應(yīng)用,2014,50(4):7-11 LIANG Z,PENG S P,ZHENG J.Self-adaptive denoising for micro seismic signal based on EMD and mutual information entropy[J].Computer Engineering and Applications,2014,50(4):7-11
[19] SHANNON C E.A mathematical theory of communication[J].The Bell System Technical Journal,1948,27(4):379-423
[20] 李輝,戴旭初,葛洪魁,等.基于互信息量的地震信號檢測和初至提取方法[J].地球物理學(xué)報,2007,50(4):1190-1197 LI H,DAI X C,GE H K,et al.Seismic signal detection and first arrival pickup based on mutual information[J].Chinese Journal of Geophysics,2007,50(4):1190-1197
[21] 秦晅,宋維琪.基于時窗能量比與互信息量的微地震初至拾取方法[J].物探與化探,2016,40(2):374-379 QIN X,SONG W Q.Automatic first arrival pickup method of microseismic event based on energy ratio and mutual information[J].Geophysical and Geochemical Exploration,2016,40(2):374-379
[22] OMITAOMU O A,PROTOPOPESCU V A,GANGULY A R.Empirical mode decomposition technique with conditional mutual information for denoising operational sensor data[J].IEEE Sensors Journal,2011,11(10):2565-2575
[23] 尚帥,韓立國,胡瑋,等.壓縮小波變換地震譜分解方法應(yīng)用研究[J].石油物探,2015,54(1):51-55 SHANG S,HAN L G,HU W,et al.Applied research of synchrosqueezing wavelet transform in seismic spectral decomposition method[J].Geophysical Prospecting for Petroleum,2015,54(1):51-55
[24] HERRERA R H,HAN J,VAN DER BAANM.Applications of the synchrosqueezing transform in seismic time-frequency analysis[J].Geophysics,2014,79(3):V55-V64
[26] HERRERA R H,TARY J B,VAN DER BAAN M,et al.Body wave separation in the time-frequency domain[J].IEEE Geoscience and Remote Sensing Letters,2015,12(2):364-368
[27] 秦晅,宋維琪.基于同步壓縮變換微地震弱信號提取方法研究[J].石油物探,2016,55(1):60-66 QIN X,SONG W Q.Weak signal extraction method of microseismic data based on synchrosqueezing transform[J].Geophysical Prospecting for Petroleum,2016,55(1):60-66
[28] 汪祥莉,王斌,王文波,等.混沌干擾中基于同步擠壓小波變換的諧波信號提取方法[J].物理學(xué)報,2015,64(10):11-20 WANG X L,WANG B,WANG W B,et al.Harmonic signal extraction from chaotic interference based on synchrosqueezed wavelet transform[J].Acta Physica Sinica,2015,64(10):11-20
[29] 王姣,李振春,王德營.基于CEEMD的地震數(shù)據(jù)小波閾值去噪方法研究[J].石油物探,2014,53(2):164-172 WANG J,LI Z C,WANG D Y.A method for wavelet threshold denoising of seismic data based on CEEMD[J].Geophysical Prospecting for Petroleum,2014,53(2):164-172
(編輯:戴春秋)
MicroseismicdatadenoisingmethodbasedonEMDmutualinformationentropyandsynchrosqueezingtransform
QIN Xuan,CAI Jianchao,LIU Shaoyong,BIAN Aifei
(HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China)
On the basis of the characteristics of randomness,non-stationarity,time-frequency coupling of microseismic data,and of the problem of modal aliasing in empirical mode decomposition (EMD),this paper proposes a microseismic data denoising method based on EMD mutual information entropy and synchrosqueezing transform (SST).First,the microseismic signal is decomposed by EMD to acquire the intrinsic mode function (IMF) sequencing from high to low frequency.Next,the mutual information entropy of adjacent IMF components is calculated to identify the boundary between the high-frequency and the low-frequency part.Finally,the effective signal of the high-frequency part is extracted by the SST and reconstructed with the low-frequency part to achieve an effective microseismic data denoising.We applied the method to synthetic data sets with different noise intensities and to field data,and the results showed that this method can better remove the aliasing noise,extract the effective signal,and improve the SNR,compared with denoising methods that directly discard the high-frequency components.
microseismic,empirical mode decomposition,synchrosqueezing transform,mutual information entropy,reconstruction,denoising
P631
:A
1000-1441(2017)05-0658-09DOI:10.3969/j.issn.1000-1441.2017.05.006
秦晅,蔡建超,劉少勇,等.基于經(jīng)驗?zāi)B(tài)分解互信息熵與同步壓縮變換的微地震信號去噪方法研究[J].石油物探,2017,56(5):666
QIN Xuan,CAI Jianchao,LIU Shaoyong,et al.Microseismic data denoising method based on EMD mutual information entropy and synchrosqueezing transform
[J].Geophysical Prospecting for Petroleum,2017,56(5):666
2017-01-09;改回日期:2017-03-14。
秦晅(1988—),男,博士在讀,主要研究方向為微觀巖石物理與信號處理。
蔡建超(1981—),男,教授,博士生導(dǎo)師,主要從事微觀巖石物理研究。
國家重點研發(fā)計劃(2016YFC060110304)、國家自然科學(xué)基金(41572116)、中央高?;究蒲袠I(yè)務(wù)費專項資金(CUG160602)聯(lián)合資助。
This research is financially supported by the State Key Program of National Natural Science of China (Grant No.2016YFC060110304),the National Nature Science Foundation of China (Grant No.41572116) and the Fundamental Research Funds for the Central Universities (Grant No.CUG160602).