渝東南南川地區(qū)五峰組—龍馬溪組頁巖氣層地應(yīng)力數(shù)值模擬及有利區(qū)預(yù)測(cè)

劉 明,楊瑞青,楊風(fēng)麗,劉昊娟,張志萍,王 瑋,戶盼盼

1.中國(guó)石化 華東油氣分公司 勘探開發(fā)研究院,南京 210000;2.同濟(jì)大學(xué) 海洋與地球科學(xué)學(xué)院,上海 200092

南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組頁巖氣資源豐富[1],但受中新生代以來多期構(gòu)造應(yīng)力的疊加作用[2-6],導(dǎo)致南川地區(qū)現(xiàn)今地應(yīng)力復(fù)雜、方位變化快,給頁巖氣開發(fā)帶來了很大的難度。因此,采用有效的應(yīng)力場(chǎng)預(yù)測(cè)模擬方法,對(duì)南川地區(qū)頁巖氣應(yīng)力場(chǎng)進(jìn)行預(yù)測(cè)研究,是提高該區(qū)頁巖氣產(chǎn)能與優(yōu)化壓裂方案制定的首要任務(wù)之一。

截至目前,地應(yīng)力預(yù)測(cè)已經(jīng)形成了多種研究方法[7-8],為應(yīng)力場(chǎng)研究提供了較好的支持。在眾多方法中,有限元數(shù)值模擬,無疑是應(yīng)力場(chǎng)研究廣泛且有效的研究手段之一,目前該方法已應(yīng)用于板塊碰撞、大陸形變、地震、頁巖氣、油氣運(yùn)移等不同構(gòu)造尺度上的多個(gè)領(lǐng)域[8-11]。在南川及周緣地區(qū),前人也采用有限元數(shù)值模擬計(jì)算方法,對(duì)頁巖氣分布進(jìn)行了一定的預(yù)測(cè)研究,如丁文龍等[12]對(duì)渝東南地區(qū)下志留統(tǒng)龍馬溪組開展了應(yīng)力場(chǎng)模擬,進(jìn)行了頁巖儲(chǔ)層裂縫的分布預(yù)測(cè);李東東[13]對(duì)渝東南周緣武陵山區(qū)晚燕山期—喜馬拉雅早期不同時(shí)期的構(gòu)造應(yīng)力場(chǎng)進(jìn)行了模擬分析;張斗中等[8]利用ANSYS軟件對(duì)南川地區(qū)龍馬溪組現(xiàn)今應(yīng)力場(chǎng)進(jìn)行了有限元數(shù)值模擬等。

雖然前人的研究取得了一定成效,但我們也注意到,由于普遍采用的有限元模擬軟件,多是以力學(xué)理論為基礎(chǔ)的模擬方法。而SHELLS有限元應(yīng)力場(chǎng)模擬方法[14-16],作為一種基于地質(zhì)力學(xué)、流變學(xué)、地殼均衡等原理來進(jìn)行構(gòu)造應(yīng)力場(chǎng)模擬的有效方法[17],針對(duì)復(fù)雜地質(zhì)背景或深埋地下的區(qū)域,往往具有較好的模擬效果。該方法目前已廣泛應(yīng)用于美國(guó)加州[17]、新西蘭[18]、非洲[19]以及中國(guó)下?lián)P子—南黃海地區(qū)[20-21]、南海[22]和遼東灣盆地[23-24]等地區(qū)的區(qū)域應(yīng)力場(chǎng)模擬研究中,并取得了很好的成果。

由此,本次研究將SHELLS有限元應(yīng)力場(chǎng)模擬方法引入到南川地區(qū)頁巖氣的研究中,通過SHELLS應(yīng)力場(chǎng)數(shù)值模擬,明確五峰組—龍馬溪組應(yīng)力場(chǎng)分布,為南川地區(qū)五峰組—龍馬溪組頁巖氣的進(jìn)一步勘探開發(fā)提供依據(jù)。

1 南川地區(qū)地質(zhì)概況

南川地區(qū)位于華南陸塊內(nèi)部,在大地構(gòu)造位置上,南川地區(qū)既屬于四川盆地盆緣過渡帶,又屬于川東—湘鄂西褶皺帶(圖1)。特殊的構(gòu)造位置,決定了研究區(qū)在地質(zhì)演化過程中構(gòu)造作用的復(fù)雜性。已有研究資料表明,研究區(qū)在現(xiàn)今整體處于太平洋板塊、印度板塊的推擠和西伯利亞、菲律賓板塊約束的地球動(dòng)力學(xué)背景控制下[25-28],其東側(cè)主要受川東—湘鄂西褶皺帶自SE向NW的遠(yuǎn)程傳導(dǎo)擠壓作用[29-32],西側(cè)受龍門山自NW向SE的褶皺—沖斷作用[33-34]、北側(cè)受大巴山構(gòu)造帶由NE向SW的擴(kuò)展作用[32,35];南側(cè)受大婁山由SW向NE的擠壓作用[36]控制。在這種多向、應(yīng)力復(fù)合疊加作用下,不但形成了現(xiàn)今隆坳相間的地質(zhì)面貌(圖2),也控制和影響了深部五峰組—龍馬溪組斷裂與構(gòu)造的發(fā)育。

圖2 南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組斷裂分布

通過對(duì)研究區(qū)3D地震資料的精細(xì)解釋可知(圖2),區(qū)內(nèi)以發(fā)育逆沖斷裂和褶皺為特征,其中大部分逆沖斷裂(如龍濟(jì)橋斷裂、平橋西斷裂、平橋東1號(hào)斷裂、平橋東2號(hào)斷裂等)和褶皺(東勝背斜、平橋背斜)主要呈NE-SW向展布,僅在研究區(qū)東部和西部地區(qū)發(fā)育2條近SN向逆沖斷裂(陽春溝斷裂、大千斷裂),揭示出研究區(qū)總體應(yīng)力方向的復(fù)雜性。

2 SHELLS有限元應(yīng)力場(chǎng)數(shù)值模擬

2.1 原理與方法

SHELLS有限元應(yīng)力場(chǎng)模擬方法,以模擬時(shí)期的斷裂體系、深部結(jié)構(gòu)為約束,遵循非線性流變學(xué)、地殼均衡原理,在地學(xué)、力學(xué)、幾何學(xué)等多學(xué)科交叉限制下驗(yàn)證推演[14-16],可適用于不同的區(qū)域板塊構(gòu)造背景(離散、會(huì)聚、走滑)下的應(yīng)力場(chǎng)、應(yīng)變場(chǎng)及速度場(chǎng)研究[17-24]。SHELLS模型的基本原理滿足以下條件和方程:

(1)構(gòu)造應(yīng)變方程。在小應(yīng)變前提下,材料形變滿足虛功原理:

δUs=δVi+δVf

(1)

式中:δUs為內(nèi)力所做的虛功,δVi為體力所作的虛功。δVf為節(jié)點(diǎn)接觸力所作的虛功。

(2)熱應(yīng)變方程。根據(jù)小應(yīng)變前提下,滿足熱應(yīng)變能守恒:

(2)

(3)流變學(xué)原理。以巖石介質(zhì)的不可壓縮性和垂向熱導(dǎo)率不變?yōu)榧僭O(shè)前提,SHELLS模型遵守流變學(xué)理論[18],在地殼及上地幔巖石中滿足[37]:

ε=Aσnexp(-Q/RT)

(3)

式中:ε為應(yīng)變速率;σ為應(yīng)力差;R為氣體常數(shù);Q為蠕變活化能;A,n為實(shí)驗(yàn)常數(shù)。

根據(jù)以上公式可以進(jìn)而推導(dǎo)出流變應(yīng)力(蠕變強(qiáng)度)公式:

(4)

在以上原理的基礎(chǔ)上,可根據(jù)SHELLS模型計(jì)算并預(yù)測(cè)研究區(qū)的應(yīng)力場(chǎng)性質(zhì)及展布特征,并與實(shí)測(cè)結(jié)果進(jìn)行對(duì)比驗(yàn)證[18]。

2.2 模擬過程及參數(shù)賦值

2.2.1 模擬過程

模擬過程總體分5步。(1)構(gòu)造圖數(shù)字化:對(duì)研究區(qū)五峰組—龍馬溪組構(gòu)造圖進(jìn)行數(shù)字化處理;(2)有限元網(wǎng)格化:運(yùn)用SHELLS軟件將數(shù)字化的構(gòu)造圖進(jìn)行有限元網(wǎng)格化處理;(3)對(duì)有限元節(jié)點(diǎn)等元素進(jìn)行模擬參數(shù)(熱流、地形、巖石物性等)賦值;(4)進(jìn)行SHELLS軟件模型計(jì)算;(5)經(jīng)反復(fù)校驗(yàn)合理后,輸出模擬成果圖件。

2.2.2 模擬參數(shù)賦值

合理正確的模擬參數(shù)賦值是成功的關(guān)鍵。SHELLS應(yīng)力場(chǎng)模擬共需要5類參數(shù)賦值,分別是斷裂、熱流值、地形、巖石物性參數(shù)和邊界條件。本次研究在對(duì)南川地區(qū)所需參數(shù)收集和詳細(xì)分析的基礎(chǔ)上,對(duì)南川地區(qū)五峰組—龍馬溪組模擬參數(shù)進(jìn)行了賦值。

(1)斷裂:以研究區(qū)3D地震資料解釋獲得的五峰組—龍馬溪組構(gòu)造圖為基礎(chǔ)(圖2),進(jìn)行斷裂參數(shù)名稱、性質(zhì)、走向、傾向等的描述和賦值(表1)。

表1 南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組主要斷裂特征

(2)熱流:基于南川及周緣地區(qū)大地?zé)崃髦?以前人研究成果為基礎(chǔ)[38],確認(rèn)出的南川地區(qū)大地?zé)崃髦禐?7 mW/m2。

(3)地形:地形參數(shù)主要根據(jù)現(xiàn)今海拔數(shù)據(jù)來確定和賦值(數(shù)據(jù)源自91衛(wèi)圖APP),南川地區(qū)地形(海拔)由西到東在433～1 690 m之間變化。

(4)巖石物性:SHELLS應(yīng)力場(chǎng)模擬需要對(duì)斷裂賦予相應(yīng)的巖石物性參數(shù)。在使用了南川地區(qū)及華南前人對(duì)巖石圈最大深度、平均巖石密度、熱導(dǎo)率巖石物性參數(shù)研究成果基礎(chǔ)上[39-41],對(duì)其他斷裂摩擦系數(shù)、孔隙水密度、平均巖石密度、熱膨脹系數(shù)等一系列參數(shù)主要采用了普遍認(rèn)可的常量賦值(表2)。

(5)邊界條件

邊界條件的設(shè)置涉及到邊界方向和大小兩方面。邊界方向基于中國(guó)大陸GPS速度場(chǎng)方向[42]和研究區(qū)所在的現(xiàn)今地球動(dòng)力學(xué)背景下受到的川東—湘鄂西、龍門山、大巴山、大婁山等多向應(yīng)力擠壓方向綜合分析后確定。邊界條件的數(shù)值大小賦值則基于中國(guó)大陸GPS速度場(chǎng)[42]和研究區(qū)平衡剖面[43]獲得數(shù)值大小基礎(chǔ)上,經(jīng)過多次的模擬實(shí)驗(yàn)對(duì)比分析后,獲得了最合適的邊界條件形變速率(0.002 36 mm/a)。具體的邊界條件設(shè)置見表3,圖3。

表3 南川地區(qū)SHELLS應(yīng)力場(chǎng)模擬邊界條件設(shè)置

圖3 南川地區(qū)SHELLS應(yīng)力場(chǎng)模擬邊界條件設(shè)置

3 模擬結(jié)果與有利區(qū)預(yù)測(cè)

SHELLS應(yīng)力場(chǎng)模擬可以獲得7類成果圖件,包括最大水平主應(yīng)力方向、塊體應(yīng)變率、斷裂滑移速率、塊體運(yùn)移速率、垂向綜合差異應(yīng)力量、地殼及巖石圈厚度等。本次主要選用對(duì)解決研究區(qū)問題比較明顯的前3類成果圖件進(jìn)行分析。

3.1 模擬結(jié)果

3.1.1 最大水平主應(yīng)力方向

最大水平主應(yīng)力方向代表了地塊的運(yùn)動(dòng)方向、板塊內(nèi)部的縮短方向以及裂谷的走向,可分為伸展、走滑和逆沖3種應(yīng)力體制,分別對(duì)應(yīng)正斷層、走滑斷層和逆斷層[44]。

研究區(qū)現(xiàn)今最大主應(yīng)力方向模擬結(jié)果表明,最大主應(yīng)力以擠壓應(yīng)力為主,總體上存在NW-SE、NE-SW、近EW和近SN四個(gè)主應(yīng)力方向和相應(yīng)所在的區(qū)域(圖4)。其中:① NW-SE向區(qū)域:主要分布于平橋背斜東翼以東地區(qū),最大主應(yīng)力方向與斷層夾角較大,夾角近90°。② NE-SW向區(qū)域:主要分布在平橋背斜西翼、東勝背斜和龍濟(jì)橋斷裂西側(cè),最大主應(yīng)力方向與斷層夾角較小。③ 近EW向區(qū)域:主要分布在龍濟(jì)橋斷裂南部和陽春溝斷裂北側(cè),最大主應(yīng)力方向與斷層夾角較大。④ 近SN向區(qū)域:主要分布在陽春溝斷裂南側(cè)以東,最大主應(yīng)力方向與斷層夾角較小。

圖4 南川地區(qū)最大水平主應(yīng)力方向模擬結(jié)果與分區(qū)

3.1.2 塊體應(yīng)變率

塊體應(yīng)變率代表單位時(shí)間應(yīng)變的變化量,通常有張應(yīng)變、扭應(yīng)變和壓應(yīng)變3種基本情況。在SHELLS應(yīng)力場(chǎng)模擬軟件的成果圖中,通常用張應(yīng)變率、扭應(yīng)變率和壓應(yīng)變率3種不同形式的圖標(biāo)指示不同的應(yīng)變率狀態(tài),圖標(biāo)大小代表應(yīng)變率大小。

南川地區(qū)塊體應(yīng)變率模擬結(jié)果顯示,雖然整體以壓應(yīng)變?yōu)橹?但應(yīng)變率的大小和展布的趨勢(shì)方向存在較大差別,基本上可以分為3個(gè)區(qū):(1)低應(yīng)變率區(qū):應(yīng)變率數(shù)量級(jí)別小于等于-18,性質(zhì)以壓應(yīng)變?yōu)橹?方向以NE-SW為主,主要分布平橋背斜東部;(2)中應(yīng)變率區(qū):應(yīng)變率數(shù)量級(jí)別在-18～-17.6之間,性質(zhì)以壓應(yīng)變和壓扭應(yīng)變?yōu)橹?方向包括NW-SE、NE-SW、EW向和SN向,主要分布在平橋東二號(hào)斷裂以東、東勝背斜、平橋背斜、龍濟(jì)橋斷裂和陽春溝斷裂之間的部分地區(qū)和陽春溝地區(qū)南部區(qū)域;(3)高應(yīng)變率區(qū):低、中應(yīng)變率區(qū)以外的地區(qū),應(yīng)變率數(shù)量級(jí)別普遍在-17.6以上,性質(zhì)為壓應(yīng)變,方向包括NW-SE、NE-SW和SN向(圖5)。

圖5 南川地區(qū)塊體應(yīng)變率模擬結(jié)果與分區(qū)

3.1.3 斷裂滑移速率

斷裂滑移速率是單位時(shí)間內(nèi)的斷裂滑動(dòng)位移量,單位為mm/a。在SHELLS應(yīng)力場(chǎng)模擬軟件成果圖中,以不同顏色代表不同的斷層性質(zhì),以顏色寬度代表滑動(dòng)速率的大小。

研究區(qū)斷裂滑移速率模擬結(jié)果顯示(圖6),總體滑移速率范圍在0～0.001 2 mm/a之間,以逆斷層為主。其中,陽春溝斷裂滑移速率為研究區(qū)最高,達(dá)到0.001 2 mm/a;其他斷裂滑移速率普遍偏低,在0～0.000 27 mm/a之間。這些斷層滑移速率的差異,揭示出雖然外邊界作用力相同,但內(nèi)部斷裂的運(yùn)動(dòng)速率是存在大小差別的。

圖6 南川地區(qū)斷裂滑移速率模擬結(jié)果與分區(qū)

總的來說,南川地區(qū)應(yīng)力場(chǎng)模擬結(jié)果顯示出了該地區(qū)應(yīng)力場(chǎng)復(fù)雜,內(nèi)部變化快的特征。

3.2 模擬結(jié)果與實(shí)測(cè)結(jié)果比較

3.2.1 最大水平主應(yīng)力方向

將模擬結(jié)果與研究區(qū)13口實(shí)測(cè)鉆井主應(yīng)力方向?qū)Ρ冉Y(jié)果顯示,總體上具有85%的基本一致性(圖5)。其中:模擬方向與鉆井實(shí)測(cè)值方向完全一致的井有3口(JY194-3、JY10-10、SY1)(表4),主要分布在平橋背斜和東勝背斜南部區(qū)域;誤差范圍在5°～10°的井有4口(JY8、JY10、SY3、SY5),主要分布在平橋背斜和陽春溝區(qū)域;誤差范圍在11°～18°的井有4口(NY1、JY11、JY201、SY13-3),主要分布在龍濟(jì)橋斷裂以西、東勝背斜北部和平橋背斜東翼區(qū)域;誤差在20°以上的井有2口(SY2、SY9-1),主要分布在龍濟(jì)橋斷裂和袁家溝斷裂之間區(qū)域,分析原因可能與這兩口井正好處于NW與NE向2種應(yīng)力轉(zhuǎn)換交接處,SHELLS軟件分辨率難以達(dá)到所致。

表4 南川地區(qū)SHELLS應(yīng)力場(chǎng)模擬結(jié)果與鉆井實(shí)測(cè)值對(duì)比

3.2.2 塊體應(yīng)變率

現(xiàn)今GPS方法獲得的華南黔渝地區(qū)應(yīng)變率性質(zhì)也以壓應(yīng)變?yōu)橹?且數(shù)量級(jí)別也在-17～-18之間[45-46],這與南川地區(qū)模擬獲得的結(jié)果相一致。不同之處在于,南川地區(qū)模擬獲得的應(yīng)變率趨勢(shì)展布方向,在GPS觀測(cè)結(jié)果中沒有顯示出變化,這可能與GPS觀測(cè)的分辨率有關(guān)。

3.2.3 斷裂滑移速率

盡管目前暫無關(guān)于研究區(qū)及周緣斷裂滑移速率的研究成果報(bào)道,但分析認(rèn)為,它與研究區(qū)(圖2)及區(qū)域斷裂發(fā)育的構(gòu)造背景[31, 43, 47]是相符的。

通過對(duì)比可知,南川地區(qū)模擬結(jié)果的總體特征與實(shí)際相符。同時(shí),以上對(duì)比分析揭示的模擬結(jié)果與實(shí)際結(jié)果的較高吻合性說明了本次SHELLS預(yù)測(cè)結(jié)果的實(shí)效性和可靠性。南川地區(qū)復(fù)雜和快速變化的地應(yīng)力特征,可能是周緣復(fù)雜的應(yīng)力環(huán)境、周緣基底的差異性、巖體自身的蠕變和應(yīng)力松弛等因素共同作用的結(jié)果,最終造成了研究區(qū)內(nèi)部最大主應(yīng)力方向、塊體應(yīng)變率、斷裂滑移速率等的復(fù)雜性和差異性。

3.3 有利區(qū)預(yù)測(cè)

模擬結(jié)果與研究區(qū)已知的鉆探開發(fā)成果綜合對(duì)比分析后發(fā)現(xiàn),應(yīng)力場(chǎng)對(duì)裂縫儲(chǔ)層發(fā)育具有一定影響。最大主應(yīng)力方向影響裂縫的開啟性:最大主應(yīng)力方向與旁側(cè)斷層夾角越小,裂縫的開啟性越好,有利于頁巖氣富集;但若與斷層距離過近,頁巖氣則易順斷層逸散;相關(guān)經(jīng)驗(yàn)表明,最大主應(yīng)力方向與旁側(cè)斷層夾角越小,且與斷層距離大于2 km,對(duì)于頁巖氣裂縫儲(chǔ)層發(fā)育越為有利。而塊體應(yīng)變率和斷裂滑移速率則影響了裂縫的數(shù)量和形態(tài):在高應(yīng)變速率和高斷層滑移速率下,較高的應(yīng)力強(qiáng)度往往會(huì)形成走向單一斷層或裂縫,造成頁巖氣逸散或裂縫儲(chǔ)層的儲(chǔ)集能力差;而在中低應(yīng)變率和低斷層滑移速率下,更易形成數(shù)量較多、形態(tài)復(fù)雜的網(wǎng)狀裂縫,有利于頁巖氣裂縫儲(chǔ)層發(fā)育?？傮w上,頁巖氣裂縫儲(chǔ)層有利發(fā)育區(qū)與構(gòu)造裂縫的開啟性程度、數(shù)量、形態(tài)等息息相關(guān)。

由此,我們對(duì)南川地區(qū)頁巖氣儲(chǔ)層有利發(fā)育區(qū)進(jìn)行了評(píng)價(jià)(表5),預(yù)測(cè)出Ⅰ、Ⅱ兩類裂縫儲(chǔ)層發(fā)育有利區(qū)塊(圖7)。Ⅰ類裂縫儲(chǔ)層發(fā)育有利區(qū)塊15個(gè),分別位于:大千斷裂和青龍鄉(xiāng)斷裂以南(1個(gè))、平橋背斜東翼(3個(gè))、平橋背斜(4個(gè))、東勝背斜(2個(gè))、龍濟(jì)橋斷裂南部(2個(gè))、龍濟(jì)橋斷裂與陽春溝斷裂之間(2個(gè))、陽春溝地區(qū)南部(1個(gè)),代表了區(qū)域裂縫開啟性好、裂縫發(fā)育保存條件好的特點(diǎn);Ⅱ類裂縫儲(chǔ)層發(fā)育區(qū)塊,則代表該區(qū)域具有裂縫開啟性較好、或裂縫密度較大的發(fā)育特征。

表5 南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組頁巖氣有利區(qū)預(yù)測(cè)原則

圖7 南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組應(yīng)力場(chǎng)頁巖氣有利區(qū)預(yù)測(cè)

該預(yù)測(cè)結(jié)果將對(duì)南川地區(qū)五峰組—龍馬溪組頁巖氣勘探及其開發(fā)過程中的產(chǎn)能與優(yōu)化壓裂方案制定具有重要的指示意義。

4 結(jié)論

(1)南川地區(qū)上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組斷裂發(fā)育,斷裂特征明顯。斷層性質(zhì)為逆斷層,主斷層以NE-SW向和近SN向?yàn)橹?次斷層總體在走向上與主斷層平行發(fā)育。

(2)研究區(qū)現(xiàn)今最大主應(yīng)力以擠壓應(yīng)力為主,總體上存在NW-SE、NE-SW、近EW和近SN四個(gè)主應(yīng)力方向;現(xiàn)今應(yīng)變率整體以壓應(yīng)變?yōu)橹?存在低應(yīng)變率(數(shù)量級(jí)別≤-18)、中應(yīng)變率(數(shù)量級(jí)別介于-18～-17.6之間)和高應(yīng)變率(數(shù)量級(jí)別≥-17.6)3個(gè)區(qū)域及其相應(yīng)的NE-SW、NW-SE、SN和EW不同展布方向;現(xiàn)今斷裂滑移速率總體以逆斷層為主,范圍介于0～0.001 2 mm/a之間,最高達(dá)到0.001 2 mm/a。

(3)基于模擬結(jié)果與已知鉆探成果綜合對(duì)比分析基礎(chǔ)上,從應(yīng)力場(chǎng)模擬成果揭示出的裂縫開啟性、裂縫發(fā)育信息,對(duì)研究區(qū)頁巖氣裂縫儲(chǔ)層有利發(fā)育區(qū)進(jìn)行了評(píng)價(jià),預(yù)測(cè)出Ⅰ、Ⅱ兩類裂縫儲(chǔ)層發(fā)育有利區(qū)塊。

利益沖突聲明/Conflict of Interests

所有作者聲明不存在利益沖突。

All authors disclose no relevant conflict of interests.

作者貢獻(xiàn)/Authors’Contributions

劉明、劉昊娟和張志萍提供資料數(shù)據(jù);劉明、楊瑞青、楊風(fēng)麗,戶盼盼和王瑋共同參與論文寫作和修改;所有作者均閱讀并同意最終稿件的提交。

The profile and data were provided by LIU Ming, LIU Haojuan and ZHANG Zhiping. The manuscript was drafted and revised by LIU Ming,YANG Ruiqing, YANG Fengli,HU Panpan and WANG Wei. All the authors have read the last version of paper and consented for submission.