孟建國(guó) 蔣其峰 盧雙苓 李惠玲 楊彬 周均太

1）泰安基準(zhǔn)地震臺(tái)，山東省泰安市羅漢崖路2號(hào) 271000

2）山東省地震局，濟(jì)南 250014

3）鄒城地震臺(tái)，山東濟(jì)寧 273500

0 引言

將應(yīng)變連續(xù)觀測(cè)數(shù)據(jù)應(yīng)用于地震活動(dòng)預(yù)測(cè)是自邢臺(tái)地震后在李四光倡導(dǎo)和組織下開(kāi)展起來(lái)的。國(guó)內(nèi)對(duì)洞體應(yīng)變資料的研究始于20世紀(jì)80年代，經(jīng)過(guò)近30年的發(fā)展，洞體應(yīng)變理論研究基礎(chǔ)不斷完善。目前，洞體應(yīng)變伸縮儀在我國(guó)很多臺(tái)站都有安裝，分2測(cè)向和3測(cè)向2種，得到的資料不僅可用于各方向的線(xiàn)應(yīng)變分析，還可經(jīng)換算后用于主應(yīng)變及剪切應(yīng)變的分析。對(duì)地表附近水平面內(nèi)的面應(yīng)變變化進(jìn)行監(jiān)測(cè)的測(cè)項(xiàng)除了洞體應(yīng)變，還有鉆孔應(yīng)變。兩者具有很強(qiáng)的可比性。同時(shí)，區(qū)域內(nèi)地下物質(zhì)的應(yīng)變積累與釋放直接影響了地下物質(zhì)密度的改變，可以通過(guò)重力觀測(cè)進(jìn)行驗(yàn)證（馬棟等，2013）。本文研究的內(nèi)容及方法參考了杜慧君等（1993）用Venedikov方法對(duì)應(yīng)變觀測(cè)資料的調(diào)和分析，計(jì)算結(jié)果表明，應(yīng)變觀測(cè) NS向固體潮幅度及均方差明顯優(yōu)于EW向。本研究中還應(yīng)用了蔣駿等（1993、1994）根據(jù)應(yīng)變花理論（尹祥礎(chǔ)，1991）對(duì)多方向潮汐線(xiàn)應(yīng)變組合觀測(cè)確定平面應(yīng)變狀態(tài)及提取潮汐面應(yīng)變、體應(yīng)變、平面剪切應(yīng)變等信息的方法和公式。張雁濱等（1996、1997）對(duì)中國(guó)部分潮汐應(yīng)變觀測(cè)點(diǎn)的理論應(yīng)變狀態(tài)和觀測(cè)應(yīng)變狀態(tài)的計(jì)算為本文的研究提供了很好的理論依據(jù)。

1 觀測(cè)點(diǎn)概況

泰安基準(zhǔn)地震臺(tái)地處泰山南麓、萊蕪弧形斷裂帶的北側(cè)附近，臺(tái)址基巖為花崗片麻巖，巖體完整，致密均勻，地脈動(dòng)水平低，波導(dǎo)性能良好，測(cè)量信噪比高。臺(tái)站形變觀測(cè)儀器和測(cè)震觀測(cè)儀器位于臺(tái)站的專(zhuān)用山洞內(nèi)。山洞洞室覆蓋厚度29m以上，總進(jìn)深130m，其中主洞進(jìn)深76m，高差33m，避免了洞室內(nèi)空氣流動(dòng)對(duì)儀器觀測(cè)造成的影響，室溫年變幅約0.6℃。泰安臺(tái)的洞體應(yīng)變觀測(cè)儀器為SS-Y短基線(xiàn)銦鋼伸縮儀，屬首批“九五”數(shù)字化試驗(yàn)儀器。觀測(cè)試驗(yàn)運(yùn)行始于1998年2月，1999年通過(guò)驗(yàn)收并正式運(yùn)行，于2008年12月進(jìn)行了“十五”儀器升級(jí)改造。其中伸縮儀改造升級(jí)部分為標(biāo)定器、傳感器以及數(shù)據(jù)采集器，原銦鋼基線(xiàn)未做改動(dòng)。伸縮儀布設(shè)在原水管儀的基線(xiàn)墩上，與水管基線(xiàn)并行。伸縮儀基線(xiàn)EW向長(zhǎng)10.00m，NS向長(zhǎng) 30.03m；方位角 EW 向?yàn)?124°37′，NS向?yàn)?39°48′。2008年 9月，在山洞旁邊距離伸縮儀直線(xiàn)距離10m處，安裝了 TJ-2型體積式應(yīng)變儀，探頭底部的實(shí)際埋深為74.80m。探頭與地層耦合采用水泥固結(jié)方式。

2 觀測(cè)資料評(píng)價(jià)

本文采用伸縮儀的整點(diǎn)值數(shù)據(jù)進(jìn)行計(jì)算，數(shù)據(jù)處理范圍從2009年1月～2014年6月。

調(diào)和分析方法采用Venedikov調(diào)和分析公式（張雁濱等，2000），固體潮的表達(dá)為各種諧波的綜和。即

式中，y（t）為t時(shí)的固體潮觀測(cè)值；Hn為某潮汐波的觀測(cè)振幅；ωn為某潮汐波的角頻率；φn為某潮汐波的初位相；y0（t）為固體潮在t時(shí)刻的零點(diǎn)漂移值。Venedikov調(diào)和分析是把日波、半日波群分開(kāi)。分離出不同的日波波群和半日波群，然后分別求出它們各自的振幅比和相位差以及每個(gè)潮汐波計(jì)算結(jié)果的均方差和總均方差。

公式（1）可改寫(xiě)成

式中，tj為自Tj起算的以小時(shí)表示的時(shí)間；Tj為tj為零時(shí)的時(shí)間；φn（T）j為某潮汐波在Tj時(shí)初相位。

根據(jù)濾波理論，利用偶數(shù)濾波C求出與Hncosφn（Ti）有關(guān)的數(shù)Mi，利用奇數(shù)濾波S求出與 Hnsinφn（Ti）有關(guān)的數(shù) N，而后得到未知數(shù)為 xk，yk的線(xiàn)性方程組

其中，τ＝1，2；I＝1，2，…，m；k＝1，2，…p；m ＞p

上式中未知數(shù)的個(gè)數(shù)為2p個(gè)，方程的個(gè)數(shù)為2m個(gè)，ak，βk為第k波群的起止波號(hào)。每連續(xù)48h的觀測(cè)數(shù)據(jù)，經(jīng)調(diào)和分析后，對(duì)日波和半日波群分別得出2個(gè)方程。最后用最小二乘法求解以上線(xiàn)性方程組，得到潮汐因子δk和相位差σk

通過(guò)MATLAB軟件編程，將2013年全年洞體應(yīng)變2分項(xiàng)整點(diǎn)值及體應(yīng)變?nèi)暾c(diǎn)值帶入運(yùn)算方程，分別求出各參數(shù)項(xiàng)。表1、表2分別為洞體應(yīng)變的2分項(xiàng)調(diào)和分析結(jié)果。

計(jì)算結(jié)果表明，泰安臺(tái)伸縮儀NS向潮汐因子保持在0.7001、EW向保持在0.3008范圍內(nèi)，應(yīng)變觀測(cè)NS向固體潮幅度及均方差明顯優(yōu)于EW向。在張渤帶及其鄰區(qū)洞體應(yīng)變儀測(cè)項(xiàng)潮汐因子及其誤差（馬棟等，2013）的對(duì)比分析中，36個(gè)臺(tái)數(shù)據(jù)統(tǒng)計(jì)發(fā)現(xiàn)72%的臺(tái)站潮汐因子取值為0～1和1～2的占27%、2～3的占11%；充分說(shuō)明了泰安臺(tái)儀器響應(yīng)靈敏度、洞室條件等因素都符合形變臺(tái)站儀器運(yùn)行管理規(guī)范所需要求。

表1 洞體應(yīng)變NS項(xiàng)應(yīng)變調(diào)和分析表

表2 洞體應(yīng)變EW項(xiàng)應(yīng)變調(diào)和分析表

3 應(yīng)變參數(shù)提取及對(duì)比

目前，泰安臺(tái)安裝的洞體應(yīng)變?yōu)?個(gè)分量的伸縮儀，研究人員將通過(guò)構(gòu)建計(jì)算模型提取應(yīng)變組合觀測(cè)信息（蔣駿等，1993；張雁濱等，1997），進(jìn)而對(duì)應(yīng)變參數(shù)進(jìn)行提取解算。為分析和驗(yàn)證該方法的可信度，我們將通過(guò)實(shí)際觀測(cè)資料與解算結(jié)果進(jìn)行對(duì)比。首先根據(jù)實(shí)際觀測(cè)的伸縮儀資料將2個(gè)分量潮汐線(xiàn)應(yīng)變觀測(cè)的方位角定為α1、α2，對(duì)地表潮汐實(shí)際的剪切應(yīng)變eθλ可用理論值進(jìn)行模擬響應(yīng)系數(shù)（潮汐因子）δk由2個(gè)不同方向的潮汐線(xiàn)應(yīng)變實(shí)測(cè)資料用調(diào)和分析擬合獲得，已由式（5）得出。將泰安臺(tái)伸縮儀觀測(cè)NS、EW向2個(gè)不同方位（方位角分別為 α1＝39.8000°、α2＝124.6167°）組成方程式（6）。

根據(jù)式（6）中兩 2個(gè)布設(shè)潮汐線(xiàn)應(yīng)變觀測(cè)值 e11、e22，（因作為已知，）即可求得另外 2個(gè)未知數(shù) eθθ、eλλ（劉序儼，1994；蔣駿等，1993；張雁濱等，1996）。

通過(guò)求得的 3個(gè)已知量 eθθ、eλλ及 eθλ進(jìn)而可以通過(guò)式（7）（8）計(jì)算得到其它的應(yīng)變矢量

求得最大主應(yīng)變e1、最小主應(yīng)變e2。最大剪應(yīng)變?chǔ)雍妥畲笾鲬?yīng)變方位角α依式（9）、（10）求得

通過(guò)MATLAB軟件編程，導(dǎo)入觀測(cè)數(shù)據(jù)，其中伸縮儀2分項(xiàng)數(shù)值提取為2013年全年日均值。通過(guò)運(yùn)算得出最大主應(yīng)變e1、最小主應(yīng)變e2、最大剪應(yīng)變?chǔ)雍妥畲笾鲬?yīng)變方位角α，并分別示于圖1～4。

通過(guò)對(duì)比可清晰看到，主應(yīng)力及剪應(yīng)力自年初開(kāi)始呈正弦波趨勢(shì)轉(zhuǎn)變，應(yīng)力開(kāi)始不斷積累，應(yīng)力方位角一致，到105天（3月15日）附近到達(dá)應(yīng)力積累頂點(diǎn)，后又不斷進(jìn)行釋放，一直到220天（8月10日）附近應(yīng)力釋放完全，最大主應(yīng)力恢復(fù)到年初水準(zhǔn)；而最小主應(yīng)力、最大剪應(yīng)力略高于年初水平，此時(shí)最大主應(yīng)力方位角突然變?yōu)榉聪?，接續(xù)了應(yīng)力的積累到釋放的過(guò)程。

值得注意的是，在第195天（7月14日）附近時(shí)，主應(yīng)力及剪應(yīng)力出現(xiàn)了4天左右的相對(duì)遲滯，隨后應(yīng)力加速釋放。從最大主應(yīng)變方位角時(shí)序圖中可以明顯看出，195～202天（7月14日～7月21日）附近的方位角出現(xiàn)了一次停滯，這應(yīng)該是某一個(gè)不明作用力對(duì)原年變趨勢(shì)的一個(gè)附加作用。從方位角及主應(yīng)力變化形態(tài)上判斷，該作用力為壓應(yīng)力。從第195天（7月14日）一直到2013年年底，整個(gè)洞體應(yīng)變擬合后的最大主應(yīng)變一直維持在呈相對(duì)固定斜率的加速積累狀態(tài)，積累量及變化速率均高于上半年；最小主應(yīng)變的變化趨勢(shì)則表現(xiàn)為相對(duì)遲滯，變化量及變化速率明顯低于上半年；而與之對(duì)應(yīng)的主應(yīng)變方位角則變化較為突出，這可能是新應(yīng)變作用力與原正常應(yīng)變之間相互作用的結(jié)果（牛安福等，2003）。

圖1 泰安臺(tái)洞體應(yīng)變擬合最大主應(yīng)變時(shí)序變化

圖2 泰安臺(tái)洞體應(yīng)變擬合最小主應(yīng)變時(shí)序變化

為了驗(yàn)證這一結(jié)論，我們比對(duì)了同樣反映地表應(yīng)變量變化的體應(yīng)變數(shù)據(jù)。首先提取體應(yīng)變2013年全年整時(shí)值，通過(guò)多項(xiàng)式曲線(xiàn)擬合對(duì)數(shù)據(jù)進(jìn)行處理。相關(guān)步驟為：將觀測(cè)序列劃分為p個(gè)數(shù)據(jù)段，每段用m次多項(xiàng)式擬合。約束條件為：相鄰數(shù)據(jù)段中擬合曲線(xiàn)在共同處的r階導(dǎo)數(shù)相等。設(shè)所要求的曲線(xiàn)方程為

通過(guò)MATLAB軟件實(shí)現(xiàn)公式的解算。將數(shù)據(jù)一次帶入方程，運(yùn)算得出體應(yīng)變擬合數(shù)據(jù)并繪成圖5。

圖3 泰安臺(tái)洞體應(yīng)變擬合最大剪應(yīng)變時(shí)序變化

圖4 泰安臺(tái)洞體應(yīng)變擬合最大主應(yīng)變方位角（弧度）時(shí)序變化

由圖5可知，體應(yīng)變對(duì)地表地應(yīng)力觀測(cè)值在170天（6月18日）左右因年變而開(kāi)始逐漸反向恢復(fù)，但在第195天（7月14日）附近時(shí)突然加速恢復(fù)，之后地應(yīng)力加速積累。該應(yīng)力在體應(yīng)變記錄上表現(xiàn)為壓應(yīng)力。截止到205天（7月25日）附近，通過(guò)計(jì)算發(fā)現(xiàn)，應(yīng)力變化量速率達(dá)到年變量速率的2倍。從250天（9月8日）起至2013年年底，擬合值變化速率放緩，明顯低于上半年同期變化量。上述結(jié)果清晰證明了體應(yīng)變應(yīng)力的加速積累與洞體應(yīng)變擬合解算結(jié)果是一致的。

綜上所述，在2013年的195天（7月14日）起，泰安臺(tái)所處區(qū)域出現(xiàn)一力源尚不明了的應(yīng)力變化，導(dǎo)致地表觀測(cè)應(yīng)變量高于正常值。該應(yīng)力為壓應(yīng)力，造成地表巖石擠壓，其日作用力應(yīng)變量為正常值的2倍。由于該應(yīng)變作用力一直與原背景應(yīng)力之間相互作用，放大了最大主應(yīng)變的變化量，對(duì)最小主應(yīng)變則起到抑制作用，同時(shí)也造成了最大主應(yīng)變方位角的不斷改變。

圖5 體應(yīng)變多項(xiàng)式曲線(xiàn)擬合時(shí)序變化

4 與重力觀測(cè)結(jié)果的對(duì)比

通過(guò)對(duì)洞體應(yīng)變2測(cè)向的擬合解算，并分析提取的相關(guān)參數(shù)后發(fā)現(xiàn)了自2013年第195天（7月14日）到年底出現(xiàn)了不明力源的應(yīng)力加速積累，導(dǎo)致了地表應(yīng)變呈加速受壓狀態(tài)。因地下物質(zhì)受壓會(huì)導(dǎo)致重力加速度的轉(zhuǎn)變，結(jié)合泰安臺(tái)PET重力儀監(jiān)測(cè)數(shù)據(jù)，我們同樣采取擬合法提取了重力參數(shù)，經(jīng)分析發(fā)現(xiàn)，該時(shí)間段內(nèi)重力加速度也發(fā)生相應(yīng)的突然轉(zhuǎn)變（見(jiàn)圖6中的第200天（7月19日）附近出現(xiàn)的變化）。隨后重力加速度出現(xiàn)了異常（如圖6所示）。重力數(shù)據(jù)的整個(gè)變化趨勢(shì)直接對(duì)應(yīng)了洞體應(yīng)變、體應(yīng)變擬合提取的解算結(jié)果，這也驗(yàn)證了這一不明力源的作用。

圖6 重力觀測(cè)K-L法最佳直線(xiàn)擬合時(shí)序變化

5 結(jié)論

對(duì)洞體應(yīng)變觀測(cè)數(shù)據(jù)進(jìn)行Venedikov調(diào)和分析計(jì)算，給出了洞體應(yīng)變固體潮觀測(cè)資料的維氏調(diào)和分析結(jié)果。計(jì)算結(jié)果表明應(yīng)變觀測(cè)NS向固體潮幅度及均方差明顯優(yōu)于EW向。泰安臺(tái)洞體應(yīng)變觀測(cè)響應(yīng)靈敏度、洞室條件等因素都符合前兆臺(tái)站儀器運(yùn)行管理規(guī)范所需要求。

通過(guò)對(duì)泰安臺(tái)2個(gè)測(cè)向潮汐線(xiàn)應(yīng)變組合觀測(cè)進(jìn)行了平面應(yīng)變狀態(tài)建模，提取了最大主應(yīng)變、最小主應(yīng)變、最大剪切應(yīng)變及最大主應(yīng)變方位角等信息，在對(duì)比研究中發(fā)現(xiàn)了該區(qū)域在2013年6月后出現(xiàn)了應(yīng)變場(chǎng)地應(yīng)力的異常，通過(guò)對(duì)比同樣反映地表應(yīng)變場(chǎng)變化的體應(yīng)變擬合值，發(fā)現(xiàn)該應(yīng)變場(chǎng)異常為客觀存在，并且表現(xiàn)突出。綜合數(shù)據(jù)判斷：2013年自195天（7月14日）起，泰安臺(tái)所處區(qū)域，出現(xiàn)一不明原因應(yīng)變力，造成了地表一定的應(yīng)變量高于正常值。主要表象為壓應(yīng)力，造成地表巖石擠壓。日作用力應(yīng)變量為正常應(yīng)變量的2倍。該應(yīng)變作用力一直與原正常形應(yīng)變之間相互作用，放大了最大主應(yīng)變的變化量，對(duì)最小主應(yīng)變則起到抑制作用，同時(shí)也造成了最大主應(yīng)變方位角的不斷改變。通過(guò)分析應(yīng)力場(chǎng)異常與重力場(chǎng)異常的對(duì)應(yīng)關(guān)系，也驗(yàn)證了這一不明作用場(chǎng)的變化。