基于數字化水位的瓊東北地區構造應力場時序特征分析＊

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

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

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

引言

區域地殼構造應力場的分析是一項較困難復雜的事情，通常通過測量應變實現應力的測量，其耗資較高。研究表明，承壓含水層井水位的變化反映了含水層孔隙壓力的變化，而含水層孔隙壓力的變化與含水層所受的應力狀態有密切的關系[1]。測量地殼深部彈性應力應變的方法有觀測記錄斷層活動區附近鉆井內流體的壓力變化，而測量鉆井水位的變化是研究流體壓力變化的方法之一[2]。流動性和難壓縮性是地殼結構中活躍組分之一的地下水所存在的主要性狀，而當地下水在含水層中形成一個封閉條件較好的承壓系統時，其水位的升降變化就能較客觀地反映出周邊區域地殼應力應變的狀態。構造應力變化將導致孔隙體積的變化，從而造成孔隙壓力的改變；井孔水位的變化量與巖石中的孔隙壓力有關，進而與構造應力有關。因此，利用已有承壓井水位連續觀測資料，研究區域現今構造應力場狀態將是一條既有效又經濟的途徑[2]。

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

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

1 資料的分析與處理

1.1 觀測井基本概況

本文選取了位于瓊東北地區的5 口數字水位觀測井，其中，海口ZK26 井、向榮村井、火山流體井，這3 口井均沿馬裊—鋪前斷裂帶分布，且為活斷裂，井距離海岸較近，井水位效應表現為海潮潮汐顯著、氣壓潮汐次之、固體潮汐最弱的復合潮汐形態；文昌潭牛井和瓊海加積井，主要沿文昌—瓊?！齺啍嗔褞Х植?，這2 口井水位潮汐效應表現為以固體潮汐顯著、氣壓潮汐次之、海潮潮汐最弱的復合潮汐形態（圖1）。各井承壓性均較好，觀測資料相對較長且連續穩定。表1 為各井經緯度、井深、含水層巖性、地下水類型、水位埋深等基本情況。

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

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

1.2 數據的整理與處理

本研究從海南省地球物理臺網中心數據庫下載井孔數字化水位和氣壓的整點值數據，選取研究資料為2008—2020 年，其中，文昌潭牛井因數據庫故障原因，2008—2012 年氣壓資料缺失，故該井研究起始時間為2013 年。首先將水位、氣壓兩測項數據進行逐月檢查，其存在的突跳或錯誤數據需做刪除處理，并按999999 進行缺數標記。水位數據整理成以米為單位的數據文件。

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

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

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

各井氣壓系數的主要獲取步驟為：首先對水位和氣壓數據進行別爾采夫濾波，其次對濾波后的數據進行一階差分處理，最后根據差分結果，構建一元回歸模型來擬合氣壓系數[8]。其大小因井而異，變化范圍一般為1—10 mm/hPa，個別井特殊時段下氣壓系數值會出現負值或較大值，本文獲取的各井氣壓系數值基本處于此合理變化范圍之內。圖4 為各井氣壓系數月值曲線。

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

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

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

經過計算得出的潮汐因子、氣壓系數和含水層參數等，均為月值數據，基本能滿足地震的中短期預測對資料間隔的需求。

2 基本理論及參數求取過程

2.1 含水層參數的反演計算

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

上兩式聯合可得到：

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

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

2.2 井水位變化量與含水層垂直向應力變化量的關系

依據前人的研究結果[8-9]，在理想的水平層狀含水層（一維）模式下，井孔含水層垂直向的應力變化量Δσz與含水層部分介質參數、井水位變化量之間存在以下定量關系，并且能在一定程度上反映出該區域的構造活動特征，即

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

從物理意義及異常性質判定來講，當井孔含水層系統所受應力增強，即Δσz＞0 時，井水位上升，水位埋深值H變小，其變化量ΔH＜0；當井孔含水層系統所受應力減弱，即Δσz＜0 時，井水位下降，水位埋深值H變大，其變化量ΔH＞0[9]。

3 應力變化及地震形勢分析

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

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

圖6 海口ZK26 井、瓊海加積井、向榮村井、火山流體井、文昌潭牛井水位變化量與含水層垂直向應力變化量月值曲線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 口深井的含水層垂直向應力變化量曲線形態相似，且2017 年以來兩井—含水層系統應力明顯較2017 年前減弱；向榮村井、火山流體井深度相差不大的2 口較深井的含水層垂直向應力變化量曲線形態也相似，井—含水層系統以壓應力為主的變化相較海口ZK26 井、瓊海加積井明顯較弱，且2016 年以來兩井—含水層系統應力變化同樣呈減弱狀態；文昌潭牛井，2018 年以前，井—含水層系統以壓應力為主的應力變化量相對集中，而2018 年以后較零散減弱。綜合分析，基于數字化水位的井水位變化量與含水層垂直向應力變化量的研究結果表明，瓊東北地區近幾年來應力變化活動水平相對減弱，需繼續跟蹤其轉折變化等情況的發生。

3.2 基于重力場累積動態及GPS 基線的區域構造活動分析

基于井水位變化量與含水層垂直向應力變化量的應力場時序變化特征研究結果，表明瓊東北地區近幾年來應力變化活動水平相對較弱，下面將結合本區重力場累積動態變化以及基于GPS 基線觀測數據的區域構造活動的多手段綜合分析，進一步增強對基于井水位變化量與含水層垂直向應力變化量研究結果的可參考性。

累積動態變化圖像，表示相對某一期或某一基準的重力變化，其可突出觀測區域不同時期重力場動態演化的相對狀態或累積信息。以2016 年9 月這一期為基準，計算2016-09—2017-10（累積1 年）、2016-09—2018-10（累積2 年）、2016-09—2019-08（累積3 年）和2016-09—2020-09（累積4 年）累積重力變化：海南島陸地區4 年期間的重力場變化相對平穩，總體呈東北部負變化，西南部正變化特征；海南島陸重力場變化較大的區域分布于瓊西南。海口ZK26 井、向榮村井、火山流體井位于海南島陸北部，瓊海加積井、文昌潭牛井位于海南島陸東部，2016 年以來海南島陸北部和東部地區重力場累積變化比較平穩，無顯著異常（圖7）。

圖7 海南島陸重力場累積變化動態分析圖像Fig.7 Dynamic analysis image of cumulative change of land gravity field in Hainan Island

GPS 觀測數據可以反映地球內部應力場變化在地面的形變響應分布，及構造應力緩慢作用過程中地殼的某些形變和運動特征。因此，我們通過利用GPS 資料來分析研究區域地殼形變特征以及對區域地下流體變化的影響。

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

圖8 （a） 瓊北地區流動GPS 觀測點分布；（b） 馬鞍嶺參考站與各流動站點基線逐年變化示意圖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 結論與討論

瓊東北地區作為海南島陸主要活動構造帶的集中區域，是海南地震監測研究的重點區域。本文借鑒前人的研究結果，利用5 口井的氣壓系數和各井水位的主要分波O1、M2潮汐因子等，滑動擬合得到不排水狀態下，各井部分含水層參數。接著在理想水平層狀含水層模式下，利用這些參數與井水位變化量、含水層垂直向應力變化量的關系式，定量分析了瓊東北地區構造應力場的時序變化特征?；跀底只坏木蛔兓颗c含水層垂直向應力變化量研究結果，再結合區域重力場累積動態變化與基于GPS 基線觀測數據構造活動的綜合分析，表明瓊東北地區近5 年來應力變化活動總體呈減弱狀態，需跟蹤其轉折變化情況的發生。

同時發現瓊東北地區井孔性質相當的井，主要表現為井深差不多的井孔，其井水位變化量與含水層垂直向應力變化量曲線形態相似。隨著本區域井孔觀測數據的積累，將繼續對比其他井孔的應力變化定量分析結果，進行更為全面的論證。含水層介質構造應力變化的定量分析，對于有效利用數字化水位資料來甄別地震地球物理異常和機理分析具有重要的意義[8]，因此，本文的研究結果將積極運用到海南省日常的會商及異常核實等分析工作中，為本區域地震形勢的分析研判提供有力依據。

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

致謝

在本文的研究過程中，非常感謝內蒙古自治區地震局丁風和高級工程師提供的幫助！