郭昱琴 龔麗文 郭瑛霞 胡久常 李 盛
1 海南省地震局,海口市美苑路49號(hào),570203 2 重慶市地震局,重慶市紅黃路339號(hào),401147 3 青海省地震局,西寧市興海路1號(hào),810001
早在20世紀(jì)80年代,地震相關(guān)分析和地質(zhì)建模研究中就引入了數(shù)值模擬方法[1],特別是有限元數(shù)值模擬方法,可用于三維非線性多場(chǎng)耦合的復(fù)雜建模分析,目前已成為主流建模分析方法之一[2-7]。本文利用有限元分析方法,基于地質(zhì)構(gòu)造資料構(gòu)建五指山臺(tái)臺(tái)址幾何模型,結(jié)合GPS和重力觀測(cè)結(jié)果,分析五指山臺(tái)定點(diǎn)形變觀測(cè)數(shù)據(jù)趨勢(shì)變化的物理意義,為該區(qū)域地殼應(yīng)力與應(yīng)變狀態(tài)的長(zhǎng)期變化研究提供資料。
五指山臺(tái)位于五指山市沖山鎮(zhèn)阿里山度假村旁太平山半山腰(109.53°E,18.79°N),臺(tái)站正南方25 km處發(fā)育有尖峰-吊羅深大斷裂,距臺(tái)站10 km處則有一條隱伏平移斷層,該隱伏斷裂極有可能是尖峰-吊羅深大斷裂的延長(zhǎng)發(fā)育斷裂,對(duì)臺(tái)站觀測(cè)有一定影響。研究區(qū)域廣泛出露二疊紀(jì)和三疊系花崗巖[8]。臺(tái)站洞體基巖堅(jiān)硬完整、致密均勻,測(cè)點(diǎn)周邊觀測(cè)環(huán)境較為穩(wěn)定,為臺(tái)站提供了較好的觀測(cè)背景場(chǎng)。臺(tái)站為無人值守臺(tái)站,“十五”期間架設(shè),觀測(cè)山洞由部隊(duì)防空洞改建。山洞內(nèi)觀測(cè)儀器有DSQ型水管傾斜儀、VP型寬頻帶垂直擺傾斜儀和SS-Y型銦瓦棒伸縮儀,其中DSQ型水管傾斜儀和SS-Y型銦瓦棒伸縮儀布設(shè)在主洞室內(nèi),兩者相同分量的儀器部件平行布設(shè)且共用同一槽體,分別監(jiān)測(cè)該區(qū)域的應(yīng)變和傾斜變化;VP型寬頻帶傾斜儀布設(shè)在主洞室南側(cè)的側(cè)洞室。
五指山地區(qū)在大地構(gòu)造上屬于西太平洋地殼構(gòu)造不同發(fā)展階段的大陸邊緣區(qū),成土母巖為花崗巖和流紋巖,地貌主要為陡坡火山地貌和侵入巖體構(gòu)造地貌。區(qū)內(nèi)出露地層有長(zhǎng)城系抱板群峨文嶺組,是海南島結(jié)晶基底的主要組成部分,侵入巖主要形成于二疊紀(jì)、三疊紀(jì)、侏羅紀(jì)和白堊紀(jì)。研究區(qū)位于昌江-瓊海構(gòu)造帶以南、尖峰-吊羅斷裂帶以北,主要發(fā)育NW向構(gòu)造[9]。
以橫貫海南島的EW向王五-文教大斷裂及九所-陵水大斷裂為界,2條斷裂之間為五指山地層分區(qū)。區(qū)內(nèi)發(fā)育前震旦系、震旦系、古生界地層,前震旦系地層以具復(fù)理石韻律結(jié)構(gòu)的砂泥質(zhì)巖及火山巖為主,少量碳酸鹽巖;震旦系、下古生界地層以具復(fù)理石韻律結(jié)構(gòu)的粉砂泥質(zhì)巖石為主,少量砂巖、碳酸鹽巖及酸性、基性、超基性火山巖,具地槽沉積特點(diǎn)。石炭紀(jì)-早三疊世為相對(duì)穩(wěn)定的準(zhǔn)地臺(tái)發(fā)展階段[10]。
觀測(cè)洞室的臺(tái)基為角閃黑云二長(zhǎng)花崗巖,厚度較大,約為40 m,下部為薄層巖層碳質(zhì)粉砂巖,有時(shí)與呈條帶狀、條帶紋白云巖互層,其下與碳質(zhì)板巖呈過渡狀態(tài),或夾碳質(zhì)板巖,厚約50 m,其下與石碌群上亞群為假整合接觸。
由于不同巖性的物質(zhì)成分和結(jié)構(gòu)不相同,所以其物理特性和巖石力學(xué)參數(shù)具有很大的差異。海南島巖漿巖分布廣泛,具有多期次活動(dòng)特征,其中侵入巖以海西-印支期斑狀花崗巖、燕山期花崗巖及花崗閃長(zhǎng)巖為主[11]。海南島內(nèi)花崗巖出露非常廣泛,樂東-五指山-萬寧地區(qū)出露面積約800 km2中二疊世同碰撞強(qiáng)變形花崗巖。結(jié)合五指山區(qū)域的地層分布特征,該區(qū)域出露的巖性主要為二疊紀(jì)花崗巖,以黑云母二長(zhǎng)花崗巖和角閃黑云二長(zhǎng)花崗巖為主。
長(zhǎng)期以來,由于植被覆蓋嚴(yán)重,對(duì)海南島的地質(zhì)研究較少。統(tǒng)計(jì)雷瓊地區(qū)MS≥4.5地震和局域小震的震源機(jī)制解可知,瓊北至雷州半島地區(qū)平均主壓應(yīng)力方向?yàn)镹NW向,瓊南地區(qū)平均主壓應(yīng)力方向?yàn)榻麼S向;斷層面走向主要為NW向和NE向,斷層面傾角大。全區(qū)主壓應(yīng)力軸(與海岸線正交)具有明顯的扇形分布特征,反映了東南沿海構(gòu)造應(yīng)力場(chǎng)的基本狀況,即來源于太平洋板塊向西的擠壓和印度洋板塊向北和東的推擠[12]。
從單個(gè)GPS站點(diǎn)數(shù)據(jù)(圖1)來看,海南三亞(HISY)站2015-11以來呈SW向偏移,海南瓊中(QION)站2010-06以來呈NW向偏移,海南???HIHK)站則一直較為平穩(wěn)。從GPS站點(diǎn)間基線數(shù)據(jù)(圖1)來看,HIHK-QION雖表現(xiàn)為緩慢縮短擠壓狀態(tài),但HIHK-HISY、QION-HISY均表現(xiàn)為快速拉伸狀態(tài),標(biāo)志著整個(gè)海南島主要處于地殼拉伸狀態(tài)。地殼拉伸可能會(huì)引起地殼下降(沉),而不均勻下降可引起地殼S傾,從而使山體出現(xiàn)S傾。
從海南重力場(chǎng)變化(圖2)來看,2015-09~2016-09海南呈現(xiàn)東北部正變化、南部和西部負(fù)變化的過程,零值線穿瓊東南而過,形成重力場(chǎng)正負(fù)異常變化高梯度區(qū),五指山臺(tái)周邊地區(qū)重力變化約為-10 μGal。2016-09~2017-10海南重力場(chǎng)變化比較平穩(wěn),基本在30 μGal以內(nèi),南部為正變化,北部為負(fù)變化,五指山臺(tái)周邊地區(qū)重力變化約為10 μGal。故2015-09以來,五指山臺(tái)周邊地區(qū)的重力場(chǎng)變化經(jīng)歷了先負(fù)后正的變化過程。
圖2 海南重力場(chǎng)動(dòng)態(tài)變化Fig.2 Dynamic variation of gravity field in Hainan
研究結(jié)果表明,重力值減小說明地表隆升或地下物質(zhì)減少或兩者疊加;重力值增大說明地表沉降或地下物質(zhì)增多或兩者疊加。比較瓊中臺(tái)重力非潮汐變化(圖3)和GPS站基線變化(圖1)可知,2017-01以來瓊中臺(tái)重力非潮汐變化處于小幅上升變化趨勢(shì),瓊中GPS基線垂直分量自2017年中期開始處于緩慢下降過程,兩者的變化均表明2017年初以來瓊中臺(tái)地區(qū)地表小幅沉降,而2017-01以來GPS基線瓊中-三亞呈伸長(zhǎng)趨勢(shì)(圖2)。結(jié)合這些變化特征認(rèn)為,2017年以來海南南部地區(qū)均處于地表沉降的過程,且北部的沉降幅度小于南部。
圖3 瓊中臺(tái)重力非潮汐變化和五指山臺(tái)水管傾斜儀NS分量變化Fig.3 Non tidal variation of gravity at Qiongzhong station and NS component variation of water pipe tiltmeter at Wuzhishan station
選用ANSYS Workblench 16.0軟件進(jìn)行建模分析。
自2015年觀測(cè)以來,五指山臺(tái)3套儀器工作正常,積累了豐富的觀測(cè)資料。經(jīng)對(duì)比分析發(fā)現(xiàn):1)在基線布設(shè)完全一致的情況下,傾斜類儀器的映震能力遠(yuǎn)高于應(yīng)變類儀器;2)2套傾斜儀能記錄到明顯的固體潮,而應(yīng)變儀記錄的固體潮幅度較小,特別是EW分量幾乎記錄不到同震響應(yīng);3)2套傾斜儀的長(zhǎng)期趨勢(shì)不一致,垂直擺趨于N傾,而水管儀出現(xiàn)S傾現(xiàn)象。
為便于建模分析,進(jìn)行以下假設(shè):1)用簡(jiǎn)單的縱向拉伸地質(zhì)剖面圖模擬實(shí)際的山體模型;2)由于模型是針對(duì)山體整體的特性進(jìn)行分析,故對(duì)巖性和構(gòu)造要素作簡(jiǎn)單的歸一化處理;3)由于大部分研究結(jié)果都是關(guān)于應(yīng)力方向和GPS速率的,為便于顯示,假設(shè)應(yīng)力大小與楊氏模量量級(jí)相同。
4.2.1 幾何模型建立
為便于模型分析、提高運(yùn)算效率,選用應(yīng)力45°方向AB剖面(圖4)進(jìn)行建模,并結(jié)合實(shí)際地質(zhì)資料,參考巖石力學(xué)參數(shù)表設(shè)定模型的材料屬性(表1)。
表1 五指山臺(tái)地殼介質(zhì)巖石力學(xué)參數(shù)匯總Tab.1 Summary of rock mechanical parameters of crustal media at Wuzhishan station
圖4 分析模型Fig.4 Analysis model
4.2.2 模型分析
1)網(wǎng)格劃分。選用Workblench中的自動(dòng)網(wǎng)格劃分工具,為保證模型網(wǎng)格密度和計(jì)算質(zhì)量,綜合考慮計(jì)算精度和計(jì)算量,選擇網(wǎng)格尺寸為1.5 m,其他參數(shù)均采用默認(rèn)值,滿足模型需求共生成95 400個(gè)單元、409 128個(gè)節(jié)點(diǎn)。
2)接觸關(guān)系。Workblench提供多種接觸類型,但由于不同巖層間的摩擦系數(shù)和粘合系數(shù)很難確定,本文采用簡(jiǎn)單的綁定接觸以達(dá)到模型需求。
3)邊界條件。載荷和約束是ANSYS軟件求解的邊界條件,它們是以所選單元的自由度形式定義的。載荷主要選擇力載荷工具的壓力載荷,使用固定底部、兩邊拉張及兩邊壓縮模型計(jì)算的位移和應(yīng)力結(jié)果一致。本文在此只根據(jù)區(qū)域構(gòu)造應(yīng)力方向計(jì)算總體變形結(jié)果和等效應(yīng)力分布結(jié)果,大小設(shè)定為1 010 Pa。
本模型主要研究AB剖面?zhèn)让?5°受壓時(shí)臺(tái)址山體的變形特征和等效應(yīng)力分布特征。但由于區(qū)域應(yīng)力大小的不確定性,本文同時(shí)模擬了兩邊壓縮和兩邊拉張工況(圖5和表2)??梢园l(fā)現(xiàn),兩邊壓縮和兩邊拉張工況下,除了變形大小和應(yīng)力影響范圍變大外,其他結(jié)果與側(cè)面45°受壓時(shí)趨于一致。等效應(yīng)力分布結(jié)果顯示,最大應(yīng)力為78.288 MPa,兩邊受壓或者兩邊受拉最大總體變形結(jié)果一致,受壓應(yīng)力45°方向后,在剖面A巖層上部應(yīng)力的方向發(fā)生明顯變化。
表2 對(duì)剖面等效應(yīng)力的模型結(jié)果Tab.2 Model results of equivalent stress of profile
圖5 兩邊壓縮、兩邊拉張、側(cè)面45°受壓情況下有限元模型分析結(jié)果Fig.5 Finite element model analysis results under compression on both sides, tension on both sides and 45° compression on one side
總體變形結(jié)果顯示,最大總體變形為19.116 mm,水平最大變形為17.881 mm,垂向最大變形為6.803 mm。黑云母二長(zhǎng)花崗巖和斷層碎裂巖更容易發(fā)生變形,而角閃黑云二長(zhǎng)花崗巖變形幅度較小,這可能與研究區(qū)的巖石占比有關(guān)。受巖層之間摩擦力和其他接觸力的影響,模型上部的變形幅度明顯大于下部,當(dāng)下部持續(xù)縮短,便容易形成上拱的背斜,這與實(shí)際的地貌相符合。由于巖層的傾向趨于SE,造成變形趨勢(shì)也趨于SE。
結(jié)合實(shí)際觀測(cè)資料(圖6)發(fā)現(xiàn),2015-01-01~2019-12-01水管傾斜儀數(shù)據(jù)曲線整體呈SE傾(NS分量趨勢(shì)下降,EW分量趨勢(shì)上升),反映觀測(cè)巖層呈SE傾,與數(shù)值模擬結(jié)果相似;垂直擺傾斜儀數(shù)據(jù)曲線長(zhǎng)期趨勢(shì)不明顯(自觀測(cè)以來,NS分量趨勢(shì)上升,EW分量數(shù)據(jù)趨勢(shì)無規(guī)律),部分時(shí)間段出現(xiàn)較大的波動(dòng),之后趨于恢復(fù),其物理意義需進(jìn)一步討論;伸縮儀NS分量呈長(zhǎng)時(shí)間張性(上升),EW分量2018-08-06急速壓性變化(下降),2018-11-26趨勢(shì)轉(zhuǎn)折緩慢張性(上升),整體呈張性趨勢(shì)。
圖6 五指山臺(tái)形變儀器觀測(cè)數(shù)據(jù)曲線Fig.6 Observation data curve of deformation instrument at Wuzhishan station
在建臺(tái)初期,水管傾斜儀與伸縮儀位于同一儀器墩、同一水平面上,由于構(gòu)造地質(zhì)作用巖層呈小角度S傾。當(dāng)受到構(gòu)造應(yīng)力擠壓及南部地表沉降作用后,地層發(fā)生塑性縮短變形并擠壓上層巖石,造成上覆巖層上拱形成一定弧度,使對(duì)應(yīng)基線拉長(zhǎng),這是伸縮儀兩分量觀測(cè)結(jié)果為張性變化的原因。而水管傾斜儀利用連通器原理觀測(cè)基線兩端水位變化,由于基線較長(zhǎng),兩側(cè)的高差容易發(fā)生變化,因此2017年起其NS分量呈緩速S傾長(zhǎng)期變化趨勢(shì)。垂直擺傾斜儀通過擺體的傾角記錄基點(diǎn)巖層的傾斜變化,NS分量長(zhǎng)期趨勢(shì)為N傾,與儀器受零漂影響有關(guān);EW分量趨勢(shì)不規(guī)律,部分時(shí)段出現(xiàn)較大波動(dòng),之后趨于恢復(fù),其物理意義需進(jìn)一步討論。
1)五指山地區(qū)構(gòu)造較復(fù)雜,因地質(zhì)資料的稀缺,僅以出露地質(zhì)圖作為山體模型,基底巖性差異未知,且長(zhǎng)期受NNW向擠壓應(yīng)力影響,在一定程度上影響著區(qū)域應(yīng)力場(chǎng)的大小和方向。
2)有限元模型分析結(jié)果顯示,在NNW向壓力作用下,模型上部的變形幅度明顯大于下部,下部持續(xù)縮短便容易形成上拱,這與GPS及重力場(chǎng)結(jié)果相符。由于巖層的傾向趨于SE,造成變形的趨勢(shì)也趨于SE傾。
3)結(jié)合五指山臺(tái)各儀器觀測(cè)數(shù)據(jù)進(jìn)行分析可知,水管傾斜儀利用連通器原理觀測(cè)基線兩端水位變化,由于基線較長(zhǎng),基線兩側(cè)的高差容易發(fā)生變化,這是2017年起水管儀NS分量出現(xiàn)緩速S傾長(zhǎng)期變化趨勢(shì)的原因。數(shù)據(jù)曲線整體呈SE傾(NS分量趨勢(shì)下降,EW分量趨勢(shì)上升),反映觀測(cè)巖層呈SE傾,這與數(shù)值模擬結(jié)果相似。垂直擺通過擺體的傾角記錄基點(diǎn)巖層的傾斜變化,其長(zhǎng)期趨勢(shì)不明顯(自觀測(cè)以來NS分量趨勢(shì)上升,EW分量無規(guī)律),部分時(shí)間段出現(xiàn)較大波動(dòng),之后趨于恢復(fù),其物理意義需進(jìn)一步討論。伸縮儀NS分量呈長(zhǎng)時(shí)間張性(上升),EW分量2018-08-06急速壓性變化(下降),2018-11-26趨勢(shì)轉(zhuǎn)折呈緩慢張性(上升)趨勢(shì),整體呈張性。當(dāng)受構(gòu)造應(yīng)力擠壓及南部地表沉降作用后,地殼發(fā)生塑性縮短變形并擠壓上層巖石,使上覆巖層上拱形成一定弧度,對(duì)應(yīng)基線拉長(zhǎng),很好地解釋了伸縮儀兩分量觀測(cè)結(jié)果整體為張性變化的原因。