頁巖油氣層地震巖石物理計(jì)算方法研究

劉喜武,劉宇巍,霍志周,劉志遠(yuǎn),張劍鋒,高紅偉

(1.中國石油化工股份有限公司石油勘探開發(fā)研究院,北京100083;2.中國石油化工股份有限公司頁巖油氣勘探開發(fā)重點(diǎn)實(shí)驗(yàn)室,北京100083;3.國家能源頁巖油研發(fā)中心,北京100083;4.中國科學(xué)院地質(zhì)與地球物理研究所,北京100029)

?

頁巖油氣層地震巖石物理計(jì)算方法研究

劉喜武1,2,3,劉宇巍1,2,3,霍志周1,2,3,劉志遠(yuǎn)1,2,3,張劍鋒4,高紅偉4

(1.中國石油化工股份有限公司石油勘探開發(fā)研究院,北京100083;2.中國石油化工股份有限公司頁巖油氣勘探開發(fā)重點(diǎn)實(shí)驗(yàn)室,北京100083;3.國家能源頁巖油研發(fā)中心,北京100083;4.中國科學(xué)院地質(zhì)與地球物理研究所,北京100029)

摘要:研究通過數(shù)值計(jì)算地震波場特征直接獲得地震波長尺度巖石物理參數(shù)及其響應(yīng)特征的計(jì)算巖石物理方法?；趲r心資料,從頁巖沉積的隨機(jī)過程出發(fā),以“層”和“紋層”為基本單元,構(gòu)建厘米—毫米(cm—mm)級頁巖小尺度數(shù)值地質(zhì)模型,并隨機(jī)加入有機(jī)質(zhì)、有機(jī)孔、層間縫、垂直縫和基質(zhì)孔隙。在小尺度網(wǎng)格應(yīng)用宏觀巖石物理等效介質(zhì)模型,充分考慮介質(zhì)的非均勻性,將小尺度地質(zhì)模型轉(zhuǎn)化為地球物理參數(shù)模型。基于小尺度地質(zhì)模型與小網(wǎng)格地球物理參數(shù)模型,進(jìn)行不同角度平面非均勻地震波傳播模擬,提取傳播時(shí)差,直接計(jì)算得到地震波長尺度的巖石物理參數(shù)。以勝利油田羅家地區(qū)頁巖油層為例,改變TOC含量、層與紋層不同厚度和組合、裂隙分布,分別計(jì)算彈性參數(shù),并甄選敏感彈性參數(shù),驗(yàn)證了方法的有效性。

關(guān)鍵詞：小尺度地質(zhì)建模;等效介質(zhì)理論;格子法地震波場模擬;計(jì)算巖石物理

頁巖油氣已成為非常規(guī)油氣勘探開發(fā)的重要領(lǐng)域之一。由于頁巖油氣層的微觀復(fù)雜性、非均勻性和各向異性,納米級尺度的儲集空間和甜點(diǎn)要素的宏觀地球物理響應(yīng)規(guī)律不清楚[1]。如何獲得泥頁巖紋層發(fā)育程度、有機(jī)質(zhì)豐度(TOC)、有機(jī)孔、基質(zhì)孔、微裂隙、裂縫等巖石參數(shù)與地震響應(yīng)特征的關(guān)系,是巖石物理和地球物理面臨的挑戰(zhàn),也是地球物理甜點(diǎn)預(yù)測的基礎(chǔ)。目前,針對頁巖油氣儲層的巖石物理模型和響應(yīng)特征研究[2-6],存在儲層微觀尺度與地震宏觀尺度不匹配的問題[7],如何獲得真實(shí)地震波場尺度和頻率段下的頁巖油氣層巖石物理參數(shù)和地震響應(yīng)規(guī)律,對于實(shí)驗(yàn)巖石物理來說面臨較大挑戰(zhàn),而且基于實(shí)驗(yàn)室物理測量的方法和基于數(shù)字巖心圖像處理的計(jì)算巖石物理方法,均不能直接獲得地震波場的動態(tài)巖石物理彈性參數(shù)[8]。

地震計(jì)算巖石物理方法是一個(gè)比較新的概念,就是試圖利用地震波場數(shù)值模擬計(jì)算方法獲取彈性參數(shù),解決巖石物理實(shí)驗(yàn)測量面臨的樣品不足、尺度效應(yīng)等問題,特別是對各向異性、衰減等特征的測量。目前地震巖石物理計(jì)算方法仍然處于小尺度地震波場模擬探索階段,尚未形成有效技術(shù)[8]。需要指出,巖石物理實(shí)驗(yàn)和計(jì)算巖石物理是相互促進(jìn)的,計(jì)算巖石物理并不能取代巖石物理實(shí)驗(yàn),地震波傳播模擬的一些基本參數(shù)還是需要由巖石物理實(shí)驗(yàn)來準(zhǔn)確測定,例如,在不同溫壓條件下不同巖石組分的彈性模量、裂縫的柔度系數(shù)等。

本文首先闡述小尺度地質(zhì)數(shù)值建模、非均勻網(wǎng)格介質(zhì)等效地球物理參數(shù)建模、非均勻非結(jié)構(gòu)網(wǎng)格地震波場數(shù)值模擬、地震波計(jì)算巖石物理方法,探索建立一套有別于實(shí)驗(yàn)巖石物理、數(shù)字巖石物理的地震波計(jì)算巖石物理方法,并形成一種新的頁巖油氣層地震巖石物理參數(shù)計(jì)算方法;最后,以勝利油田羅家地區(qū)頁巖油儲層為例,應(yīng)用提出的地震波計(jì)算巖石物理方法提取彈性參數(shù),分析響應(yīng)特征和響應(yīng)規(guī)律,甄別敏感彈性參數(shù),為頁巖油氣層巖石物理分析和地震甜點(diǎn)預(yù)測研究提供新的思路與手段。

1頁巖小尺度地質(zhì)與地球物理參數(shù)建模方法

從頁巖沉積的隨機(jī)過程出發(fā),采用二維馬爾科夫鏈模型[9-11],以“層”和“紋層”為基本單元,構(gòu)造地震波長尺度(100m左右)的頁巖模型。首先,通過分析已知巖心樣品,定義不同類型巖相的“層”或“紋層”基本單元;然后,估計(jì)轉(zhuǎn)移概率矩陣,生成非均勻波長尺度頁巖模型。具體步驟為：①根據(jù)典型巖心樣品分析結(jié)果,決定巖相、層結(jié)構(gòu)和礦物組分含量;②確定“層”和“紋層”基本單元及其成分;③綜合考慮巖心樣品統(tǒng)計(jì)分析結(jié)果和成分含量,確定轉(zhuǎn)移概率矩陣;④確定不同深度區(qū)間“層”和“紋層”的比例和分布模式;⑤隨機(jī)加入有機(jī)質(zhì)(包括有機(jī)孔),根據(jù)“層”和“紋層”泥質(zhì)含量大小,確定有機(jī)質(zhì)的發(fā)生概率,獲得非均勻的有機(jī)質(zhì)分布;⑥隨機(jī)加入層間縫、垂直縫;⑦依據(jù)巖心樣品統(tǒng)計(jì)得到的孔隙率,加入孔隙。

在小尺度地質(zhì)建模的基礎(chǔ)上,劃分小尺度網(wǎng)格(0.001m×0.001m),對網(wǎng)格內(nèi)包括的非均勻礦物,考慮孔隙,采用自洽理論(self-consistent approximation,SCA)、Brown-Korring廣義Gassmann理論、Hudson裂縫等效理論等生成等效的物理參數(shù)[12]。業(yè)界通用的巖石物理建模方法是將上述理論直接應(yīng)用于整體模型,計(jì)算宏觀響應(yīng)參數(shù)。本文研究將其應(yīng)用于各個(gè)小尺度網(wǎng)格上進(jìn)行等效,這就充分考慮了實(shí)際頁巖的強(qiáng)非均勻性。具體等效步驟與通用巖石物理建模方法一致：①采用SCA計(jì)算不同成熟度干酪根的物理參數(shù)(考慮有機(jī)孔);②采用SCA計(jì)算含泥、干酪根、石英、方解石、孔隙的干巖石模量;③采用Wood公式計(jì)算孔隙中氣、水、油構(gòu)成的混合流體的密度和體積模量;④采用Brown-Korring廣義Gassmann方程進(jìn)行孔隙混合流體替換;⑤采用Hudson公式引入層間縫,對介質(zhì)加入各向異性特征。

2地震巖石物理計(jì)算方法原理

圖1　格子法波場模擬網(wǎng)格計(jì)算示意圖解

(1)

式中：Mi為虛線輪廓線包圍的質(zhì)量,它是i點(diǎn)周圍三角形質(zhì)量(面積乘密度)和的1/3(在本項(xiàng)目計(jì)算中僅用三角形離散);(bi)l和(ci)l為對應(yīng)三角形的形狀系數(shù);(σx)l,(σz)l和(τxz)l是各個(gè)三角形形心處的應(yīng)力;ut和wt分別為x方向和z方向的速度分量。

根據(jù)縱、橫波不同方向入射的波場傳播時(shí)間差可以計(jì)算彈性參數(shù)。如圖2所示,通過計(jì)算縱波通過尺寸為L的模型的起始時(shí)間與結(jié)束時(shí)間,可計(jì)算波場傳播時(shí)間差Δt。(2)式為垂直入射的縱波速度計(jì)算公式,(3)式為各向異性參數(shù)計(jì)算公式。

(2)

式中：vP為縱波速度;L為模型尺寸;Δt為波前記錄的時(shí)差;v0為給定的均勻介質(zhì)縱波速度;x1和x2分別為波前位置。

圖2　小尺度地質(zhì)模型縱波波場傳播快照

(3)

式中：v(θ)為不同入射角θ的P波速度;ε和δ為各向異性參數(shù)。

計(jì)算得到縱、橫波速度后,相應(yīng)的可以采用(4)式 和(5)式計(jì)算楊氏模量、泊松比等彈性參數(shù)(各向同性介質(zhì))。

(4)

(5)

式中：ρ為密度;vP為縱波速度;vS為橫波速度。

3實(shí)際應(yīng)用效果分析

以勝利油田羅家地區(qū)頁巖油層為例,地震巖石物理計(jì)算方法實(shí)際應(yīng)用流程如圖3所示。頁巖小尺度地質(zhì)建模以羅69井為例確定層、紋層結(jié)構(gòu)以及礦物成分,孔隙添加依據(jù)巖心樣品統(tǒng)計(jì)結(jié)果,層間縫和垂直縫隨機(jī)加入;將宏觀等效巖石物理方法用于小尺度網(wǎng)格(0.001m×0.001m)進(jìn)行地球物理參數(shù)建模;為提高地震波場模擬時(shí)差提取的精度,將地質(zhì)模型進(jìn)行了疊置,并放在一個(gè)均勻背景介質(zhì)中,分別計(jì)算TOC變化和裂隙密度變化的彈性參數(shù),分析敏感彈性參數(shù)。

圖3　地震巖石物理計(jì)算流程

3.1小尺度地質(zhì)與地球物理參數(shù)建模

根據(jù)實(shí)際巖心分析得到的主要巖相組合和礦物成分,認(rèn)為頁巖儲層非均勻的關(guān)鍵因素為：紋層結(jié)構(gòu)、層結(jié)構(gòu)、礦物成分;根據(jù)地質(zhì)錄井和測井評價(jià)結(jié)果,將“層”的尺度定義在“10cm”量級,將“紋層”的尺度定義在“mm”量級,而將橫向非均勻的尺度定義在“10m”量級;礦物組分主要考慮：粘土、方解石、石英、有機(jī)質(zhì)。確定的3種“層”單元為：含泥質(zhì)灰?guī)r層、泥質(zhì)灰?guī)r層、灰質(zhì)泥巖層(0.1m×10.0m),確定兩種“紋層”單元：灰質(zhì)紋層、泥質(zhì)紋層(0.001m×10.000m)。孔隙空間考慮基質(zhì)孔隙、有機(jī)孔和裂縫。

需指出的是,盡管“層”單元的尺度為0.1m×10.0m,但這一單元在地質(zhì)模型中實(shí)際上是由100×10000個(gè)離散點(diǎn)構(gòu)成的,依據(jù)單元的礦物成分,這些離散點(diǎn)將分別對應(yīng)方解石、粘土或石英這幾種基本礦物;同理,“紋層”單元是由1×10000個(gè)離散點(diǎn)構(gòu)成的,這些離散點(diǎn)將分別對應(yīng)方解石、粘土、石英或有機(jī)質(zhì)等幾種基本礦物,即可產(chǎn)生0.001m×0.001m分辨率的小尺度模型。依據(jù)上述策略,利用地層沉積的馬爾科夫鏈過程和礦物成分含量來隨機(jī)生成小尺度、精細(xì)地質(zhì)模型,并在小尺度網(wǎng)格等效計(jì)算地球物理參數(shù)。圖4和圖5分別給出了頁巖油工區(qū)的層狀和紋層狀模型數(shù)值計(jì)算結(jié)果。圖中黃色樣點(diǎn)代表方解石;深藍(lán)色樣點(diǎn)代表粘土;天藍(lán)色樣點(diǎn)代表石英;紅色樣點(diǎn)代表TOC。右側(cè)巖心掃描照片為示意圖,說明小尺度數(shù)值地質(zhì)建模能夠反映頁巖層結(jié)構(gòu)特征。

圖6為小尺度巖石物理等效方法宏觀測試結(jié)果。將給出的巖石物理等效方法用于預(yù)測測井曲線。圖中藍(lán)色曲線為預(yù)測的縱波曲線;紅色曲線為預(yù)測的橫波曲線;右側(cè)為相應(yīng)的預(yù)測誤差曲線。縱波曲線預(yù)測誤差很小,橫波曲線預(yù)測的誤差最大為20m/s,表明采用的介質(zhì)等效方法和步驟是有效的。介質(zhì)等效建模考慮了頁巖的各種要素,能夠滿足小尺度地球物理參數(shù)建模的需要。基于巖石物理等效方法在小尺度網(wǎng)格計(jì)算得到小尺度地球物理參數(shù)。

3.2地震計(jì)算巖石物理參數(shù)提取與響應(yīng)規(guī)律研究

以小尺度地質(zhì)模型為基礎(chǔ),進(jìn)行基于波場模擬的地震彈性參數(shù)計(jì)算與提取分析。為充分反映非均勻頁巖儲層的尺度效應(yīng),提高地震波傳播時(shí)差提取的精度,基于上述參數(shù)生成厚度100m(接近一個(gè)波長)、寬1000m的地質(zhì)模型,生成的頁巖模型重復(fù)鋪設(shè),并且放置到一個(gè)均勻介質(zhì)中,得到更大的模型(圖7)。采用格子法進(jìn)行地震波場模擬,圖8 為傾斜平面波穿過儲層的波場快照(圖中兩條黑線所夾的區(qū)域?yàn)閮?。

圖4　數(shù)值計(jì)算得到的頁巖層狀小尺度地質(zhì)模型(右側(cè)為巖心掃描照片)

圖5　數(shù)值計(jì)算得到的頁巖紋層狀小尺度地質(zhì)模型(右側(cè)為巖心掃描照片)a 生成的紋層狀巖相局部(1); b 生成的紋層狀巖相局部(2)

圖6　某頁巖油井縱、橫波速度等效介質(zhì)模擬結(jié)果

圖7　地震波場模擬采用的模型

在地震波場模擬的基礎(chǔ)上,提取計(jì)算彈性參數(shù),研究地震響應(yīng)特征與規(guī)律。計(jì)算表明：隨TOC的增加速度降低,隨裂縫密度增加速度也降低,與實(shí)驗(yàn)室及其理論巖石物理研究規(guī)律和結(jié)論相一致[3]。在研究縱、橫波速度和各向異性參數(shù)變化規(guī)律的基礎(chǔ)上,綜合考慮儲層各向異性參數(shù)和地震波速度,建議TOC和裂縫密度預(yù)測的敏感參數(shù)為：ε/vP,δ/vP。圖9和圖10分別給出了敏感彈性參數(shù)與TOC以及裂縫密度(主要是水平順層縫)的變化關(guān)系;圖11和圖12分別給出了楊氏模量、泊松比等彈性參數(shù)與TOC和水平順層縫裂隙密度變化的關(guān)系。由圖9到圖12可以看出,TOC和裂縫密度預(yù)測的敏感參數(shù)分別為脆性(負(fù)相關(guān))、ε/vP(正相關(guān))和δ/vP(正相關(guān))。注意到當(dāng)水平順層縫足夠大時(shí),出現(xiàn)負(fù)泊松比現(xiàn)象,解釋為：軸向拉伸、裂縫導(dǎo)致橫向膨脹。

圖8　傾斜平面P波穿過儲層的波場快照a 傾斜入射波場; b 傾斜出射波場

圖9　敏感彈性參數(shù)ε/vP(a)和δ/vP(b)與TOC含量變化的關(guān)系

圖10　敏感彈性參數(shù)ε/vP(a)和δ/vP(b)與裂縫密度變化的關(guān)系

圖11　楊氏模量(a)、泊松比(b)與TOC變化的關(guān)系

圖12　楊氏模量(a)、泊松比(b)與水平裂縫密度變化的關(guān)系

4結(jié)束語

1) 探索建立了一套頁巖油氣層地震巖石物理計(jì)算方法,包括小尺度地質(zhì)與地球物理參數(shù)建模、非均勻地震波場模擬和參數(shù)計(jì)算,為明確頁巖油氣層地球物理響應(yīng)規(guī)律,提供了一種全新的計(jì)算巖石物理手段。

2) 針對勝利油田羅家地區(qū)頁巖油儲層,充分考慮頁巖巖相、層和紋層結(jié)構(gòu)、礦物成分,考慮基質(zhì)孔隙、有機(jī)孔、裂隙特征,應(yīng)用地震計(jì)算巖石物理方法,建立微觀復(fù)雜性與宏觀地震彈性參數(shù)之間的關(guān)系,綜合考慮儲層各向異性參數(shù)和地震波速度,建議TOC和裂縫密度預(yù)測的敏感參數(shù)選為ε/vP和δ/vP;可進(jìn)一步研究單參數(shù)、組合參數(shù)變化下的彈性參數(shù)變化規(guī)律,甄選甜點(diǎn)預(yù)測敏感彈性參數(shù)。

3) 地震巖石物理計(jì)算方法的優(yōu)勢在于克服了實(shí)驗(yàn)條件和巖心樣品的限制,通過數(shù)值模擬獲得地震波長尺度的巖石物理參數(shù),得到物理實(shí)驗(yàn)很難直接獲得的參數(shù),如衰減特征、各向異性特征;不足之處在于過程和參數(shù)復(fù)雜,計(jì)算量大,需要并行計(jì)算,也需要巖石物理實(shí)驗(yàn)提供基本參數(shù)。

參考文獻(xiàn)

[1]劉振武,撒利明,楊曉,等．頁巖氣勘探開發(fā)對地球物理技術(shù)的需求[J]．石油地球物理勘探,2011,46(5)：810-818

LIU Z W,SA L M,YANG X,et al.Needs of geophysical technologies for shale gas exploration[J].Oil Geophysical Prospecting,2011,46(5):810-818

[2]MAVKO G,MUKERJI T,DVORKIN J.The rock physics handbook[M].Cambridge:Cambridge University Press,2009:1-10

[3]ZHU Y P,LIU E R,ALEX M,et al.Understanding geophysical responses of shale-gas plays[J].The Leading Edge,2011,30(3):332-338

[4]GUO Z Q,LI X Y,LIU C.Anisotropy parameters estimate and rock physics analysis for the Barnett Shale[J].Journal of Geophysics and Engineering,2014,11(6):1-11

[5]LI Y,GUO Z Q,LIU C,et al.A rock physics model for the characterization of organic-rich shale from elastic properties[J].Petroleum Science,2015,12(2):264-272

[6]WU X,CHAPMAN M,LI X Y,et al.Anisotropic elastic modeling for organic shales[J].Expanded Abstracts of 74thEAGE Conference & Exhibition,2012:314-318

[7]BAYUK I O,AMMERMAN M,CHESNOKOV E M.Upscaling of elastic properties of anisotropic sedimentary rocks[J].Geophysical Journal International,2008,172(2):842-860

[8]陳颙,黃庭芳,劉恩儒.巖石物理學(xué)[M].合肥：中國科學(xué)技術(shù)大學(xué)出版社,2009:400-475

CHEN R,HUANG T F,LIU E R.Rock physics[M].Hefei:China University of Science and Technology Press,2009:400-475

[9]WANG G C,TIMOTHY R.Methodology of organic-rich shale lithofacies identification and prediction:a case study from Marcellus Shale in the appalachian basin[J].Computers & Geosciences,2012,49:151-163

[10]KALE S V,RAI C S,SONDERGELD C H.Petrophysical characterization of Barnett Shale[C].SPE Unconventional Gas Conference Expanded Abstracts.Pittsburgh,Pennsylvania,USA:SPE,2010:17

[11]MICHELENA R J,GODBEY K S,ANGOLA O.Constraining 3D facies modeling by seismic-derived facies probabilities:example from the tight-gas Jonah field[J].The Leading Edge,2009,28(12):10-19

[12]劉喜武,董寧,劉宇巍.裂縫性孔隙介質(zhì)頻變AVAZ反演方法研究進(jìn)展[J].石油物探,2015,54(2)：210-217

LIU X W,DONG N,LIU Y W.Progress on frequency-dependent AVAZ inversion for characterization of fractured porous media[J].Geophysical Prospecting for Petroleum,2015,54(2):210-217

[13]ZHANG J F,LIU T L.P-SV-wave propagation in heterogeneous media:grid method[J].Geophysical Journal International,1999,136(2):431-438

[14]ZHANG J F,VERSCHUUR D J.Elastic wave propagation in heterogeneous anisotropic media using the lumped finite element method[J].Geophysics,2002,67(2):625-638

[15]ZHANG J F,LIU T L.Elastic wave modelling in 3-D heterogeneous media:3-D grid method[J].Geophysical Journal International,2002,150(3):780-799

[16]GAO H W,ZHANG J F.Parallel 3-D simulation of seismic wave propagation in heterogeneous anisotropic media:a grid method approach[J].Geophysical Journal International,2006,165(3):875-888

[17]ZHANG J F.Elastic wave modeling in fractured media with an explicit approach[J].Geophysics,2005,70(5):T75-T85

[18]ZHANG J F,GAO H W.Elastic wave modelling in 3D fractured media:an explicit approach[J].Geophysical Journal International,2009,177(3):1233-1241

(編輯：陳杰)

Study on seismic rock physics computational method for shale reservoir

LIU Xiwu1,2,3,LIU Yuwei1,2,3,HUO Zhizhou1,2,3,LIU Zhiyuan1,2,3,ZHANG Jianfeng4,GAO Hongwei4

(1.PetroleumExplorationandProductionResearchInstitute,SINOPEC,Beijing100083,China;2.SinopecKeyLaboratoryofShaleOil/GasExplorationandProductionTechnology，Beijing100083,China;3.NationalEnergyResearch&DevelopmentCenterofShaleOil，Beijing100083，China;4.InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China)

Abstract:A new scheme for shale rock physics study is presented.Firstly,Small-scale shale geological model is computed by considering layer,fine-layer,TOC,crack and porosity.Then,Geophysical properties are computed on non-uniform small grids by effective medium theory with heterogeneity being considered.Finally,seismic modeling is utilized to compute seismic elastic parameters by traveltime difference based on the small-scale geological model and the small grid geophysical parameters model.The methods are applied to shale oil area in Shengli Oilfield and seismic responses are variation of TOC and crack density.The results demonstrate the effectiveness of the proposed method,which have the similar trends with lab data and theoretical data by cross plot for variation of TOC and crack density with elastic parameters,especially with anisotropic parameters.

Keywords:small-scale geological modeling,effective medium theory,grid method seismic modeling,computational rock physics

文章編號：1000-1441(2016)01-0010-08

DOI:10.3969/j.issn.1000-1441.2016.01.002

中圖分類號：P631

文獻(xiàn)標(biāo)識碼：A

基金項(xiàng)目：國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2014CB239104和2012CB214806)資助。

作者簡介：劉喜武(1970—),男,副教授,主要從事頁巖油氣地球物理與地震各向異性裂縫預(yù)測研究工作。

收稿日期：2015-07-24;改回日期：2015-08-05。

劉喜武,劉宇巍,霍志周,等.頁巖油氣層地震巖石物理計(jì)算方法研究[J].石油物探,2016,55(1):-17

LIU Xiwu,LIU Yuwei,HUO Zhizhou,et al.Study on seismic rock physics computational method for shale reservoir[J].Geophysical Prospecting for Petroleum,2016,55(1):-17

This research is financially supported by the National Key Basic Research Program of China (973 Program) (Grant Nos.2014CB239104,2012CB214806).