王 杰,胡光輝,劉定進,邵文潮,王鵬燕

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

陸上地震資料全波形反演策略研究

王 杰,胡光輝,劉定進,邵文潮,王鵬燕

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

全波形反演已成功應(yīng)用于海上地震資料，而陸上地震資料全波形反演應(yīng)用還面臨諸多難點。研究給出了一種針對陸上地震資料的全波形反演策略：首先采用相位擬合互相關(guān)目標(biāo)泛函增強全波形反演的適用性；其次構(gòu)建動態(tài)波數(shù)域濾波梯度預(yù)處理方法壓制梯度噪聲；最后采用構(gòu)造傾角信息約束的自適應(yīng)高斯平滑算子改善反演結(jié)果質(zhì)量。二維陸上地震資料的反演結(jié)果表明，該策略可有效實現(xiàn)陸上地震資料高波數(shù)全波形反演速度建模。

全波形反演;互相關(guān)目標(biāo)泛函;分區(qū)波數(shù)域濾波;自適應(yīng)各向異性高斯平滑

全波形反演(FWI)具有求取高分辨率、高精度地下物性參數(shù)的能力。TARANTOLA[1]于20世紀80年代首次提出了基于廣義最小二乘的時間域全波形反演技術(shù)。近10年來,隨著計算機技術(shù)的發(fā)展,全波形反演成為國內(nèi)外地震勘探領(lǐng)域研究的熱點之一。目前全波形反演理論研究日臻成熟,但相較于線性射線層析,全波形反演由于匹配走時、振幅和相位等信息,非線性程度高,極易陷入局部極值產(chǎn)生周期跳問題,此外,所需計算資源巨大,因此,全波形反演的工業(yè)化應(yīng)用面臨諸多挑戰(zhàn)。

傳統(tǒng)全波形反演是一高度非線性反問題,數(shù)據(jù)缺失低頻信息或初始模型不準時,模型在迭代更新過程中容易陷入局部極小值[2]。為降低反演的非線性性,BUNKS等[3]、SIRGUE等[4]分別從時間空間域和頻率空間域給出從低頻到高頻的多尺度反演策略。BROSSIER等[5]采用凸性更好、性態(tài)更優(yōu)的互相關(guān)目標(biāo)泛函完成全波形反演,相較常規(guī)L2范數(shù)而言,其局部極小值更少。WU等[6]從原始地震數(shù)據(jù)中提取包含豐富低頻信息的信號包絡(luò)反演宏觀背景速度以解決低頻缺失造成的局部極小值問題。LI等[7]采用相位追蹤和頻率外推方法從地震波高頻數(shù)據(jù)中人工模擬出低頻信息并完成了低波數(shù)背景速度重構(gòu),削弱了傳統(tǒng)全波形反演的非凸性。

巨大的計算成本是全波形反演在實際應(yīng)用中的重要制約因素。為提高全波形反演計算效率,KREBS等[8]和BEN-HADJ-ALI等[9]分別在時間域和頻率域采用超級震源技術(shù)顯著減少每次迭代正演模擬的計算時間,并采用隨機相位編碼有效壓制反演結(jié)果中的串?dāng)_噪聲,在保證反演質(zhì)量的前提下大大提高了全波形反演計算效率。

目前國際上海洋地震資料全波形反演成功應(yīng)用實例越來越多。在SIRGUE等[10]采用三維寬方位海底電纜數(shù)據(jù)對北海油田成功構(gòu)建全波形反演高精度速度模型之后,大偏移距、寬方位、寬頻海上拖纜數(shù)據(jù)、海底電纜數(shù)據(jù)、常規(guī)拖纜采集數(shù)據(jù)及寬頻數(shù)據(jù)、窄方位海底電纜數(shù)據(jù)等淺水環(huán)境和深水環(huán)境全波形反演成功應(yīng)用案例越來越多。

相較海上資料,陸上地震資料全波形反演的成功應(yīng)用還存在不小的挑戰(zhàn)[11]。在陸上地震數(shù)據(jù)采集過程中,受震源強度、震源激發(fā)響應(yīng)、表層吸收衰減、各種干擾噪聲、復(fù)雜近地表傳播機制及檢波器耦合效應(yīng)等影響,不同震源、不同檢波點地震信號振幅能量嚴重不均衡,局部地區(qū)存在震源稀疏,數(shù)據(jù)不完備現(xiàn)象,且受采集儀器限制,地震信號往往缺失低頻及大偏移距數(shù)據(jù)信息。PLESSIX等[12-13]采用專門為全波形反演采集的低頻、大偏移距、高密度陸地可控震源數(shù)據(jù)完成了內(nèi)蒙古地區(qū)二維高精度速度模型建模,并取得較好應(yīng)用效果。HE等[14]指出,當(dāng)陸上地震資料數(shù)據(jù)不完備時,反演結(jié)果中存在較強噪聲,迭代收斂性和穩(wěn)健性都有一定程度影響,為此他們給出了一種有效的噪聲算子濾波方法,然而當(dāng)?shù)刭|(zhì)條件較為復(fù)雜、構(gòu)造起伏劇烈時,該方法仍存在一定的局限性。我國以陸上探區(qū)為主,研發(fā)適合于陸上地震資料的全波形反演技術(shù)迫在眉睫。

本文給出適用于不完備陸上地震資料的穩(wěn)健的全波形反演策略:采用相位匹配的非零延遲歸一化互相關(guān)目標(biāo)泛函提高反演問題的凸性和穩(wěn)健性,降低不同震源和檢波器振幅能量不均衡對反演造成的影響;基于已知工區(qū)構(gòu)造先驗信息將反演工區(qū)分塊并對不同區(qū)塊在波數(shù)域完成反演噪聲壓制;對整體模型更新量采用自適應(yīng)高斯平滑技術(shù)減弱不同區(qū)塊模型更新量的耦合效應(yīng),從而使該反演策略適用于任意地質(zhì)構(gòu)造條件下陸上地震資料的高精度速度建模。最后采用東北某探區(qū)實際二維地震資料進行了應(yīng)用測試。

1 方法原理

1.1 互相關(guān)目標(biāo)泛函反演框架

L2范數(shù)全波形反演在迭代過程中需要與地震數(shù)據(jù)的振幅、波形、相位信息匹配,其對數(shù)據(jù)振幅信息非常敏感。由于陸上地震資料激發(fā)和接收條件等因素的限制,原始采集的地震記錄存在炮間、道間能量不均衡,在局部區(qū)域非常稀疏。當(dāng)采用此類數(shù)據(jù)進行凸性相對較差的L2范數(shù)全波形反演時,不均衡的殘差能量會引起模型更新量在不同地區(qū)存在較大差異,從而導(dǎo)致反演不收斂或不穩(wěn)定,使反演結(jié)果陷入局部極值。為提高反演的穩(wěn)健性,在進行陸上地震資料反演時需要給旅行時更多的權(quán)重,以降低振幅的影響?；谙辔粩M合的互相關(guān)目標(biāo)泛函更強調(diào)匹配數(shù)據(jù)的相位信息,凸性更優(yōu),且相對于L2范數(shù)具有更廣的波谷,因而反演結(jié)果更加穩(wěn)健[5,15]。本文采用BROSSIER等[5]給出的非零延遲歸一化互相關(guān)目標(biāo)泛函理論框架實現(xiàn)陸上地震資料的全波形反演。非零延遲歸一化互相關(guān)目標(biāo)泛函為:

(1)

式中:m為要反演的模型參數(shù);τ為延遲時間;xr為檢波器坐標(biāo);xs為震源坐標(biāo);W(τ)為懲罰函數(shù);dcal(t,xr;xs)為人工合成地震記錄;dobs(t+τ,xr;xs)為時延觀測地震記錄。

即使觀測地震記錄和模擬地震記錄相位差大于半個地震波波長,該目標(biāo)泛函亦可穩(wěn)定收斂,而且歸一化的引入也減小了不同震源和檢波器振幅能量不均衡引起的噪聲干擾問題。梯度表達式可用伴隨狀態(tài)法求解,但與L2范數(shù)不同,非零延遲歸一化互相關(guān)目標(biāo)泛函需修改反傳震源項[5]。采用目標(biāo)泛函對模型參數(shù)的導(dǎo)數(shù)可求解梯度表達式:

(2)

1.2 分區(qū)波數(shù)域濾波算子

NANGOO[16]采用稀疏震源子集對常規(guī)缺乏低頻和大偏移距的窄方位海上拖纜數(shù)據(jù)進行了全波形反演研究,采用水平光滑預(yù)處理技術(shù)壓制梯度大傾角噪聲,保留反演速度中的地質(zhì)構(gòu)造信息。當(dāng)工區(qū)地質(zhì)背景較為復(fù)雜時,水平光滑預(yù)處理技術(shù)不再適用。相對海上拖纜數(shù)據(jù),陸上資料的稀疏性更強,而且信噪比低,數(shù)據(jù)質(zhì)量相對更差。當(dāng)工區(qū)震源較為稀疏,且覆蓋次數(shù)較低時,全波形反演會在梯度中引入強烈的反演噪聲[14,16]。為得到較為準確的、符合地質(zhì)構(gòu)造背景的速度模型,需要壓制迭代過程中產(chǎn)生的反演噪聲。HE等[14]提出了在波數(shù)域壓制噪聲的方法,但該方法的不足之處在于若選擇的構(gòu)造傾角較小,則高陡構(gòu)造也被壓制;若選擇的構(gòu)造傾角過大,雖然在濾波后的反演結(jié)果中保留了高陡構(gòu)造特征,但會帶來保留較強噪聲干擾的結(jié)果。本文為突破上述單角度波數(shù)域濾波的限制,基于已有地質(zhì)構(gòu)造信息,對地下介質(zhì)分區(qū),采用分而治之思想,對不同區(qū)塊選取相應(yīng)角度的波數(shù)域濾波算子壓制噪聲,在有效壓制噪聲的前提下,盡可能保留地質(zhì)構(gòu)造特征信息。繼而采用整體自適應(yīng)高斯平滑濾波技術(shù)減小不同區(qū)塊梯度波數(shù)域濾波后的不一致性。

分區(qū)波數(shù)域濾波需要在每次迭代過程中按照成像剖面沿地質(zhì)構(gòu)造傾角分區(qū)。局部構(gòu)造傾角場可采用局部平面波解構(gòu)濾波法[17]或圖像梯度結(jié)構(gòu)張量算法[18]求解,本文采用局部平面波解構(gòu)法求解構(gòu)造傾角場。平面波解構(gòu)濾波器利用局部平面波的疊加表征地震數(shù)據(jù),可以較好地估計平滑連續(xù)同相軸的局部傾角信息。根據(jù)局部平面波的物理模型,用局部有限差分方程定義平面波解構(gòu)濾波器為:

(3)

式中:P(x,t)為局部平面波波場;σ(x,t)為平面波局部斜率場。(3)式可化簡為基于最小二乘的傾角估計

問題:

(4)

式中:σ0為斜率初始值;Δσ為斜率增量;C為局部傾角的函數(shù);d為原始數(shù)據(jù)。利用(4)式可構(gòu)建局部傾角場。根據(jù)構(gòu)造的局部傾角場將整個地質(zhì)工區(qū)分為不同區(qū)塊,每個區(qū)塊對應(yīng)不同的主要構(gòu)造傾角。利用不同區(qū)塊的主要構(gòu)造傾角和各自相應(yīng)的波數(shù)域濾波算子對不同區(qū)塊梯度壓制噪聲。各區(qū)塊采用的噪聲壓制算子公式為:

(5)

1.3 自適應(yīng)各向異性高斯平滑算子

分區(qū)濾波梯度場在整合為整體梯度時不同區(qū)塊交界處存在邊界效應(yīng),這將在反演結(jié)果中引入明顯的區(qū)塊邊界。為削弱這種邊界效應(yīng),本文對整體梯度場采用自適應(yīng)各向異性高斯平滑濾波算子[19]沿傾角方向濾波。各向異性高斯濾波算子在x,y平面上的投影為橢圓(如圖1a所示),該濾波算子表達式為:

(6)

針對不同的構(gòu)造傾角θ,采用旋轉(zhuǎn)的方式得到如圖1b所示的投影,由坐標(biāo)變換可得旋轉(zhuǎn)后的坐標(biāo)u-v和x-y坐標(biāo)的關(guān)系為:

(7)

旋轉(zhuǎn)后的高斯濾波算子為:

(8)

圖1 各向異性高斯平滑算子無旋轉(zhuǎn)(a)和旋轉(zhuǎn)θ(b)在x-y平面上的投影

采用構(gòu)建的自適應(yīng)各向異性高斯濾波算子可對整合后的梯度場沿著構(gòu)造傾角平滑:

(9)

式中:sn為整合梯度場;dn表示高斯濾波梯度場;G為自適應(yīng)各向異性高斯濾波算子。

在構(gòu)建出構(gòu)造較為連續(xù)、噪聲污染較少的梯度場后,采用預(yù)條件擬牛頓(L-BFGS)迭代最優(yōu)化算法對速度模型迭代更新,迭代步長采用非精確線性搜索技術(shù)求取。

2 實際資料應(yīng)用

為驗證本文方法的有效性,選取東北某探區(qū)一條二維線進行全波形反演。沿探區(qū)東西方向抽取44炮單炮地震記錄,每炮180道左右。震源在x方向的分布范圍為13000m,炮點距350m左右,相對稀疏。圖2為該二維線覆蓋次數(shù)圖。其中,最大覆蓋次數(shù)為16次。圖3a為去除面波、線性噪聲、異常振幅,并進行地表一致性振幅補償、球面擴散補償后的原始單炮地震記錄。圖3b為6Hz以上頻率段觀測地震記錄。該記錄中面波得到了有效壓制,信噪比得以提高。圖3c 為圖3b的低頻段信息,頻率范圍為6～20Hz。本文采用6～20Hz頻率段數(shù)據(jù)進行反演測試。圖4為二維深度域克?；舴虔B前偏移剖面。對偏移剖面采用局部平面波解構(gòu)濾波可獲取地質(zhì)構(gòu)造傾角信息(圖5)。從構(gòu)造傾角場中可以看出,淺層的主要反射層構(gòu)造傾角較小,大部分集中在10°之內(nèi),中、深層構(gòu)造傾角較大,為-25°～35°。

地震子波是傳統(tǒng)基于L2范數(shù)實際資料全波形反演成敗的一個關(guān)鍵因素[11]。由于實際資料的地震子波是時變和空變的,而且在相同條件下,不同子波正演模擬出的人工合成地震記錄在振幅、相位和波形上均存在較大差別,因而如何做好地震子波估計是全波形反演的重要環(huán)節(jié)。本文采用基于相位擬合的非零延遲歸一化互相關(guān)目標(biāo)反演理論框架完成速度場的重構(gòu),在反演迭代過程中削弱了原始地震記錄和觀測地震記錄振幅和波形的影響,更加強調(diào)相位的作用,在一定程度上減小了對地震子波估計的要求,因此本文采用合適主頻的雷克子波完成正演模擬及速度重構(gòu),非零延遲的引入在一定程度上允許地震子波具有一定的相位延遲。本文采用的初始速度模型為反射波旅行時層析速度反演模型(圖6),該模型具有較為充分的低波數(shù)信息及宏觀地質(zhì)構(gòu)造背景信息。

圖2 研究區(qū)某二維線覆蓋次數(shù)分布

圖3 研究區(qū)某二維線單炮地震記錄a 去噪和振幅補償后的觀測地震記錄; b 消除面波后的觀測地震記錄; c 6～20Hz觀測地震記錄

采用原始稀疏觀測地震記錄進行互相關(guān)全波形反演時,由于覆蓋次數(shù)較少,反演的非線性程度較高,反演結(jié)果存在較強的噪聲。圖7顯示了6～20Hz數(shù)據(jù)迭代20次的全波形反演結(jié)果。由圖7可見,反演結(jié)果信噪比較低,有效地質(zhì)構(gòu)造信息被反演噪聲湮沒。圖8給出了在每次迭代過程中采用單一角度(10°,15°,25°,35°)波數(shù)域濾波迭代20次得到的反演結(jié)果。相較于原始資料直接重構(gòu)的速度模型,采用單一角度波數(shù)域濾波后構(gòu)建的反演速度有效壓制了串?dāng)_噪聲,增強了有效地質(zhì)構(gòu)造的連續(xù)性。然而單一角度波數(shù)域濾波本身存在缺陷,當(dāng)采用的濾波算子角度過小時,雖然反演結(jié)果信噪比顯著提高,但濾波算子消除了大角度的地質(zhì)構(gòu)造信息(圖8a,圖8b)。由于該工區(qū)的構(gòu)造傾角范圍為-25°～35°,當(dāng)采用大角度波數(shù)域濾波算子時,雖然大傾角的地質(zhì)構(gòu)造信息得以保留,但是反演結(jié)果中也相應(yīng)引入了較為嚴重的大角度線性噪聲(圖8c,圖8d)。

圖4 深度域克?；舴虔B前偏移剖面

圖5 疊前偏移剖面構(gòu)造傾角場

圖6 反射波旅行時層析速度模型

圖7 互相關(guān)目標(biāo)泛函全波形反演結(jié)果

該工區(qū)構(gòu)造傾角在淺層集中在10°以下,深層主要集中在25°左右,因此本文采用淺、深層分區(qū)波數(shù)域濾波算子完成全波形反演建模。圖9給出了分區(qū)波數(shù)域濾波建模結(jié)果。由圖9可見,淺層和深層都保留了有效地質(zhì)構(gòu)造信息,與單一大角度波數(shù)域濾波反演結(jié)果(圖8c,圖8d)相比,分辨率有一定改善,但是在紅色箭頭處,出現(xiàn)低速異常體,反演結(jié)果陷入局部極值,且白色箭頭的位置速度出現(xiàn)不耦合現(xiàn)象。為此,每次迭代在分區(qū)波數(shù)域濾波的基礎(chǔ)上引入各向異性自適應(yīng)高斯平滑算子。圖10顯示了最終的重構(gòu)模型,該重構(gòu)結(jié)果不僅信噪比提高而且構(gòu)造更加連續(xù),局部極值問題得到削弱。

圖8 采用單一角度波數(shù)域濾波后互相關(guān)目標(biāo)泛函全波形反演結(jié)果a 濾波角度10°; b 濾波角度15°; c 濾波角度25°; d 濾波角度35°

圖11給出了初始模型、反演速度模型正演模擬歸一化地震記錄和歸一化觀測地震記錄的對比結(jié)果。

圖9 分區(qū)波數(shù)域濾波互相關(guān)全波形反演結(jié)果(淺層10°,深層15°)

圖10 分區(qū)波數(shù)域濾波后采用各向異性自適應(yīng)高斯平滑處理得到的互相關(guān)全波形反演結(jié)果

圖11 觀測地震記錄與模擬地震記錄結(jié)果對比a 初始模型正演模擬地震記錄(左)與觀測地震記錄(右); b 反演模型正演模擬地震記錄(左)與觀測地震記錄(右)

與初始模型正演模擬地震記錄相比,反演模型的模擬地震記錄主要反射層的相位與觀測地震記錄匹配得更好,從相位擬合的角度驗證了反演的準確性與穩(wěn)健性。

3 結(jié)論

本文給出陸上地震資料不完備條件下采用全波形反演進行高精度速度建模的方案。采用基于相位擬合的非零延遲歸一化互相關(guān)目標(biāo)泛函取代傳統(tǒng)的基于數(shù)據(jù)波形、振幅、相位匹配的L2范數(shù)目標(biāo)泛函,提高了反演的穩(wěn)健性,且降低了對地震子波精度和原始地震資料振幅能量均衡的要求。采用分區(qū)波數(shù)域濾波算子和各向異性自適應(yīng)高斯平滑技術(shù)有效減少了原始地震資料稀疏時反演結(jié)果中較強的反演噪聲并保留了有效地質(zhì)構(gòu)造信息,反演結(jié)果更加合理、準確,而且反演過程較為穩(wěn)健。二維陸上地震資料反演測試結(jié)果表明本文給出的陸上地震資料反演策略不僅有效消除了反演結(jié)果中強烈的噪聲,而且構(gòu)建了具有地質(zhì)意義的高分辨率速度模型,波形的匹配表明本文反演方法相對層析反演結(jié)果更優(yōu)。

[1] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[2] VIRIUEX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[3] BUNKS C,SALECK F.Multi-scale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

[4] SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(1):231-248

[5] BROSSIER R,OPERTO S.Velocity model building from seismic reflection data by full-waveform inversion[J].Geophysical Prospecting,2015,63(2):354-367

[6] WU R S,LUO J,WU B.Seismic envelope inversion and modulation signal model[J].Geophysics,2014,79(3):WA13-WA24

[7] LI Y Y,DEMANET L.Full-waveform inversion with extrapolated low-frequency data[J].Geophysics,2016,81(6):R339-R348

[8] KREBS J R,ANDERSON J E,HINKLEY D,et al.Fast full-wavefield seismic inversion using encoded sources[J].Geophysics,2009,74(6):WCC177-WCC188

[9] BEN-HADJ-ALI H,OPERTO S,VIRIEUX J.An efficient frequency-domain full waveform inversion method using simultaneous encoded sources[J].Geophysics,2011,76(4):R109-R124

[10] SIRGUE L,BARKVEL O.3D waveform inversion in Valhall wide-azimuth OBC[J].Expanded Abstracts of 71stAnnual Internat EAGE Mtg,2009:U038

[11] 楊勤勇,胡光輝,王立歆.全波形反演研究現(xiàn)狀及發(fā)展趨勢[J].石油物探,2014,53(1):77-83 YANG Q Y,HU G H,WANG L X.Research status and development trend of full waveform inversion[J].Geophysical Prospecting for Petroleum,2014,53(1):77-83

[12] PLESSIX R E,BAETEN G.Application of acoustic full waveform inversion to a low-frequency large-offset land dataset[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:930-934

[13] PLESSIX R E,BAETEN G.Full waveform inversion and distance separated simultaneous sweeping:a study with a land seismic data set[J].Geophysical Prospecting,2012,60(4):733-747

[14] HE B H,WU G C.Robust full-waveform inversion with incomplete refraction[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:150-153

[15] MASONI I,BROSSIER R.Alternative misfit functions for FWI applied to surface wave[J].Expanded Abstracts of 75thAnnual Internat EAGE Mtg,2013:Th P10 13

[16] NANGOO T P.Seismic full-waveform inversion of 3D field data——from the near surface to the reservoir [D].London:University of Imperial College,2013

[17] FOMEL S.Applications of plane-wave destruction filters[J].Geophysics,2002,67(10):1946-1960

[18] HALE D.Local dip filtering with directional Laplacians[R].CWP report,2007

[19] 王懷野,張科,李言俊.一種自適應(yīng)各向異性高斯濾波方法[J].計算機工程與應(yīng)用,2004,40(10):18-19 WANG H Y,ZHANG K,LI Y J.An adaptive anisotropic Gaussian filter for noise reduction[J].Computer Engineering and Applications,2004,40(10):18-19

(編輯：陳 杰)

Strategy study on full waveform inversion for the land seismic data

WANG Jie,HU Guanghui,LIU Dingjin,SHAO Wenchao,WANG Pengyan

(SinopecGeophysicalResearchInstitute,Nanjing211103,China)

Full waveform inversion of marine seismic data has been successfully applied,however,the application of this technology for onshore seismic data faces many difficulties.In this study,a robust strategy of full waveform inversion for onshore seismic data is presented.Firstly,the phase-matched cross-correlation cost function is adopted to improve the applicability of full waveform inversion.Secondly,dynamic wavenumber domain gradient preprocessing method is used to suppress gradient noise.Finally,an adaptive Gaussian smoothing operator is applied to improve the quality of inversion results.The inversion case of two-dimensional onshore seismic data shows that this strategy presented in this study can effectively realize high wavenumber onshore seismic data velocity modeling.

fullwaveform inversion,cross-correlation cost function,sub-regional wavenumber domain filtering,adaptive anisotropic Gaussian smoothing

2016-11-10;改回日期：2016-12-01。

王杰(1988—),男,碩士,工程師,主要從事全波形反演速度建模研究。

國家科技重大專項(2016ZX05014-001-002)資助。

P631

A

1000-1441(2017)01-0081-08

10.3969/j.issn.1000-1441.2017.01.010

This research is financially supported by National Science and Technology Major Project of China (Grant No.2016ZX05014-001-002).