基于數(shù)字化水位的瓊東北地區(qū)構(gòu)造應(yīng)力場(chǎng)時(shí)序特征分析＊

解曉靜 張藝瀟 張 帆 李 盛 鄭在壯

1) 海南省地震局，海南?？?570203

2) 中國(guó)地震局地球物理研究所，北京 100081

引言

區(qū)域地殼構(gòu)造應(yīng)力場(chǎng)的分析是一項(xiàng)較困難復(fù)雜的事情，通常通過(guò)測(cè)量應(yīng)變實(shí)現(xiàn)應(yīng)力的測(cè)量，其耗資較高。研究表明，承壓含水層井水位的變化反映了含水層孔隙壓力的變化，而含水層孔隙壓力的變化與含水層所受的應(yīng)力狀態(tài)有密切的關(guān)系[1]。測(cè)量地殼深部彈性應(yīng)力應(yīng)變的方法有觀測(cè)記錄斷層活動(dòng)區(qū)附近鉆井內(nèi)流體的壓力變化，而測(cè)量鉆井水位的變化是研究流體壓力變化的方法之一[2]。流動(dòng)性和難壓縮性是地殼結(jié)構(gòu)中活躍組分之一的地下水所存在的主要性狀，而當(dāng)?shù)叵滤诤畬又行纬梢粋€(gè)封閉條件較好的承壓系統(tǒng)時(shí)，其水位的升降變化就能較客觀地反映出周邊區(qū)域地殼應(yīng)力應(yīng)變的狀態(tài)。構(gòu)造應(yīng)力變化將導(dǎo)致孔隙體積的變化，從而造成孔隙壓力的改變；井孔水位的變化量與巖石中的孔隙壓力有關(guān)，進(jìn)而與構(gòu)造應(yīng)力有關(guān)。因此，利用已有承壓井水位連續(xù)觀測(cè)資料，研究區(qū)域現(xiàn)今構(gòu)造應(yīng)力場(chǎng)狀態(tài)將是一條既有效又經(jīng)濟(jì)的途徑[2]。

數(shù)字化觀測(cè)以來(lái)，積累了較豐富的井水位觀測(cè)資料，而利用承壓井水位觀測(cè)資料進(jìn)行區(qū)域地殼應(yīng)力應(yīng)變狀態(tài)的反演，從而進(jìn)行地震形勢(shì)的分析研究等，前人已做了較多定性的分析工作[3]：李永善[4]研究了構(gòu)造應(yīng)力引起地下水位變化的特征，并推導(dǎo)出了含水層水位變化與其垂直向應(yīng)力變化之間的關(guān)系式；張昭棟等[5-6]利用含水層參數(shù)、固體潮效應(yīng)和氣壓效應(yīng)等3 種方法反演含水層應(yīng)力應(yīng)變，并進(jìn)行了分析比較；黃輔瓊等[2]利用華北地區(qū)40 多口深井水位動(dòng)態(tài)變化資料，定性地分析了大華北地區(qū)的構(gòu)造應(yīng)力場(chǎng)狀態(tài)；孫小龍等[7]利用華北地區(qū)63 口井的水位變化和井含水層水文地質(zhì)參數(shù)資料，反演出華北地區(qū)構(gòu)造應(yīng)力場(chǎng)變化圖像；而丁風(fēng)和等[8]和其他研究人員[9-13]利用氣壓系數(shù)和維尼迪科夫潮汐調(diào)和分析結(jié)果（O1、M2波潮汐因子）等參數(shù)，反演得到井孔含水層部分介質(zhì)參數(shù)，進(jìn)而在水平層狀含水層模式下，定量地計(jì)算出觀測(cè)井的垂直向應(yīng)力場(chǎng)的動(dòng)態(tài)變化過(guò)程。

鑒于前人的研究，既有定性的分析也有定量的計(jì)算。本文擬借鑒丁風(fēng)和等[8]定量計(jì)算的方法，選擇資料可靠性較好的數(shù)字化井水位資料，利用氣壓系數(shù)和O1、M2波潮汐因子，反演井孔含水層部分介質(zhì)參數(shù)，定量分析海南省區(qū)域構(gòu)造應(yīng)力場(chǎng)的時(shí)序變化過(guò)程，從而分析探尋與地震有關(guān)的部分地球物理信息。

1 資料的分析與處理

1.1 觀測(cè)井基本概況

本文選取了位于瓊東北地區(qū)的5 口數(shù)字水位觀測(cè)井，其中，?？赯K26 井、向榮村井、火山流體井，這3 口井均沿馬裊—鋪前斷裂帶分布，且為活斷裂，井距離海岸較近，井水位效應(yīng)表現(xiàn)為海潮潮汐顯著、氣壓潮汐次之、固體潮汐最弱的復(fù)合潮汐形態(tài)；文昌潭牛井和瓊海加積井，主要沿文昌—瓊?！齺啍嗔褞Х植迹@2 口井水位潮汐效應(yīng)表現(xiàn)為以固體潮汐顯著、氣壓潮汐次之、海潮潮汐最弱的復(fù)合潮汐形態(tài)（圖1）。各井承壓性均較好，觀測(cè)資料相對(duì)較長(zhǎng)且連續(xù)穩(wěn)定。表1 為各井經(jīng)緯度、井深、含水層巖性、地下水類型、水位埋深等基本情況。

表1 開展研究的瓊東北地區(qū)5 口井的基本情況Table 1 Basic situation of 5 wells in northeast Hainan

圖1 研究的5 口觀測(cè)井空間位置分布圖Fig.1 Spatial distribution of 5 observation wells studied

1.2 數(shù)據(jù)的整理與處理

本研究從海南省地球物理臺(tái)網(wǎng)中心數(shù)據(jù)庫(kù)下載井孔數(shù)字化水位和氣壓的整點(diǎn)值數(shù)據(jù)，選取研究資料為2008—2020 年，其中，文昌潭牛井因數(shù)據(jù)庫(kù)故障原因，2008—2012 年氣壓資料缺失，故該井研究起始時(shí)間為2013 年。首先將水位、氣壓兩測(cè)項(xiàng)數(shù)據(jù)進(jìn)行逐月檢查，其存在的突跳或錯(cuò)誤數(shù)據(jù)需做刪除處理，并按999999 進(jìn)行缺數(shù)標(biāo)記。水位數(shù)據(jù)整理成以米為單位的數(shù)據(jù)文件。

本文研究的5 口水位井地處海南，屬于我國(guó)沿海地區(qū)，其水位除了受固體潮、氣壓潮的影響，同時(shí)還受海潮的影響，表現(xiàn)出復(fù)合潮汐效應(yīng)。井水位受海潮荷載影響的程度，與井孔距海岸的距離大小的關(guān)系十分明顯，一般距離在5 km 以內(nèi)時(shí)，海潮的影響顯著；超過(guò)5 km 后，海潮的影響明顯減弱。井水位可以分離成363 個(gè)潮汐分波，主要有Q1、O1、M1、S1K1、J1、OO1、2N2、N2、M2、L2、S2K2、M3。其中，O1、S1K1分波主要受海潮影響，S2K2分波主要受溫度和氣壓影響，M2分波主要受固體潮影響。瓊北地區(qū)的海口ZK26 井、向榮村井、火山流體井，這3 口井距離海岸均約5 km，再根據(jù)各井水位潮汐分波振幅特征的分析（圖2），可看出該3 口井振幅較大的分波主要是O1、S1K1波，其次是S2K2波，因此，主要受海潮潮汐影響，進(jìn)而選擇其中的主要分波O1波進(jìn)行維尼迪科夫潮汐調(diào)和分析得出3 口井的潮汐因子。而瓊海加積井和文昌潭牛井，距離海岸約為20 km，由圖3 可看出，該2 口井水位的潮汐分波主要是M2波，受固體潮汐影響顯著，因此，選擇M2分波進(jìn)行維尼迪科夫潮汐調(diào)和分析計(jì)算2 口井的潮汐因子。圖3 為各井潮汐因子月值曲線。

圖2 瓊東北地區(qū) 5 口井水位潮汐分波振幅特征圖Fig.2 Characteristics of tidal wave amplitude of water level of 5 wells in northeast Hainan

圖3 瓊東北地區(qū) 5 口井水位潮汐因子月值曲線圖Fig.3 Monthly value curve of water level tidal factor of 5 wells in northeast Hainan

各井氣壓系數(shù)的主要獲取步驟為：首先對(duì)水位和氣壓數(shù)據(jù)進(jìn)行別爾采夫?yàn)V波，其次對(duì)濾波后的數(shù)據(jù)進(jìn)行一階差分處理，最后根據(jù)差分結(jié)果，構(gòu)建一元回歸模型來(lái)擬合氣壓系數(shù)[8]。其大小因井而異，變化范圍一般為1—10 mm/hPa，個(gè)別井特殊時(shí)段下氣壓系數(shù)值會(huì)出現(xiàn)負(fù)值或較大值，本文獲取的各井氣壓系數(shù)值基本處于此合理變化范圍之內(nèi)。圖4 為各井氣壓系數(shù)月值曲線。

圖4 瓊東北地區(qū)5 口井水位氣壓系數(shù)月值曲線圖Fig.4 Monthly value curve of water level pressure coefficient of 5 wells in northeast Hainan

因多數(shù)井水位觀測(cè)資料的趨勢(shì)性下降變化，直接受區(qū)域地下水的開采等影響，且降水量的年周期變化也對(duì)井水位觀測(cè)資料的趨勢(shì)性變化具有干擾影響，為此需采用去趨勢(shì)的分析方法，剔除地下水開采和降水補(bǔ)給等長(zhǎng)周期干擾對(duì)水位的影響[8]。但長(zhǎng)周期干擾（趨勢(shì)性變化）剔除后，降雨和氣壓的年周期變化對(duì)大部分水位的影響依然存在。接著對(duì)剔除趨勢(shì)性變化后的水位進(jìn)行一般矩平去周期分析，進(jìn)一步剔除降水和氣壓對(duì)其的影響，然后再進(jìn)行井水位變化量的計(jì)算[8]。地下水開采及降水量的補(bǔ)給，是絕大多數(shù)地區(qū)多數(shù)井孔水位的主要影響因素，因此，本文采用多次剔除這些干擾因素的方法，再進(jìn)行水位變化量的計(jì)算，結(jié)果才更可靠合理。圖5 為各井水位變化量月值曲線。

圖5 瓊東北地區(qū)5 口井水位變化量月值曲線圖Fig.5 Monthly value curve of water level change of 5 wells in northeast Hainan

經(jīng)過(guò)計(jì)算得出的潮汐因子、氣壓系數(shù)和含水層參數(shù)等，均為月值數(shù)據(jù)，基本能滿足地震的中短期預(yù)測(cè)對(duì)資料間隔的需求。

2 基本理論及參數(shù)求取過(guò)程

2.1 含水層參數(shù)的反演計(jì)算

依據(jù)丁風(fēng)和等[8-9]研究結(jié)果，在理想的承壓含水層模式下，井水位的氣壓系數(shù)和潮汐因子可分別表示為：

上兩式聯(lián)合可得到：

式中，n為含水層的孔隙度，α為固體骨架的體積壓縮系數(shù)，β為含水層內(nèi)水的體積壓縮系數(shù)。ρg為水的重度，在研究過(guò)程中取ρg=9.8×103hPa/m。

潮汐因子Bg通過(guò)維尼迪科夫（Venedikov）潮汐調(diào)和分析獲??；氣壓系數(shù)Bp主要利用水位數(shù)據(jù)和同井觀測(cè)的氣壓數(shù)據(jù)經(jīng)過(guò)別爾采夫?yàn)V波和一階差分，再通過(guò)一元回歸模型進(jìn)行擬合后獲?。蛔詈蠼?jīng)過(guò)滑動(dòng)推算得到含水層孔隙度n與水的體積壓縮系數(shù)β，然后再利用氣壓系數(shù)Bp和潮汐因子Bg的公式進(jìn)行推導(dǎo)計(jì)算得出固體骨架的體積壓縮系數(shù)α。

2.2 井水位變化量與含水層垂直向應(yīng)力變化量的關(guān)系

依據(jù)前人的研究結(jié)果[8-9]，在理想的水平層狀含水層（一維）模式下，井孔含水層垂直向的應(yīng)力變化量Δσz與含水層部分介質(zhì)參數(shù)、井水位變化量之間存在以下定量關(guān)系，并且能在一定程度上反映出該區(qū)域的構(gòu)造活動(dòng)特征，即

式中，Δσz為含水層垂直向應(yīng)力變化量，E為含水層固體骨架的楊氏模量（E=1/α），ΔH指剔除地下水開采、降雨影響后的含水層應(yīng)力變化引起的壓力水頭變化量，即井水位變化量。

從物理意義及異常性質(zhì)判定來(lái)講，當(dāng)井孔含水層系統(tǒng)所受應(yīng)力增強(qiáng)，即Δσz＞0 時(shí)，井水位上升，水位埋深值H變小，其變化量ΔH＜0；當(dāng)井孔含水層系統(tǒng)所受應(yīng)力減弱，即Δσz＜0 時(shí)，井水位下降，水位埋深值H變大，其變化量ΔH＞0[9]。

3 應(yīng)力變化及地震形勢(shì)分析

3.1 基于井水位變化量和含水層垂直向應(yīng)力變化量的定量分析

本文利用瓊東北地區(qū)5 口觀測(cè)井的氣壓系數(shù)和O1、M2波潮汐因子等參數(shù)，反演得到井孔含水層孔隙度、水的體積壓縮系數(shù)、含水層固體骨架的體積壓縮系數(shù)等參數(shù)，在理想承壓含水層水平層狀模式下，定量地計(jì)算了2008 年以來(lái)5 口觀測(cè)井的垂直向應(yīng)力場(chǎng)的動(dòng)態(tài)變化過(guò)程（圖6）。圖6 中各井井水位變化量和含水層垂直向應(yīng)力變化量的時(shí)序月值曲線顯示，井孔含水層系統(tǒng)所受應(yīng)力增強(qiáng)時(shí)即Δσz＞0，ΔH＜0。5 口井水位變化量ΔH＜0 時(shí)，其相應(yīng)的各井含水層垂直向應(yīng)力變化量Δσz＞0 對(duì)應(yīng)地很好。

圖6 ?？赯K26 井、瓊海加積井、向榮村井、火山流體井、文昌潭牛井水位變化量與含水層垂直向應(yīng)力變化量月值曲線Fig.6 Monthly curve of water level variation and vertical stress variation of aquifer Haikou ZK26 well，Qionghai Jiaji well，Xiangrong village well，Volcanic fluid well and Wenchang Tanniu well

其中，?？赯K26 井、瓊海加積井2 口深井的含水層垂直向應(yīng)力變化量曲線形態(tài)相似，且2017 年以來(lái)兩井—含水層系統(tǒng)應(yīng)力明顯較2017 年前減弱；向榮村井、火山流體井深度相差不大的2 口較深井的含水層垂直向應(yīng)力變化量曲線形態(tài)也相似，井—含水層系統(tǒng)以壓應(yīng)力為主的變化相較?？赯K26 井、瓊海加積井明顯較弱，且2016 年以來(lái)兩井—含水層系統(tǒng)應(yīng)力變化同樣呈減弱狀態(tài)；文昌潭牛井，2018 年以前，井—含水層系統(tǒng)以壓應(yīng)力為主的應(yīng)力變化量相對(duì)集中，而2018 年以后較零散減弱。綜合分析，基于數(shù)字化水位的井水位變化量與含水層垂直向應(yīng)力變化量的研究結(jié)果表明，瓊東北地區(qū)近幾年來(lái)應(yīng)力變化活動(dòng)水平相對(duì)減弱，需繼續(xù)跟蹤其轉(zhuǎn)折變化等情況的發(fā)生。

3.2 基于重力場(chǎng)累積動(dòng)態(tài)及GPS 基線的區(qū)域構(gòu)造活動(dòng)分析

基于井水位變化量與含水層垂直向應(yīng)力變化量的應(yīng)力場(chǎng)時(shí)序變化特征研究結(jié)果，表明瓊東北地區(qū)近幾年來(lái)應(yīng)力變化活動(dòng)水平相對(duì)較弱，下面將結(jié)合本區(qū)重力場(chǎng)累積動(dòng)態(tài)變化以及基于GPS 基線觀測(cè)數(shù)據(jù)的區(qū)域構(gòu)造活動(dòng)的多手段綜合分析，進(jìn)一步增強(qiáng)對(duì)基于井水位變化量與含水層垂直向應(yīng)力變化量研究結(jié)果的可參考性。

累積動(dòng)態(tài)變化圖像，表示相對(duì)某一期或某一基準(zhǔn)的重力變化，其可突出觀測(cè)區(qū)域不同時(shí)期重力場(chǎng)動(dòng)態(tài)演化的相對(duì)狀態(tài)或累積信息。以2016 年9 月這一期為基準(zhǔn)，計(jì)算2016-09—2017-10（累積1 年）、2016-09—2018-10（累積2 年）、2016-09—2019-08（累積3 年）和2016-09—2020-09（累積4 年）累積重力變化：海南島陸地區(qū)4 年期間的重力場(chǎng)變化相對(duì)平穩(wěn)，總體呈東北部負(fù)變化，西南部正變化特征；海南島陸重力場(chǎng)變化較大的區(qū)域分布于瓊西南。海口ZK26 井、向榮村井、火山流體井位于海南島陸北部，瓊海加積井、文昌潭牛井位于海南島陸東部，2016 年以來(lái)海南島陸北部和東部地區(qū)重力場(chǎng)累積變化比較平穩(wěn)，無(wú)顯著異常（圖7）。

圖7 海南島陸重力場(chǎng)累積變化動(dòng)態(tài)分析圖像Fig.7 Dynamic analysis image of cumulative change of land gravity field in Hainan Island

GPS 觀測(cè)數(shù)據(jù)可以反映地球內(nèi)部應(yīng)力場(chǎng)變化在地面的形變響應(yīng)分布，及構(gòu)造應(yīng)力緩慢作用過(guò)程中地殼的某些形變和運(yùn)動(dòng)特征。因此，我們通過(guò)利用GPS 資料來(lái)分析研究區(qū)域地殼形變特征以及對(duì)區(qū)域地下流體變化的影響。

目前海南在瓊北地區(qū)活動(dòng)斷裂帶兩側(cè)布設(shè)有GPS 流動(dòng)觀測(cè)站（圖8a），以期通過(guò)斷裂兩側(cè)的GPS流動(dòng)觀測(cè)點(diǎn)對(duì)斷裂活動(dòng)進(jìn)行形變觀測(cè)，從而監(jiān)測(cè)區(qū)域內(nèi)火山及地震活動(dòng)狀態(tài)。海口ZK26 井、火山流體井、向榮村井、瓊海加積井、文昌潭牛井均位于該監(jiān)測(cè)區(qū)域周邊，根據(jù)瓊北地區(qū)流動(dòng)GPS 的解算結(jié)果（圖8b）可知，2008—2020 年以來(lái)參考站馬鞍嶺與各流動(dòng)觀測(cè)站之間的6 條基線逐年變化基本在15 mm以內(nèi)，總體變化較小且穩(wěn)定，即瓊北及周邊地區(qū)沒(méi)有明顯的拉張和擠壓趨勢(shì)，表明塊體內(nèi)部較穩(wěn)定。

圖8 （a） 瓊北地區(qū)流動(dòng)GPS 觀測(cè)點(diǎn)分布；（b） 馬鞍嶺參考站與各流動(dòng)站點(diǎn)基線逐年變化示意圖Fig.8 （a） Distribution of mobile GPS observation points in northern Hainan；（b） Schematic diagram of annual change of baseline of Ma’anling reference station and mobile stations

4 結(jié)論與討論

瓊東北地區(qū)作為海南島陸主要活動(dòng)構(gòu)造帶的集中區(qū)域，是海南地震監(jiān)測(cè)研究的重點(diǎn)區(qū)域。本文借鑒前人的研究結(jié)果，利用5 口井的氣壓系數(shù)和各井水位的主要分波O1、M2潮汐因子等，滑動(dòng)擬合得到不排水狀態(tài)下，各井部分含水層參數(shù)。接著在理想水平層狀含水層模式下，利用這些參數(shù)與井水位變化量、含水層垂直向應(yīng)力變化量的關(guān)系式，定量分析了瓊東北地區(qū)構(gòu)造應(yīng)力場(chǎng)的時(shí)序變化特征?；跀?shù)字化水位的井水位變化量與含水層垂直向應(yīng)力變化量研究結(jié)果，再結(jié)合區(qū)域重力場(chǎng)累積動(dòng)態(tài)變化與基于GPS 基線觀測(cè)數(shù)據(jù)構(gòu)造活動(dòng)的綜合分析，表明瓊東北地區(qū)近5 年來(lái)應(yīng)力變化活動(dòng)總體呈減弱狀態(tài)，需跟蹤其轉(zhuǎn)折變化情況的發(fā)生。

同時(shí)發(fā)現(xiàn)瓊東北地區(qū)井孔性質(zhì)相當(dāng)?shù)木?，主要表現(xiàn)為井深差不多的井孔，其井水位變化量與含水層垂直向應(yīng)力變化量曲線形態(tài)相似。隨著本區(qū)域井孔觀測(cè)數(shù)據(jù)的積累，將繼續(xù)對(duì)比其他井孔的應(yīng)力變化定量分析結(jié)果，進(jìn)行更為全面的論證。含水層介質(zhì)構(gòu)造應(yīng)力變化的定量分析，對(duì)于有效利用數(shù)字化水位資料來(lái)甄別地震地球物理異常和機(jī)理分析具有重要的意義[8]，因此，本文的研究結(jié)果將積極運(yùn)用到海南省日常的會(huì)商及異常核實(shí)等分析工作中，為本區(qū)域地震形勢(shì)的分析研判提供有力依據(jù)。

與傳統(tǒng)的通過(guò)現(xiàn)場(chǎng)抽水試驗(yàn)和室內(nèi)實(shí)驗(yàn)等不同，利用數(shù)字化水位等資料，結(jié)合氣壓系數(shù)和維尼迪科夫潮汐調(diào)和分析結(jié)果，而獲取含水層介質(zhì)的孔隙度、固體骨架的體積壓縮系數(shù)和水的體積壓縮系數(shù)是簡(jiǎn)便易行的[8]。相較定性的分析，經(jīng)參與了含水層垂直向應(yīng)力變化量的分析而計(jì)算得到的各井含水層的孔隙度、固體骨架的體積壓縮系數(shù)和水的體積壓縮系數(shù)是動(dòng)態(tài)變化的，因此，得到的瓊東北地區(qū)構(gòu)造應(yīng)力場(chǎng)的變化結(jié)果是可靠、合理的[8]。另，本文均是在假定理想的水平層狀含水層，介質(zhì)是各向同性、均質(zhì)的彈性體模型情況下完成的，實(shí)際應(yīng)用過(guò)程中還需深入研究。

致謝

在本文的研究過(guò)程中，非常感謝內(nèi)蒙古自治區(qū)地震局丁風(fēng)和高級(jí)工程師提供的幫助！