張 威,韓立國(guó),李洪建

(1.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,吉林長(zhǎng)春130026;2.中國(guó)石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

基于起伏海水表面的拖纜鬼波壓制方法

張 威1,韓立國(guó)1,李洪建2

(1.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,吉林長(zhǎng)春130026;2.中國(guó)石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

基于平坦海水表面假設(shè)的鬼波壓制方法在起伏海水表面情況下鬼波壓制效果不佳。為此,基于最小平方線性拉東反演去鬼波算法,反演得到海水表面處上行波場(chǎng),正向和反向延拓海水表面處上行波場(chǎng)至檢波器處,采用最小平方成像條件成像海水表面,以校正檢波器實(shí)際沉放深度,消除起伏海水表面的影響,實(shí)現(xiàn)了基于起伏海水表面的鬼波壓制。該方法考慮了起伏海水表面對(duì)鬼波壓制的影響,改善了基于平坦海水表面假設(shè)的方法壓制鬼波不完全的情況。數(shù)值模擬數(shù)據(jù)和海上實(shí)際拖纜數(shù)據(jù)測(cè)試結(jié)果表明,該方法能更好地壓制鬼波,最大程度拓寬地震記錄頻帶。

鬼波壓制;最小平方成像條件;起伏海水表面;最小平方線性拉東反演;拓寬頻帶

海洋拖纜地震數(shù)據(jù)采集時(shí),檢波器接收到的記錄包括來(lái)自地下反射界面的一次反射波和一次反射波經(jīng)過(guò)海水表面反射到達(dá)檢波器的虛反射(鬼波)[1]。鬼波緊隨一次波之后,與一次波相互干涉,降低了地震剖面的分辨率。同時(shí),鬼波的陷頻特性使得地震資料頻帶變窄,給地震資料反演和解釋帶來(lái)了困難[2-3]。

為了獲得寬頻帶的地震數(shù)據(jù),人們提出了各種各類的寬頻采集方式[4]。上下纜采集技術(shù)[5-7]通過(guò)對(duì)上下纜數(shù)據(jù)進(jìn)行合并,達(dá)到去鬼波的目的。雙檢采集技術(shù)[8]通過(guò)對(duì)壓力和速度傳感器采集的記錄進(jìn)行波場(chǎng)分離,得到上行波記錄。變深度纜采集技術(shù)[9-10]利用常規(guī)偏移和鏡像偏移聯(lián)合反褶積方法來(lái)壓制鬼波。但是,這些寬頻采集技術(shù)相對(duì)于常規(guī)水平拖纜采集技術(shù)而言,采集成本較高,并且需要特定的去鬼波方法。因此,基于常規(guī)水平拖纜采集技術(shù)的鬼波壓制方法仍然是研究的熱點(diǎn)。JOVANOVICH等[11]利用確定性反褶積壓制鬼波,其壓制效果很大程度上取決于鬼波延遲時(shí)間估計(jì)的精確性,但精確的鬼波延遲時(shí)間往往難以獲取,因此該方法并不能有效壓制鬼波。WANG等[12-13]基于bootstrap方法,分別在f-xy域和τ-p域?qū)崿F(xiàn)了去鬼波,該方法需要通過(guò)反演得到最優(yōu)鬼波延遲時(shí)間,但為了保證反演的穩(wěn)定性,需要進(jìn)行加窗處理,增加了計(jì)算的復(fù)雜性。POOLE[14]結(jié)合平面波傳播理論,對(duì)最小平方線性拉東方程進(jìn)行改進(jìn),在重建基準(zhǔn)面的同時(shí)壓制鬼波,該方法利用最小平方反演直接分離出上行波場(chǎng),在提高計(jì)算效率的同時(shí)提升了鬼波的壓制效果。在國(guó)內(nèi),王芳芳等[15]通過(guò)構(gòu)造頻率-波數(shù)域鬼波壓制算子,實(shí)現(xiàn)了數(shù)據(jù)驅(qū)動(dòng)的鬼波壓制;管西竹等[16]采用高效率的Fourier變換波場(chǎng)解析延拓方法,有效解決了水平上下纜鬼波壓制中有效信號(hào)畸變的問(wèn)題;李洪建等[17]基于格林函數(shù)理論,通過(guò)對(duì)壓力波場(chǎng)和垂直速度波場(chǎng)的預(yù)測(cè)壓制檢波點(diǎn)鬼波;ZHANG等[18]通過(guò)考慮實(shí)際粗糙海面反射系數(shù)對(duì)鬼波壓制的影響,利用高斯統(tǒng)計(jì)公式精確估計(jì)反射系數(shù),從而得到最優(yōu)化鬼波壓制算子。以上這些去鬼波算法均基于平坦海水表面的假設(shè),而實(shí)際的海水表面是起伏變化的,忽略起伏海水表面的影響使得鬼波延遲時(shí)間和海水表面反射系數(shù)的估計(jì)存在誤差,導(dǎo)致鬼波壓制不完全。

近年來(lái),ROBERTSSON等[19]利用壓力波場(chǎng)信息和近似的壓力波場(chǎng)梯度對(duì)海水表面進(jìn)行成像,實(shí)現(xiàn)起伏海水表面去鬼波,但是該方法對(duì)檢波器沉放深度和空間采樣間隔有嚴(yán)格的限制。KRAGH等[20]從壓力檢波器中提取低頻信息實(shí)現(xiàn)起伏海水表面的成像,成像質(zhì)量很大程度上取決于低頻信息提取的準(zhǔn)確性。ORJI等[21]將雙檢數(shù)據(jù)分解為上行波場(chǎng)和下行波場(chǎng)并分別延拓到海水表面,再對(duì)海水表面進(jìn)行成像,該方法受壓力和速度傳感器各自不同的靈敏度影響較大[22]。

本文基于POOLE[14]提出的最小平方線性拉東反演去鬼波方法,反演求解線性拉東方程得到海水表面處上行波波場(chǎng),將海水表面處上行波場(chǎng)正向和反向延拓至檢波器處,利用最小平方成像條件對(duì)海水表面進(jìn)行成像,以校正檢波器實(shí)際沉放深度;利用校正后的深度進(jìn)行去鬼波處理,并與POOLE提出的去鬼波方法進(jìn)行對(duì)比。合成數(shù)據(jù)和實(shí)際海上水平拖纜采集數(shù)據(jù)處理結(jié)果表明,基于起伏海水表面的鬼波壓制方法能更好地壓制鬼波,達(dá)到拓寬地震記錄頻帶的目的。

1 方法原理

1.1 海水表面處上行波波場(chǎng)的求取

根據(jù)頻率域最小平方線性拉東變換原理[23],檢波器接收到的總波場(chǎng)道集D(x,ω)的頻率拉東反變換可表示為:

(1)

式中:V(p,ω)為D(x,ω)拉東正變換的結(jié)果;ω為角頻率;p為慢度參數(shù);x為偏移距。

對(duì)于每一個(gè)頻率分量ω,可將公式(1)表示為矩陣形式:

(2)

式中:L=e-jωτn,m;τn,m=pm·xn;xn為炮集中第n道的偏移距;pm為拉東域第m個(gè)慢度道所對(duì)應(yīng)的慢度參數(shù)。模型V(p,ω)可通過(guò)最小平方反演求解線性方程(2)得到。

對(duì)于海洋水平拖纜采集,二維平面波傳播如圖1所示,第n道zn處的檢波器接收到總波場(chǎng),可以看作海水表面xn+Δxn,m處上行波場(chǎng)回傳到檢波器的波場(chǎng)與海水表面xn-Δxn,m處上行波場(chǎng)經(jīng)海水表面反射到達(dá)檢波器的下行波場(chǎng)之和,可將(2)式中延遲時(shí)間τn,m分解為τu和τg兩項(xiàng)[14],其表達(dá)式如下:

(3)

(3)式中Δxn,m,Δτn,m的幾何關(guān)系如下:

圖1 二維平面波傳播示意

(4)

式中:vw為海水速度;θm為出射角;zn為第n道檢波器的沉放深度。

結(jié)合海水表面處上行波場(chǎng)和檢波器處總波場(chǎng)之間的關(guān)系,通過(guò)對(duì)頻率域最小平方線性拉東方程進(jìn)行改造,得到新的線性方程,表達(dá)式如下:

(5)

(6)

式中:R為海水表面反射系數(shù);Lud為拉東反變換算子;Lu=e-jωτu為上行波拉東反變換算子;Ld=e-jωτg為鬼波拉東反變換算子;V0(pm)為海水表面處上行波波場(chǎng)拉東變換的結(jié)果。

1.2 起伏海水表面成像

對(duì)于水平拖纜采集數(shù)據(jù),通過(guò)求解公式(5)得到拉東域海水表面上行波波場(chǎng)V0(p,ω),利用頻率域拉東反變換得到時(shí)間-空間域海水表面上行波波場(chǎng)U0(x,t)。

(7)

(8)

(9)

(10)

在實(shí)際資料去鬼波處理中,海水表面的起伏形態(tài)隨時(shí)間變化,可通過(guò)選定時(shí)窗,使得時(shí)窗內(nèi)海水表面影響一致或者相似,準(zhǔn)確估計(jì)海水表面反射系數(shù)和檢波器實(shí)際沉放深度,然后對(duì)時(shí)窗內(nèi)的數(shù)據(jù)進(jìn)行去鬼波處理。

2 數(shù)值模擬算例

為了驗(yàn)證基于起伏海水表面的鬼波壓制方法的鬼波壓制效果,對(duì)層狀模型(圖2)的水平拖纜單炮記錄進(jìn)行去鬼波處理。正演模擬時(shí),采用主頻為80Hz的雷克子波,記錄道數(shù)為61,道間距為4m,采樣間隔為0.5ms,記錄長(zhǎng)度為0.3s,拖纜沉放深度為15m。為了驗(yàn)證本文方法的抗噪性,加入了信噪比為15的隨機(jī)噪聲。

圖2 水平層狀速度模型

模擬平坦海水表面和起伏海水表面下的單炮記錄分別如圖3a和圖3b所示,圖3c為圖3a和圖3b的差剖面?？梢?jiàn),由于受到起伏海水表面的影響,海水表面反射的鬼波同相軸不再是標(biāo)準(zhǔn)的雙曲線,呈現(xiàn)出與起伏海水表面波動(dòng)吻合的彎曲(圖3中紅色方框所示)。對(duì)比圖3a和圖3b中各道的旅行時(shí)差發(fā)現(xiàn)起伏海面引起的鬼波同相軸最大旅行時(shí)差為1.5ms,與所設(shè)計(jì)的起伏海水表面完全吻合。相反,一次反射波經(jīng)過(guò)海底反射直接傳播到檢波器處,并不受起伏海水表面的影響,其同相軸符合標(biāo)準(zhǔn)的雙曲線。

圖4a,圖4b為利用POOLE提出的方法對(duì)平坦海水表面和起伏海水表面下模擬的單炮記錄進(jìn)行鬼波壓制的結(jié)果?？梢?jiàn),POOLE提出的方法能完全壓制平坦海水表面下模擬單炮記錄中的鬼波;但對(duì)于起伏海水表面,由于該方法忽略了起伏海水表面的影響,鬼波壓制結(jié)果中存在殘余的鬼波(圖4b中藍(lán)色箭頭所指處),同時(shí)會(huì)引入較強(qiáng)的竄擾噪聲(圖4b中紅色箭頭所指處)。

圖3 平坦海水表面下模擬的單炮記錄(a)、起伏海水表面下模擬的單炮記錄(b)以及兩者的差剖面(c)

圖4 POOLE方法對(duì)平坦海水表面(a)、起伏海水表面(b)下模擬單炮記錄的去鬼波結(jié)果以及本文方法對(duì)起伏海水表面下模擬單炮記錄的去鬼波結(jié)果(c)

利用本文方法對(duì)起伏海水表面條件下模擬的單炮記錄(圖3b)的海水表面形態(tài)進(jìn)行成像,結(jié)果如圖5中藍(lán)色曲線所示,黑色曲線表示平坦海水表面,紅色曲線表示模擬的起伏海水表面?？梢?jiàn),成像的海水表面與模擬的起伏海水表面在形態(tài)上匹配良好,部分差異主要是由主頻80Hz子波的分辨率限制所引起的。

基于成像的海水表面對(duì)起伏海水表面下模擬的單炮記錄進(jìn)行鬼波壓制,得到的去鬼波結(jié)果如圖4c所示。與POOLE方法的去鬼波結(jié)果(圖4b)相比,由于考慮了起伏海水表面對(duì)鬼波的影響,殘余鬼波得到了壓制,竄擾噪聲也隨之減少,壓制效果得到明顯提升。但去鬼波結(jié)果中仍存在少量的竄擾噪聲,這是由海水表面成像誤差所導(dǎo)致的。

圖5 平坦海水表面、成像海水表面以及模擬海水表面

對(duì)去鬼波前后記錄進(jìn)行頻譜分析,如圖6所示,POOLE方法(藍(lán)色曲線)和本文方法(紅色曲線)去鬼波后的振幅譜中陷波現(xiàn)象得到消除,能量得到補(bǔ)償,頻帶得到拓寬。本文方法對(duì)于低頻端能量提升更加明顯(綠色方框所示)。

圖6 去鬼波前后振幅譜對(duì)比

3 實(shí)際水平拖纜數(shù)據(jù)測(cè)試

利用海上某工區(qū)水平拖纜采集資料進(jìn)一步驗(yàn)證本文方法的適用性(圖7),其采樣間隔2ms,道間距12.5m,拖纜沉放深度為25m,采樣總時(shí)間為10.572s,震源沉放深度為15m。圖7a為經(jīng)過(guò)直達(dá)波切除、濾波去噪、多次波壓制等預(yù)處理后的原始單炮記錄(選取2.6～3.4s,10～190道)。首先,利用本文方法對(duì)實(shí)際單炮記錄的海水表面形態(tài)進(jìn)行成像,結(jié)果如圖8所示。然后,基于成像的海水表面進(jìn)行鬼波壓制,得到鬼波壓制后的單炮記錄,如圖7c所示。與由POOLE方法壓制鬼波后的單炮記錄(圖7b)進(jìn)行對(duì)比可見(jiàn),2種方法壓制鬼波后的單炮記錄中反射同相軸相位個(gè)數(shù)減少,炮記錄變得更加清晰,分辨率得到提高。從兩種方法差剖面(圖7d)中可以看出,本文方法壓制鬼波的記錄中殘余鬼波更少,壓制效果更佳。

圖7 原始單炮記錄(a)、POOLE方法壓制鬼波后的單炮記錄(b)、本文方法壓制鬼波后的單炮記錄(c)以及兩種方法處理結(jié)果的差剖面(d)

圖8 成像海水表面

從成像的海水表面形態(tài)上看,30～80道和120～170道受海水表面起伏影響較大,兩種方法在相應(yīng)區(qū)域的鬼波壓制效果存在微小的差異(由圖7d差剖面可以明顯看出)。為了驗(yàn)證本文方法的優(yōu)越性,選取去鬼波前后相應(yīng)的區(qū)域進(jìn)行頻譜分析,如圖9a和圖9b所示。在去鬼波前的炮記錄中,在30,50,60,90Hz處存在明顯的陷波現(xiàn)象,經(jīng)過(guò)鬼波壓制,2種方法都能有效消除檢波器鬼波所引起的陷頻特征,拓寬地震記錄的頻帶。其中50Hz為震源鬼波引起的陷波點(diǎn),因此,能量變化不大。與POOLE方法相比,本文方法無(wú)論是對(duì)低頻能量還是高頻能量的提升都具有明顯優(yōu)勢(shì)(圖9黑色箭頭所指處)。

圖9 30～80道去鬼波前后振幅譜對(duì)比(a)以及120～170道去鬼波前后振幅譜對(duì)比(b)

4 結(jié)論

本文通過(guò)對(duì)起伏海水表面進(jìn)行成像,實(shí)現(xiàn)了基于起伏海水表面的拖纜鬼波壓制。通過(guò)對(duì)模擬單炮數(shù)據(jù)和實(shí)際海上數(shù)據(jù)處理結(jié)果的分析,可得以下結(jié)論和認(rèn)識(shí):

1) 基于起伏海水表面的拖纜鬼波壓制方法能有效地壓制拖纜采集數(shù)據(jù)中的鬼波,恢復(fù)陷波點(diǎn)能量,提高了地震資料的分辨率,拓寬了地震資料的頻帶。

2) 相對(duì)于基于平坦海水表面的鬼波壓制方法,基于起伏海水表面的拖纜鬼波壓制方法考慮了起伏海水表面對(duì)鬼波的影響,鬼波壓制效果更佳,低頻和高頻能量的提升更加明顯,能夠適應(yīng)實(shí)際起伏海水表面情況。

3) 起伏海水表面的拖纜鬼波壓制方法基于最小平方成像條件實(shí)現(xiàn)了對(duì)海水表面的成像,但起伏海水表面的成像精度在一定程度上取決于地震數(shù)據(jù)的信噪比。

[1] HAMARBITAN N S,MARGRAVE G F.Spectral analysis of a ghost[J].Geophysics,2001,66(4):1267-1273

[2] 李套山,高洪強(qiáng),劉振夏,等.虛反射信息的采集及其形成機(jī)制、頻率響應(yīng)的理論探討[J].石油物探,1997,36(4):38-44 LI T S,GAO H Q,LIU Z X,et al.The acquisition of ghost information and the study of its generation and frequency response[J].Geophysical Prospecting for Petroleum,1997,36(4):38-44

[3] 李洪建,韓立國(guó),鞏向博.頻率域起伏雙纜數(shù)據(jù)合并策略[J].石油物探,2015,54(6):674-680 LI H J,HAN L G,GONG X B.Merge strategy for fluctuate double streamers in frequency domain[J].Geophysical Prospecting for Petroleum,2015,54(6):674-680

[4] 吳志強(qiáng).海洋寬頻帶地震勘探技術(shù)新進(jìn)展[J].石油地球物理勘探,2014,49(3):421-430 WU Z Q.New advances in marine broadband seismic exploration[J].Oil Geophysical Prospecting,2014,49(3):421-430

[5] POSTHUMUS B J.Deghosting of twin streamer configuration[J].Geophysical Prospecting,1993,41(3):267-286

[6] ?ZDEMIR A K,CAPRIOLI P,?ZBEK A,et al.Optimized deghosting of over/under towed-streamer data in the presence of noise[J].The Leading Edge,2008,27(2):190-199

[7] 劉春成,劉志斌,顧漢明.利用上/下纜合并算子確定海上上/下纜采集的最優(yōu)沉放深度組合[J].石油物探,2013,52(6):623-629 LIU C C,LIU Z B,GU H M.The determination of optimal sinking depths of over/under streamers in offshore survey by merge operator[J].Geophysical Prospecting for Petroleum,2013,52(6):623-629

[8] CARLSON D,S?LLNER W,TABTI H,et al.Increased resolution of seismic data from a dual-sensor streamer cable[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007,994-998

[9] SOUBARAS R.Deghosting by joint deconvolution of a migration and a mirror migration[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:3406-3410

[10] 許自強(qiáng),方中于,顧漢明,等.海上變深度纜數(shù)據(jù)最優(yōu)化壓制鬼波方法及其應(yīng)用[J].石油物探,2015,54(4):404-413 XU Z Q,FANG Z Y,GU H M,et al.The application of optimal deghosting algorithm on marine variable-depth streamer data[J].Geophysical Prospecting for Petroleum,2015,54(4):404-413

[11] JOVANOVICH D B,SUMNER R D,AKINS-EASTERLIN S L.Ghosting and marine signature deconvolution:a prerequisite for detailed seismic interpretation[J].Geophysics,1983,48(11):1468-1485

[12] WANG P,PENG C.Premigration deghosting for marine towed streamer data using a bootstrap approach[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[13] WANG P,RAY S,PENG C,et al.Premigration deghosting for marine streamer data using a bootstrap approach in tau-p domain[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:4221-4225

[14] POOLE G.Premigration receiver deghosting and redatuming for variable depth streamer data[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:4216-4220

[15] 王芳芳,李景葉,陳小宏.基于逆散射級(jí)數(shù)法的鬼波壓制方法[J].地球物理學(xué)報(bào),2013,56(5):1628-1636 WANG F F,LI J Y,CHEN X H.Deghosting method based on inverse scattering series[J].Chinese Journal of Geophysics,2013,56(5):1628-1636

[16] 管西竹,陳寶書(shū),符立耘,等.基于波動(dòng)方程的上下纜地震數(shù)據(jù)鬼波壓制方法研究[J].地球物理學(xué)報(bào),2015,58(10):3746-3757 GUAN X Z,CHEN B S,FU L Y,et al.The study of a deghosting method of over/under streamer seismic data based on wave equation[J].Chinese Journal of Geophysics,58(10):3746-3757

[17] 李洪建,韓立國(guó),鞏向博,等.基于格林函數(shù)理論的波場(chǎng)預(yù)測(cè)和鬼波壓制方法[J].地球物理學(xué)報(bào),2016,59(3):1113-1124 LI H J,HAN L G,GONG X B,et al.A wavefield prediction and ghost suppression method based on Green function theory[J].Chinese Journal of Geophysics,2016,59(3):1113-1124

[18] ZHANG X Y,PAN D M,SHI W Y,et al.Deghosting towed streamer data inτ-pdomain based on rough sea surface reflectivity[J].Applied Geophysics,2015,12(4):573-584

[19] ROBERTSSON J O A,KRAGH E.Rough-sea deghosting using a single streamer and a pressure gradient approximation[J].Geophysics,2002,67(6):2005-2011

[20] KRAGH E,LAWS R.Rough seas and statistical deconvolution[J].Geophysical Prospecting,2006,54(4):475-485

[21] ORJI O,S?LLNER W,GELIUS L J.Imaging the sea surface using a dual-sensor towed streamer[J].Geophysics,2010,75(6):V111-V118

[22] TENGHAMN R,VAAGE S,BORRESEN C.A dual-sensor towed marine streamer:its viable implementation and initial results[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:989-993

[23] KOSTOV C.Toeplitz structure in slant-stack inversion[J].Expanded Abstracts of 60thAnnual Internat SEG Mtg,1990:1618-1621

[24] SCHLEICHER J,COSTA J C,NOVAIS A.A comparison of imaging conditions for wave-equation shot-profile migration[J].Geophysics,2008,73(6):S219-S227

(編輯:朱文杰)

Deghosting method based on a variable sea surface for conventional streamer seismic data

ZHANG Wei1,HAN Liguo1,LI Hongjian2

(1.CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130026,China;2.SinopecGeophysicalResearchInstitute,Nanjing211103,China)

Currently,the deghosting algorithms based on a flat sea surface assumption can not provide a robust deghsoing result,caused by the misestimation of sea surface reflection and ghost delay time.A new method of deghsoting has been developed by solving the least-square linear Radon equation.The up-going wavefield at the sea surface can be obtained and then extrapolate the up-going wavefield forward and backward to the receivers where a least-square imaging condition is applied so that the variation of the sea surface can be estimated.The effect of the variable sea surface on the ghosts can be taken into account by imaging the sea surface.The sea surface is then provided to the correction of the actual depth of the receiver.Finally,we compare the deghosting results using a flat sea surface and a variable sea surface.Less ghost residuals and reduced ringing artifacts are observed when deghosting was performed using a variable sea surface.Test results on synthetic and field data show that the proposed approach can better remove ghost wave and achieve the goal of broadening the bandwidth of seismic data.

deghosting,least-square imaging condition,variable sea surface,least squares linear Radon inversion,broaden the bandwidth

2016-06-30;改回日期：2017-01-24。

張威(1992—),男,碩士在讀,現(xiàn)主要從事寬頻地震采集、處理方面的研究。

韓立國(guó)(1961—),男,教授,博士生導(dǎo)師,主要從事復(fù)雜地震波場(chǎng)正反演與成像、地震屬性與地球物理信號(hào)分析研究。

國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2014AA06A605)資助。

P631

A

1000-1441(2017)04-0500-07

10.3969/j.issn.1000-1441.2017.04.005

The research is financially supported by the Nationl High-tech R & D Program of China (863 Program) (Grant No.2014AA06A605).