張淑亮 李宏偉 呂 芳 陳 慧

（山西省地震局，太原 030021）

引言

靜樂井水位是山西乃至華北地區(qū)映震能力較強的地下流體觀測井之一，井水位異常能客觀反映區(qū)域應(yīng)力場的變化，全球7 級以上地震均能記錄到水震波。1983 年始測以來出現(xiàn)的幾次高值異常與山西北部至晉冀蒙交界地區(qū)強震活動時間具有準同步性，在1989 年大同-陽高Ms6.1 級地震、1991 年忻州Ms5.1 級地震、1991 年大同-陽高Ms5.8 級地震、1998 年張北Ms6.2 級地震、1999 年大同-陽高Ms5.6 級等中強地震發(fā)生前有較好的異常顯示。因此將靜樂井水位高值異常作為晉冀蒙交界地區(qū)地震趨勢預(yù)測的重要判據(jù)（張淑亮等，1997；車用太等，1999；王愛英等，1998）。 2017 年8 月靜樂井水位又出現(xiàn)大于多年平均值的高值異常，到目前為止異常仍在持續(xù)。按照以往震例時間指示意義，目前已接近發(fā)震時間，未來晉冀蒙交界區(qū)震情形勢如何發(fā)展？目前的異常是否真實反映了構(gòu)造活動加劇所引起的區(qū)域應(yīng)力場和深部地球物理場的改變？是否與中強震孕育過程有關(guān)？靜樂井水位目前的異常特征不足以解析異常區(qū)的應(yīng)力狀態(tài)，從而也直接影響對晉冀蒙交界地區(qū)地震形勢的認識和未來地震趨勢的判定。

多年的觀測表明，影響靜樂井水位動態(tài)變化的主要干擾因素為短時強降雨與河流水體的荷載作用。無論是短時間的強降雨還是洪水期河流的荷載作用，均可在短時間內(nèi)造成井水位的快速變化（張淑亮等，2005；車用太等，2004；魚金子等，1994；張昭棟等，1989；張昭棟等，1990）。目前對這2 種干擾因素的識別與排除基本上以定性分析為主，定量分析方法較少。特別是對于同時受2 種因素影響的靜樂井水位而言，在異常分析與判定工作中，僅采用與同期降雨量對比的方法分析井水位變化與降雨量的關(guān)系，而河流水位變化對井水位的影響分析工作至今尚未開展。為此，本文利用數(shù)值模擬與含水層垂向應(yīng)力反演的方法，對降雨與河流荷載作用下的靜樂井區(qū)應(yīng)力進行定量計算，并從區(qū)域應(yīng)力場的角度分析2 種方法計算結(jié)果的差異，以期為地震趨勢判定提供較可靠的預(yù)測依據(jù)。

1 井區(qū)概況

靜樂流體觀測井位于呂梁山斷塊隆起區(qū)的靜樂向斜，位于東碾河斷裂帶上，地處東碾河南岸，與東碾河相距約130m，其北部和南部均為山區(qū)（見圖1）。靜樂井深362.92m，井區(qū)為河谷地區(qū)，主要發(fā)育有下古生界灰?guī)r、白云巖喀斯特含水層及第四系沖洪積砂礫石含水層，2 個含水層間有較強的水力聯(lián)系。觀測層巖性時代為中奧陶系灰?guī)r裂隙溶洞水，其頂板埋深為46.3m，單位涌水量84.5L/（s·m），滲透系數(shù)為2.41—11.5m/d。觀測層具有厚度大、透水性強的特點。含水層巖石溶洞發(fā)育，靜樂井在230—290m 處通過破碎帶，地下水由北西經(jīng)井區(qū)向東南方向排泄，井區(qū)處于強徑流區(qū)，地下水補給主要來自北、南山區(qū)基巖裸露區(qū)的大氣降水補給和東碾河河水沿斷裂破碎帶入滲補給。

圖1 井區(qū)地質(zhì)構(gòu)造 Fig. 1 Well area tectonic map

2 強降雨與河流荷載作用下靜樂井水位動態(tài)特征

多年的觀測表明，大于15mm 的集中降雨將引起靜樂井水位的快速上升（張淑亮等，2005），可能由于降雨在地表形成荷載，使得含水層受力狀態(tài)發(fā)生變化，引起水井水位上升。東碾河水位變化是影響靜樂井水位變化的又一重要因素。當大面積洪水作用通過井區(qū)時，井水位迅速大幅上升，幾乎沒有時間滯后，洪水期水位上升，洪水過后水位下降，水位升降幅度取決于洪水深度。

圖2 靜樂井水位與降雨量和東碾河流量的對比 Fig. 2 Comparison of Jingle well level with rainfall and Dongnian river discharge

圖2 所示為靜樂井水位與井區(qū)降雨量和東碾河流量的對比。由圖2 可知，2017 年8 月20 日左右單日降雨量為82.8mm 時，東碾河流量由原來的0.777m3/s 增至9.8m3/s，增加了9.023m3/s； 靜樂井水位由9.074m 升至7.747m，水位上升幅度達1.327m。由此可知，降雨與河流荷載作用是影響靜樂井水位動態(tài)變化的主要因素，但降雨結(jié)束后井水位并未迅速回落至上升前的水平。張淑亮等（1997）曾對靜樂井水位高值變化與降雨的關(guān)系進行了較系統(tǒng)的研究，認為降雨量的增大雖可在短時間內(nèi)造成井水位升降變化，但不能改變水位的年變動態(tài)，只有在區(qū)域應(yīng)力場增強的背景條件下，降雨引起的附加力源才可能造成井水位年動態(tài)的異常變化，降雨荷載是區(qū)域應(yīng)力場變化的誘發(fā)因素，也可能是在2017 年雨季結(jié)束后井水位仍處于高值的主要原因。

3 降雨與河流荷載作用下井區(qū)應(yīng)力數(shù)值模擬分析

利用井區(qū)周圍地殼結(jié)構(gòu)、不同地層介質(zhì)彈性參數(shù)建立三維地質(zhì)模型，預(yù)測降雨與河流荷載作用下的井區(qū)應(yīng)力。由于本文基于數(shù)值模擬與含水層垂向應(yīng)力反演的方法進行定量計算，因此，依據(jù)靜樂井井區(qū)地層巖性特征和井-含水層特征參數(shù)選取計算參數(shù)。因靜樂井位于河谷地區(qū)，故計算時以變質(zhì)巖系中大理巖力學性質(zhì)指標經(jīng)驗數(shù)據(jù)為依據(jù)。

3.1 三維地質(zhì)模型建立

根據(jù)靜樂井區(qū)地層結(jié)構(gòu)、河流位置及研究目標，設(shè)定有限元尺寸為1km×1km×360m，河流寬76m。采用SOLID185 單元，模型采用由面拉伸成體的建模方式，取井孔周圍圍巖彈性模量為5×1010MPa，泊松比為0.25，密度為2600kg/m3。

網(wǎng)格劃分的目的是對模型實現(xiàn)離散化，并用適當數(shù)量的網(wǎng)格單元得到最精確的解。本文選用Workblench 中的自動網(wǎng)格劃分工具，以滿足模型要求。網(wǎng)格疏密程度直接影響計算結(jié)果的精度，但網(wǎng)格加密將增加CPU 計算時間，且需更大的存儲空間。為保證模型網(wǎng)格密度和計算精度，本文選擇的網(wǎng)格尺寸為1.5m，同時適當增加井區(qū)附近區(qū)域的網(wǎng)格密度。其他參數(shù)均采用默認值，可滿足模型的需求。

載荷和約束是ANSYS 軟件求解計算的邊界條件，由所選單元的自由度形勢定義。由于本文主要分析河流與降雨荷載的影響，因此選擇垂向應(yīng)力加載。為防止模型移動，模型底端與側(cè)面采用固定約束，限制發(fā)生移動和轉(zhuǎn)動。

考慮數(shù)值模擬結(jié)果的多樣性及模型造成的誤差，建模時選擇不同的邊界條件，并反復測試計算，以可反映降雨與最大河流荷載作用對井區(qū)的最大影響范圍作為模型的邊界條件。由于降雨與河流荷載的作用主要表現(xiàn)為垂向應(yīng)力，因此模型初始地應(yīng)力只考慮重力。因本文通過數(shù)值模擬與反演的方法分析靜樂井水位在降雨與河流荷載作用下的垂向應(yīng)力在量級上是否具有相似性，進而確定影響水位變化的主要因素，對計算精度的要求相對較低，故選取物理力學參數(shù)時以研究區(qū)內(nèi)地層巖性的平均值作為模型計算參數(shù)。

3.2 結(jié)果分析

由于2017 年是近年來靜樂井水位年變幅最大的年份，也是降雨與河流荷載影響最大的年份，故本文以2017 年8 月洪水期東碾河水位、靜樂井區(qū)降雨量和靜樂井水位為研究對象，計算降雨和河流荷載作用下的靜樂井井區(qū)的垂向位移和垂向應(yīng)力。

（1）單降雨荷載的影響

以2017 年8 月降雨量作為加載體，將降雨量等效為荷載，計算得到靜樂井垂向位移和垂向應(yīng)力如圖3 所示。由圖3 可知，降雨荷載對靜樂井產(chǎn)生的垂向位移約為0.6mm，垂向應(yīng)力約為12hPa。

（2）單河流荷載的影響

以洪水期東碾河水位最大變化量作為加載體，并將河水水位等效為荷載計算得到靜樂井垂向位移和垂向應(yīng)力如圖3 所示。由圖3 可知，洪水期河流荷載對靜樂井產(chǎn)生的垂向位移約為0.3mm，垂向應(yīng)力約為15hPa。

圖3 單降雨荷載作用下靜樂井垂向位移和垂向應(yīng)力 Fig.3 Vertical displacement and vertical stress of Jingle well under single rainfall load

圖4 單河流荷載作用下靜樂井垂向位移和垂向應(yīng)力 Fig.4 Vertical displacement and vertical stress of Jingle well under single river load

（3）降雨與河流荷載的綜合影響

由于靜樂井水位變化是在因降雨與河流荷載的共同作用，單河流和單降雨的影響僅能反映井水位的部分變化，因此，需考慮二者的共同作用。以2017 年8 月東碾河水位最大變化量和降雨量作為共同加載體，并將二者等效為荷載，計算得到靜樂井垂向位移和垂向應(yīng)力如圖5 所示。由圖5 可知，降雨和河流荷載共同作用對靜樂井產(chǎn)生的垂向位移約為0.7mm，垂向應(yīng)力約為30hPa，與單河流和單降雨影響的二者之和基本一致，表明降雨與河流荷載作用是影響靜樂井水位變化的主要因素。

圖5 河流與降雨荷載共同作用下靜樂井垂向位移和垂向應(yīng)力 Fig.5 Vertical displacement and vertical stress of Jingle well under combined action of river and rainfall loads

4 靜樂井-含水層垂向應(yīng)力反演

數(shù)值模擬結(jié)果表明，降雨和河流荷載是影響靜樂井水位變化的主要因素。為進一步驗證數(shù)值模擬結(jié)果能否真實反映井區(qū)實際應(yīng)力變化，還需將其與實際水位觀測值所反映的井區(qū)應(yīng)力進行對比。在由井水位變化反演井區(qū)應(yīng)力場變化方面，目前取得了一些進展，如黃輔瓊等（2004）利用華北地區(qū)40 多個深井水位動態(tài)變化資料，探討研究了華北地區(qū)現(xiàn)今構(gòu)造應(yīng)力場狀態(tài)；孫小龍等（2011）運用小波分析法提取出能反映水位多年動態(tài)變化的趨勢信息，反演出華北地區(qū)多年構(gòu)造應(yīng)力場的變化特征；丁風和等（2015）利用部分含水層介質(zhì)參數(shù)、井水位變化量與含水層垂向應(yīng)力變化量的關(guān)系式，定量分析了張渤帶地區(qū)構(gòu)造應(yīng)力場的動態(tài)變化過程。

本文根據(jù)靜樂井-含水層介質(zhì)參數(shù)及井水位變化量等數(shù)據(jù)，反演井水位變化引起的含水層垂向應(yīng)力。應(yīng)力反演所用數(shù)據(jù)時段為2016—2017 年（根據(jù)東碾河水位與流量始測時間確定）。反演前首先對觀測數(shù)據(jù)進行預(yù)處理（丁風和等，2015），去除長周期干擾因素、降雨、氣壓等年際變化對水位的影響；然后采用去趨勢法和傅里葉周期法，剔除地下水開采和降水補給等長周期干擾對水位的影響。

4.1 反演所需各類參數(shù)計算

（1）氣壓系數(shù)與潮汐效率計算

利用卡爾曼濾波方法消除地球固體潮對所選時段水位、氣壓整點值數(shù)據(jù)的影響。利用濾波后的日均值數(shù)據(jù)，采用高階差分和一元線性回歸等方法獲得靜樂井氣壓系數(shù)。利用維尼迪科夫潮汐調(diào)和分析方法獲得靜樂井水位M2波潮汐效率。氣壓系數(shù)和潮汐效率計算結(jié)果如表1 所示。

表1 靜樂井-含水層多種參數(shù)計算結(jié)果 Table 1 Calculation results of various parameters of Jingle well-aquifer

（2）靜樂井-含水層介質(zhì)參數(shù)計算

不排水狀態(tài)下，井水位氣壓系數(shù)和潮汐效率可分別表示為（Bredehoeft，1967；李春洪等，1990；張昭棟等，1993）：

式中BP為井水位氣壓系數(shù)；Bg為M2波潮汐效率；n為含水層孔隙度；α為含水層固體骨架體積壓縮系數(shù)；β為含水層內(nèi)水的體積壓縮系數(shù)；ρ為水的密度；g為重力加速度（ρg=0.098hPa/mm）。因此，n、β可依據(jù)式（3）與表1 中的潮汐效率和氣壓系數(shù)得到。根據(jù)式（1）或式（2）可求α。

4.2 靜樂井-含水層垂向應(yīng)力反演

利用表1 數(shù)據(jù)，根據(jù)水平層狀含水層模式下，含水層介質(zhì)參數(shù)（孔隙度、水和固體骨架體積壓縮系數(shù)等）、井水位變化量與含水層垂向應(yīng)力變化量的關(guān)系式（張昭棟等，1987）等來計算含水層垂向應(yīng)力，含水層垂向應(yīng)力計算如下：

式中 Δσz為含水層垂向應(yīng)力變化量；E為含水層固體骨架的楊氏模量（E= 1/α）；△H為剔 除地下水開采、降雨和氣壓影響后的含水層應(yīng)力變化引起的壓力水頭變化量，即井水位變化 量。當井-含水層系統(tǒng)所受應(yīng)力增強，即 Δσz>0 時，井水位上升，水位埋深H變小，其變化量ΔH<0；當井-含水層系統(tǒng)所受應(yīng)力減弱，即 Δσz<0 時，井水位下降，水位埋深H變 大，其變化量ΔH＞0。利用式（4）與表1 中有關(guān)參數(shù)計算靜樂井2016—2017 年含水層垂向應(yīng)力，結(jié)果如表1、圖6 所示?？芍o樂井含水層垂向應(yīng)力高值時段與降雨量高值時段、東碾河流量高值時段較吻合，再次表明降雨與河流荷載作用是造成含水層垂向應(yīng)力高值變化的主要因素，井水位的大幅突升變化與垂向應(yīng)力增大密切相關(guān)。

5 靜樂井水位異常分析

5.1 模型預(yù)測值與井水位反演值對比分析

由數(shù)值模擬結(jié)果可知，無論單河流荷載還是單降雨荷載效應(yīng)均對洪水期靜樂井水位造成了一定影響，二者造成的垂向位移累計0.9mm，垂向應(yīng)力累計27hPa。河流和降雨綜合影響模型結(jié)果也得到類似結(jié)果（垂向位移累計0.7mm，垂向應(yīng)力累計30hPa）。由含水層垂向應(yīng)力反演結(jié)果可知，洪水期含水層的應(yīng)力顯著增強，由正常時段的4.2hPa 增至79hPa 左右。2 種方法計算結(jié)果均表明靜樂井水位上升變化與河流荷載效應(yīng)和降雨荷載效應(yīng)密切相關(guān)。但2 種方法計算結(jié)果存在一定差異，降雨與河流荷載作用下井區(qū)垂向應(yīng)力數(shù)值模擬結(jié)果小于含水層垂向應(yīng)力反演值，相差39hPa。

圖6 東碾河流量、靜樂井-含水層垂向應(yīng)力與降雨量的對比 Fig. 6 Flow rate of Dongnian river, vertical stress of Jingle well-aquifer and rainfall comparison chart

由于數(shù)值模擬計算得到的垂向應(yīng)力是研究區(qū)自身重力、降雨與河流荷載作用下產(chǎn)生的垂向應(yīng)力之和，而由井水位反演的垂向應(yīng)力不僅包含上述因素產(chǎn)生的垂向應(yīng)力，可能還含有一些因構(gòu)造活動引起的垂向應(yīng)力變化。事實上2 種方法的計算結(jié)果的確存在一定 差異。雖然2 種計算方法均受建模條件和參數(shù)選取的影響而存在計算誤差，但根據(jù)以往的研究，誤差值一般不會大于計算值，因此2 種計算結(jié)果的差異性是客觀存在的。按照表1 中靜樂井水位平均氣壓系數(shù)4.3521mm/hPa 可推算出2 種方法垂向應(yīng)力計算結(jié)果差值引起的井水位變化幅度約為16.97cm，大于正常時段水位平均月變化幅度0.9cm 的水平。表明影響靜樂井水位快速上升的因素除降雨與河流荷載作用外，可能還受其他因素的影響。

5.2 水位高值異常與構(gòu)造活動的關(guān)系

為進一步探討靜樂井水位2017年8月以來的高值異常除受降雨與河流荷載影響外是否受構(gòu)造活動的影響，本文從可反映區(qū)域應(yīng)力狀態(tài)變化的參數(shù)和測量值（缺震和b 值、地震矩釋放、GPS 基線等）及同一構(gòu)造區(qū)其他前兆測項變化特征與靜樂井水位高值異常進行對比分析。

大量研究結(jié)果表明，b 值與其對應(yīng)區(qū)域應(yīng)力狀態(tài)、地殼破裂強度有關(guān)（劉雁冰等，2017；朱艾斕等，2005；吳小平，1990）。1983 年以來山西北部至晉冀蒙交界區(qū)3 級以上地震缺震和b 值曲線有4 次異常過程（見圖7），第1 次異常出現(xiàn)在1988 年，異常持續(xù)過程中在山西地震帶出現(xiàn)一組中強地震活動，先后發(fā)生了大同-陽高6.1 級地震、侯馬4.9 級和忻州5.1 級地震、大同-陽高5.8 級地震等，其中大同-陽高6.1 級地震被認為是華北地區(qū)在經(jīng)歷唐山大地震后進入新活躍幕的標志，是華北應(yīng)力場趨于增強的結(jié)果；第2 次異常出現(xiàn)在1995 年初，異常持續(xù)過程中先后發(fā)生了1996 年5 月包頭6.4 級地震、1998 年1 月張北6.2 級地震、1999年3 月張北5.6 級地震、1999 年11 月大同-陽高5.6 級地震等中強地震活動，既是1989 年開始的華北新一輪活躍幕的繼續(xù)，又是華北應(yīng)力場增強的結(jié)果；第3 次異常始于2007 年，這次異常過程中在山西北部至晉冀蒙交界區(qū)發(fā)生了幾次4 級以上中等地震，部分學者認為，這組中等地震的發(fā)生是受汶川地震的影響，鄂爾多斯塊體與華北平原塊體相對擠壓和扭錯增強，導致山西帶形變場與構(gòu)造應(yīng)力場由原來的構(gòu)造拉張轉(zhuǎn)為構(gòu)造擠壓，是汶川地震延遲觸發(fā)的結(jié)果 （劉峽等，2013，朱艾斕等，2010；劉瑞春等，2014）；2017 年又出現(xiàn)同步異常變化，目前異常仍在發(fā)展，且缺震和b 值的幾次異常與靜樂井水位高值異常較為同步，因此，可認為靜樂井水位高值異常反映的是華北區(qū)域應(yīng)力場增強的過程，是山西地震帶地震活動增強的重要判據(jù)（張淑亮等，1997）。2017 年靜樂井水位高值異常也可能與區(qū)域應(yīng)力場增強有一定關(guān)系。

圖7 靜樂井水位月均值與山西地震帶3 級地震缺震和b 值的對比 Fig. 7 Comparison chart of monthly average water level of jingle well with magnitude 3 earthquake deficiency and b value in Shanxi seismic belt

由跨過靜樂井的山西岢嵐-山西太原GPS 基線時序曲線圖8 可知，2016 年以來基線由原來的下降趨勢（基線縮短）變?yōu)榫徛仙ɑ€拉長），應(yīng)力狀態(tài)在2016 年發(fā)生了變化，由靜樂井所在區(qū)域地震矩加速釋放時圖9 可知，2017 年井區(qū)存在明顯的加速釋放低值區(qū)。GPS基線和地震矩釋放均表明2017 年靜樂井所在區(qū)域應(yīng)力場在增強。

圖8 山西岢嵐-山西太原GPS 基線時序曲線 Fig.8 GPS baseline timing curve of Kelan, Shanxi-Taiyuan, Shanxi

靜樂井構(gòu)造上位于山西西部呂梁山隆起區(qū)，位于同一構(gòu)造區(qū)距靜樂井較近的前兆測點有寧武鉆孔應(yīng)變（60km）和神池鉆孔應(yīng)變（95km），在靜樂井水位突升前后這2 個測點的應(yīng)變觀測值也出現(xiàn)大幅度反向變化，所以說靜樂井水位高值變化不是獨立事件（見圖10），可能反映了所在構(gòu)造區(qū)應(yīng)力場的改變。

圖9 靜樂井所在區(qū)域地震矩加速釋放時空圖 Fig.9 Accelerated release diagram of seismic moment in jingle well area

綜上所述，靜樂井水位的高值異常除與降雨、河流荷載作用有關(guān)外，井區(qū)構(gòu)造活動增強可能也是一種因素。

圖10 靜樂井水位與同一構(gòu)造區(qū)其它前兆測項對比圖 Fig.10 Comparison between Jingle well water level and other precursor measurements in the same structure area

6 結(jié)論與認識

通過數(shù)值模擬、井-含水層垂向應(yīng)力反演、區(qū)域應(yīng)力場變化特征等對靜樂井水位2017 年以來出現(xiàn)的高值異常進行初步探討，得出以下結(jié)論：

（1）根據(jù)井區(qū)周圍地殼結(jié)構(gòu)、不同地層介質(zhì)彈性參數(shù)建立三維地質(zhì)模型，并以2017 年洪水期東碾河水位最大變化量和降雨量作為共同加載體，模擬得到靜樂井垂向位移約為0.7mm，垂向應(yīng)力約為30hPa。

（2）利用靜樂井含水層孔隙度、固體骨架體積壓縮系數(shù)和水的體積壓縮系數(shù)等參數(shù)，根據(jù)水平層狀含水層模式下含水層介質(zhì)參數(shù)、井水位變化量與含水層垂向應(yīng)力變化量的關(guān)系式，得到靜樂井-含水層垂向應(yīng)力約為79hPa。

（3）數(shù)值模擬與井-含水層垂向應(yīng)力反演結(jié)果均表明，靜樂井水位高值異常與同時段降雨量增多、河流荷載效應(yīng)增強關(guān)系密切，但2 種結(jié)果存在一定差異，含水層垂向應(yīng)力反演值大于降雨與河流荷載效應(yīng)，相差39hPa，這種差異可能反映靜樂井水位除受降雨與河流荷載作用外，還可能受其他因素的影響。

（4）通過與同一構(gòu)造區(qū)距靜樂井較近的前兆測點寧武鉆孔應(yīng)變和神池鉆孔應(yīng)變對比分析可知，靜樂井水位突升前后2 個測點應(yīng)變觀測值出現(xiàn)大幅度反向變化，表明靜樂井水位高值變化不是獨立事件，可能反映了所在構(gòu)造區(qū)應(yīng)力場發(fā)生了改變。

（5）可反映區(qū)域應(yīng)力場變化特征的幾種參數(shù)對2 種方法計算結(jié)果產(chǎn)生差異的原因表明，2017 年靜樂井水位高值異常期間山西地震帶3 級地震缺震和b 值、穿過靜樂井的GPS 基線和地震矩釋放均存在顯著的異常，推測靜樂井水位高值異常除受降雨與河流荷載作用的影響外，也可能受構(gòu)造活動增強的影響。