蘇淑娟　鄒春紅　王合乾　孫　豪　鄒鐘毅　吳月波

基于Matlab的井水位反演含水層體應(yīng)變可視化系統(tǒng)的研制1

蘇淑娟1，2）鄒春紅2）王合乾2）孫豪2）鄒鐘毅2）吳月波2）

1）山東省地震局，濟(jì)南 250000 2）煙臺(tái)地震監(jiān)測(cè)中心臺(tái)，山東煙臺(tái) 264000

基于反演含水層體應(yīng)變的數(shù)學(xué)模型，在Matlab GUI界面下研發(fā)1款高效、便捷的體應(yīng)變反演軟件，實(shí)現(xiàn)了規(guī)定時(shí)間閾內(nèi)連續(xù)體應(yīng)變值和實(shí)時(shí)體應(yīng)變曲線的一鍵獲得功能。軟件的用戶界面操作簡(jiǎn)便，經(jīng)實(shí)例驗(yàn)證，其運(yùn)算速度和結(jié)果均可滿足計(jì)算和分析需求。

井水位數(shù)據(jù)含水層體應(yīng)變反演軟件研制

引言

目前，利用井水位固體潮效應(yīng)和同震變化反演含水層參數(shù)的計(jì)算軟件多使用地震前兆信息系統(tǒng)EIS2000（蔣駿，2000），該系統(tǒng)可實(shí)現(xiàn)觀測(cè)井水位的預(yù)處理，后經(jīng)維涅第科夫（Venedikevo）調(diào)和分析求得井水位的潮汐因子。由于EIS2000缺少由潮汐因子反演含水層體應(yīng)變值的計(jì)算模塊，因此在反演體應(yīng)變方面，諸多學(xué)者（史浙明等，2012；楊柳等，2014；秦雙龍等，2014）依賴人工讀取、逐個(gè)統(tǒng)計(jì)井水位的差值，經(jīng)EIS2000求得潮汐因子，再利用井水位變量與體應(yīng)變關(guān)系式（張昭棟等，1999；劉序儼等，2009，2013）得出體應(yīng)變。上述計(jì)算方法僅顯示最大同震振幅時(shí)間節(jié)點(diǎn)或長(zhǎng)周期體應(yīng)變均值，無(wú)法獲取連續(xù)體應(yīng)變變化曲線，不能滿足井含水層體應(yīng)變實(shí)時(shí)、連續(xù)的研究需求。

基于上述研究現(xiàn)狀，本文整合孔隙彈性固體潮效應(yīng)和耦合效應(yīng)數(shù)學(xué)模型，利用Matlab強(qiáng)大的數(shù)學(xué)運(yùn)算能力，在Matlab GUI界面下研發(fā)1款高效、便捷的體應(yīng)變反演軟件，實(shí)現(xiàn)分鐘值與小時(shí)值時(shí)間閾內(nèi)體應(yīng)變值的一鍵獲得功能。承壓井水位本身可看作靈敏的含水層體應(yīng)變儀，若能實(shí)現(xiàn)井體應(yīng)變值的實(shí)時(shí)可視化導(dǎo)出，則可視為在此處布設(shè)了一個(gè)體應(yīng)變儀，進(jìn)而可密切關(guān)注地震構(gòu)造斷裂附近的體應(yīng)變變化情況，及時(shí)獲取地震前兆信息，為地震分析預(yù)報(bào)服務(wù)。

1　研究現(xiàn)狀及軟件開(kāi)發(fā)條件

1.1　研究現(xiàn)狀

1.1.1定量研究現(xiàn)狀

利用承壓井水位固體潮效應(yīng)和同震變化反演含水層體應(yīng)變，國(guó)內(nèi)外學(xué)者做了大量定量研究。Bredehoeft（1967）基于井-含水層系統(tǒng)水位對(duì)潮汐的響應(yīng)，計(jì)算孔隙度和儲(chǔ)水系數(shù)等含水層參數(shù)；Rhoads等（1979）通過(guò)理論公式推導(dǎo)建立了潮汐頻率和氣壓對(duì)水位的影響；Roeloffs（1996）的研究中考慮了井水位與含水層孔隙壓力的差別和滯后；尹京苑等（2000）經(jīng)擬合實(shí)際觀測(cè)井水位變化得到觀測(cè)井附近含水層內(nèi)的平均應(yīng)力變化率；車(chē)用太等（2006）通過(guò)井水位固體潮和推導(dǎo)解析式得到含水層參數(shù)；Narasimhan等（1984）、張昭棟等（1999）、劉序儼等（2009，2013）基于孔隙彈性介質(zhì)理論，得到井水位變化量與體應(yīng)變之間的關(guān)系式。19世紀(jì)，人們發(fā)現(xiàn)在潮汐和重荷載等因素影響下，井水水位會(huì)發(fā)生波動(dòng)（Wangle，2000）；20世紀(jì)初，固體和流體之間的相互影響的關(guān)系被分析研究并提出有效應(yīng)力的原理（Terzaghi，1923）；此后，小應(yīng)變情況下各向同性介質(zhì)的一般控制方程及三維各向異性多孔介質(zhì)固結(jié)方程被提出（Biot，1941，1956），Rice等（1976）在此基礎(chǔ)上推導(dǎo)出孔隙彈性耦合方程式，并提出不同數(shù)值模型研究水庫(kù)誘發(fā)地震（RIS）。

1.1.2實(shí)踐研究現(xiàn)狀

利用井水位固體潮效應(yīng)和同震變化反演含水層參數(shù)的實(shí)踐研究成果頗豐，如史浙明等（2012）利用同震水位階變資料和水位固體潮效應(yīng)，反演了含水層對(duì)汶川S8.0地震產(chǎn)生的體應(yīng)變響應(yīng)；楊柳等（2014）基于華北地區(qū)承壓井水位動(dòng)態(tài)，反演了華北地區(qū)含水層體應(yīng)變場(chǎng)的等值線變化圖；秦雙龍等（2014）利用福建井水位同震階變資料，反演了汶川S8.0地震和日本W(wǎng)9.0地震產(chǎn)生的體應(yīng)變變化；針對(duì)山東省的井水位，耿杰等（2008）、王學(xué)聚等（2013，2017）、蘇淑娟等（2016）曾對(duì)部分井孔同震效應(yīng)進(jìn)行不同程度的研究。

1.1.3反演軟件研究現(xiàn)狀

地震前兆信息系統(tǒng)EIS2000可實(shí)現(xiàn)原始井水位觀測(cè)數(shù)據(jù)的擬合、濾波、去趨勢(shì)和固體潮改正等，被廣泛用于反演含水層參數(shù)。EIS2000輔以BETCO程序（Toll等，2007）去氣壓效應(yīng)和去海潮效應(yīng)（李艷蕓等，2006；曹井泉等，2010），進(jìn)行Venedikevo調(diào)和分析求得井水位的潮汐因子。然而，EIS2000并未研發(fā)由潮汐因子反演含水層體應(yīng)變值的計(jì)算模塊，需要基于國(guó)內(nèi)外反演及正演數(shù)學(xué)模型，先統(tǒng)計(jì)水位差值，再利用井水位變化量與體應(yīng)變之間的關(guān)系式計(jì)算得出體應(yīng)變波峰值，且只能演示最大振幅時(shí)間節(jié)點(diǎn)或者幾個(gè)月的平均值，無(wú)法繪制體應(yīng)變曲線圖，故反演的體應(yīng)變值不能滿足Bodvarsson（1970）提出的“通達(dá)深部含水層的井可被視為天然的應(yīng)力計(jì)”的實(shí)時(shí)觀測(cè)與研究需求。

1.2　軟件開(kāi)發(fā)條件

1.2.1Matlab優(yōu)勢(shì)

Matlab在數(shù)值計(jì)算方面優(yōu)勢(shì)明顯，可進(jìn)行矩陣運(yùn)算、繪制函數(shù)和數(shù)據(jù)、實(shí)現(xiàn)算法、創(chuàng)建界面、連接其它程序等，主要應(yīng)用于工程計(jì)算、控制設(shè)計(jì)、信號(hào)處理與通訊、圖像處理、信號(hào)檢測(cè)、金融建模等領(lǐng)域。其多層面、多功能集成的平臺(tái)便于分析預(yù)報(bào)人員開(kāi)發(fā)所需的算法和程序。

1.2.2觀測(cè)井?dāng)?shù)據(jù)支撐

山東地震地下流體井網(wǎng)建設(shè)始于1979年，截至1984年11月，共有21口井取得國(guó)家地震局的驗(yàn)收證書(shū)。隨著“九五”時(shí)期井網(wǎng)的數(shù)字化改造，以聊古一井為代表的部分觀測(cè)井的部分測(cè)項(xiàng)于1998年實(shí)現(xiàn)數(shù)字化觀測(cè)，2006—2007年“十五”建設(shè)時(shí)期，井網(wǎng)觀測(cè)基本實(shí)現(xiàn)數(shù)字化。原始水位數(shù)據(jù)觀測(cè)周期長(zhǎng)、連續(xù)性好，為軟件研發(fā)提供了數(shù)據(jù)支持。

2　反演含水層體應(yīng)變數(shù)學(xué)模型的選取

2.1　固體潮效應(yīng)模型

對(duì)于封閉性較好的承壓含水層，可理想化地假設(shè)水位變化僅由含水層所受體應(yīng)變的變化引起。因此，根據(jù)孔隙彈性介質(zhì)理論可得到井水位的變化量與體應(yīng)變之間關(guān)系（張昭棟等，1999；劉序儼等，2009，2013）：

其中，表示井水位的潮汐因子，為含水層內(nèi)水的密度，g為地球表面重力加速度，為含水層孔隙度，E、E分別為巖石固體顆粒和孔隙流體的體積模量，d表示井水位的變化量，D表示含水層的體應(yīng)變量。對(duì)于封閉性較好的水位井，可通過(guò)井水位潮汐因子反演井孔含水層的體應(yīng)變。

2.2 　彈性耦合效應(yīng)模型

根據(jù)Biot（1941，1955，1956）提出的理論，在孔隙連通的彈性骨架中，由流體和固體骨架構(gòu)成具有守恒性質(zhì)的流-固耦合彈性系統(tǒng)。固體骨架具有壓縮性和剪切強(qiáng)度，流體壓縮后可多方向流動(dòng)。假設(shè)流-固耦合彈性系統(tǒng)中單位立方體的變形是完全可逆的，則流體與固體應(yīng)變分量分別遵守經(jīng)典的彈性理論。Biot（1955）認(rèn)為，各向同性巖石在應(yīng)力作用下具有對(duì)稱(chēng)特性，孔隙彈性耦合模型可采用如下幾何方程式：

其中，為應(yīng)變張量，、為位移張量。

對(duì)稱(chēng)二階張量（固體應(yīng)變張量）的分量在形式上應(yīng)與應(yīng)力張量保持一致，同時(shí)也體現(xiàn)了總角度變化等于2個(gè)角度相加的角應(yīng)變特征。應(yīng)力應(yīng)變的變化會(huì)引起流-固耦合彈性系統(tǒng)中固體骨架的孔隙被壓實(shí)或拉伸，即公式（2）反映的應(yīng)變張量引起位移張量的變化。水相對(duì)周?chē)鷰r體不可壓縮并流動(dòng)，由此引起孔隙壓的變化。若承壓井含水層較為封閉，深部孔隙壓的變化將傳至地表，表現(xiàn)為井水位的升降。由孔隙彈性耦合方程（公式（2））可以估算，在不排水條件下，1mm的水位變化可由含水層約4.9×10-10的體應(yīng)變引起（趙永紅等，2017）。因此，將實(shí)際觀測(cè)水位的變化量乘以4.9×10-10，可反演出井含水層流-固耦合系統(tǒng)的體應(yīng)變值。

3　軟件設(shè)計(jì)及功能實(shí)現(xiàn)

3.1　軟件設(shè)計(jì)

3.1.1操作界面

操作界面主要由groundwater.fig和groundwater.m文件完成。groundwater.fig文件負(fù)責(zé)軟件中各控件的布局、屬性的設(shè)定、工具欄及菜單欄的設(shè)計(jì)；groundwater.m文件由Matlab自動(dòng)生成，在其基礎(chǔ)上完成用戶所定義的所有函數(shù)功能。

3.1.2后臺(tái)處理程序

后臺(tái)處理程序可分為數(shù)據(jù)導(dǎo)入、分析類(lèi)型設(shè)定、分析計(jì)算、結(jié)果保存及繪圖，其執(zhí)行順序依次排列。數(shù)據(jù)導(dǎo)入與處理部分負(fù)責(zé)從外部導(dǎo)入數(shù)據(jù)，判斷數(shù)據(jù)類(lèi)型；設(shè)定部分根據(jù)用戶的研究需求，將分析類(lèi)型設(shè)定為分鐘值和小時(shí)值；分析計(jì)算部分在前2部分的基礎(chǔ)上，按式（1）、式（2）計(jì)算；數(shù)據(jù)保存及繪圖則用于保存分析結(jié)果和完成可視化繪圖。

3.1.3數(shù)據(jù)導(dǎo)入與處理

水位數(shù)據(jù)文件可通過(guò)地震行業(yè)內(nèi)網(wǎng)運(yùn)行的前兆處理系統(tǒng)、地震前兆信息系統(tǒng)（EIS2000）等軟件下載。軟件利用uigetfile函數(shù)識(shí)別所選擇的數(shù)據(jù)文件類(lèi)型后，將在后臺(tái)依據(jù)回調(diào)函數(shù)（callback_function）讀入并直接繪出水位埋深曲線圖，展示在繪圖區(qū)第一坐標(biāo)軸Axes1，將句柄文件保存在Axes1。數(shù)據(jù)導(dǎo)入功能分“導(dǎo)入水位分鐘值”和“導(dǎo)入水位小時(shí)值”按鈕?；卣{(diào)函數(shù)（callback_function）為：

［filename，pathname］ = igetfile{'*.mat；*.txt；*.xls；*.xlsx；*.xlsb'，'DataFiles￡¨*.mat；*.txt；*.xls；*.xlsx；*.xlsb）'；...*.*'， 'All Files（*.*）'}，...'Select the Data file'）；

file=fullfile（pathname，filename）；

data=load（file）；

y=data；

axes（handles.axes1）；

fpath=［pathname filename］；

plot（y）；

axis tight；

legend（'水位埋深曲線'）；

xlabel（'小時(shí)'）；

ylabel（'水位/m'）；

3.1.4潮汐因子的輸入

采用Edit Text設(shè)計(jì)潮汐因子的輸入。根據(jù)需求，在給定的類(lèi)型框下輸入潮汐因子（圖1），采用get_function函數(shù)實(shí)現(xiàn)數(shù)據(jù)傳遞，進(jìn)行后續(xù)計(jì)算。

input = str2num（get（hObject，'String'））；

handles.input=input；

guidata（hObject， handles）；

num=get（handles.edit1，'string'）；

num=str2num（num）

潮汐因子的正確取值是該軟件取得可靠結(jié)果的關(guān)鍵。借鑒前人基于大震同震響應(yīng)反演含水層體應(yīng)變的實(shí)踐經(jīng)驗(yàn)（史浙明等，2012；楊柳等，2014；秦雙龍等，2014），首先通過(guò)EIS2000對(duì)水位數(shù)據(jù)進(jìn)行擬合、濾波，再利用BETCO程序（Toll等，2007）去除水位的氣壓效應(yīng)，去除干擾因素后剩余的信息可用來(lái)計(jì)算潮汐因子值。采用Venedikov調(diào)和分析法分析大震前幾個(gè)月處理后的小時(shí)值水位數(shù)據(jù)計(jì)算潮汐因子。調(diào)和分析結(jié)果表明，固體潮分波中M2波的精度最高、振幅最大，可靠程度相對(duì)較高，故采用M2波潮汐因子。據(jù)觀測(cè)井距海岸線距離與體應(yīng)變響應(yīng)關(guān)系理論（曹井泉等，2010），對(duì)于距海岸線50km內(nèi)的水位井，在計(jì)算水位M2潮汐波振幅和相位的基礎(chǔ)上，還需進(jìn)行海潮負(fù)荷改正（李艷蕓等，2006）。

3.1.5反演體應(yīng)變值

輸入潮汐因子值后，選擇“固體潮效應(yīng)模型”或“耦合效應(yīng)模型”（圖1、圖2），將執(zhí)行不同的數(shù)學(xué)運(yùn)算方程，在后臺(tái)對(duì)讀取的數(shù)據(jù)進(jìn)行運(yùn)算，2種模型的體應(yīng)變值曲線圖分別繪制在第2和第3坐標(biāo)軸（Axes2、Axes3）。設(shè)計(jì)的初衷為2種計(jì)算模型得出的體應(yīng)變值可相互參照，客觀反應(yīng)體應(yīng)變值，避免單一計(jì)算方法得出的結(jié)果因缺乏參照而誤差偏大。固體潮效應(yīng)模型反演體應(yīng)變值按鈕的回調(diào)函數(shù)（callback_function）為：

m=num；

A=1/m；

與去年相比，十大石油公司名次沒(méi)有明顯變動(dòng)。中國(guó)三大石油公司排名變化不大。中國(guó)石油的綜合排名依然保持在第3 位，中國(guó)石化名列第20 位，中國(guó)海油名列第32 位。

B1=str2double（B1）；

B1=getappdata（handles.m_file_open，'string'）；

for i = 2 : 1440；

f=B1（1， i） - B1（1，（i-1））

C（i）=f；

g = C*1000；

h=A*g；

axes（handles.axes2）；

plot（h）

圖1 　汶川MS 8.0地震時(shí)商河魯09井“反演體應(yīng)變值”界面（分鐘值）

3.1.6體應(yīng)變值的導(dǎo)出

計(jì)算結(jié)果隨反演的結(jié)束而自動(dòng)產(chǎn)出groundwater.txt文本文件，打開(kāi)后可瀏覽計(jì)算結(jié)果，也可用于數(shù)據(jù)分析。儲(chǔ)存函數(shù)的回調(diào)函數(shù)（callback_function）為：

handles.y=y；

guidata（hObject， handles）；

save groundwater.txt y -ascii；

3.1.7圖片的導(dǎo)出

Axes坐標(biāo)軸產(chǎn)出的圖片直觀易懂。本文開(kāi)發(fā)的軟件設(shè)計(jì)了“FIG.1”、“FIG.2”和“FIG.3”3個(gè)按鈕，分別代表Axes1、Axes2和Axes3產(chǎn)出的圖件，使用時(shí)可根據(jù)需求下載圖片并自由選擇下載路徑。圖片以figure type類(lèi)型保存在指定路徑，下載函數(shù)的回調(diào)函數(shù)（callback_ function）為：

axes（handles.axes1）；

if isempty（handles.axes1）；

return；

end

newFig = figure；

set（newFig，'Visible'，'off'）；

newAxes = copyobj（handles.axes1，newFig）； set（newAxes，'Units'，'default'，'Position'，'default'）；［filename，pathname］ = uiputfile（{'*.jpg'，'figure type（*.jpg）'}， '保存水位原始曲線'）； %axes2設(shè)為“保存固體潮效應(yīng)模型所得曲線”，axes3設(shè)為“保存耦合效應(yīng)模型所得曲線”

if isequal（filename，0）||isequal（pathname，0）

return；

else

fpath=fullfile（pathname，filename）；

end

f = getframe（newFig）；

f = frame2im（f）；

imwrite（f，fpath）；

圖2 　汶川MS 8.0地震時(shí)棲霞魯07井“反演體應(yīng)變值”界面（小時(shí)值）

Matlab的圖像文件（Figure file）自帶編輯器（圖3），點(diǎn)擊鼠標(biāo)即可讀取最大振幅值，進(jìn)行圖片縮放、打印、色彩更改、文本標(biāo)記、儲(chǔ)存、新建、去網(wǎng)格、修改坐標(biāo)軸等操作，簡(jiǎn)單便捷。系統(tǒng)運(yùn)行界面的底色設(shè)為白色，同時(shí)將3個(gè)坐標(biāo)軸的底色設(shè)為白色，產(chǎn)出圖片后可直接截屏保存3張豎向圖片。將原始水位數(shù)據(jù)、潮汐效應(yīng)模型反演結(jié)果和耦合效應(yīng)反演結(jié)果置于1張圖中對(duì)比分析，截屏后可直接使用。

圖3 　Figure圖片的顯示及編輯界面

3.2　分析功能的實(shí)現(xiàn)

程序編寫(xiě)完成后，需將其打包為可在Windows操作系統(tǒng)中獨(dú)立運(yùn)行的程序，既可封裝編程函數(shù)，使其免遭破壞，又利于推廣應(yīng)用。在Matlab中實(shí)現(xiàn)groundwater.m函數(shù)文件的封裝打包，組建后形成獨(dú)立的“GroundWater.exe”程序。程序運(yùn)行的主界面（圖4）包括水位數(shù)據(jù)選擇導(dǎo)入?yún)^(qū)、潮汐因子輸入?yún)^(qū)、反演體應(yīng)變值模型選擇區(qū)、繪圖展示區(qū)、圖片導(dǎo)出區(qū)及程序退出區(qū)。

圖4 “GroundWater.exe”程序運(yùn)行界面

4 　實(shí)例驗(yàn)證

以2008年5月12日汶川S8.0地震、2011年3月11日日本本州東海岸W9.0地震和2015年4月25日尼泊爾S8.1地震為時(shí)間節(jié)點(diǎn)，選取山東省9口（聊古1井、魯02井、魯03井、魯04井、魯07井、魯09井、魯15井、魯27井、魯32井）封閉性好、固體潮效應(yīng)清晰、同震響應(yīng)明顯、觀測(cè)時(shí)間較長(zhǎng)、水位資料完整的數(shù)字化井，進(jìn)行實(shí)例驗(yàn)證和分析。

4.1 　2種計(jì)算模型所得體應(yīng)變值的相互驗(yàn)證

以商河魯09井和棲霞魯07井為例，運(yùn)行界面分別如圖1、圖2、圖5和圖6所示。通過(guò)分析可知，“固體潮效應(yīng)模型”和“耦合效應(yīng)模型”2種計(jì)算方式在同一地震、同一觀測(cè)井中反演的體應(yīng)變的數(shù)值、量級(jí)以及應(yīng)力變化曲線基本一致，計(jì)算過(guò)程實(shí)現(xiàn)了水位和體應(yīng)變的協(xié)同分析，使得由水位反演的體應(yīng)變實(shí)時(shí)、連續(xù)、可視化，且井周?chē)鷰r體所受的應(yīng)力可形象地呈現(xiàn)。分鐘值模塊可為解釋水震波及含水層應(yīng)力變化提供實(shí)時(shí)觀測(cè)描述，水位埋深曲線中水震波的震蕩、體應(yīng)變的劇烈變化、水位震蕩的振幅和體應(yīng)變量均可直接讀?。▓D1、圖5和圖6）。小時(shí)值模塊可獲得體應(yīng)變隨朔望月的月相盈虧而產(chǎn)生的周期性變化（圖2），便于獲取宏觀水位與體應(yīng)變的周期性更迭特征，捕捉臨震前兆信息。

4.2　軟件反演結(jié)果與實(shí)際體應(yīng)變觀測(cè)值對(duì)比

對(duì)基于井水位的同震響應(yīng)數(shù)據(jù)反演的含水層體應(yīng)變值與實(shí)際體應(yīng)變觀測(cè)值進(jìn)行對(duì)比驗(yàn)證。地應(yīng)力的增加產(chǎn)生壓縮區(qū)和張性區(qū)，壓縮區(qū)水位將逐漸抬升，張性區(qū)水位下降（車(chē)用太等，2000）。由圖5的水位數(shù)據(jù)可知，尼泊爾S8.1地震發(fā)生時(shí)，棲霞井巖體骨架先壓縮致使水位迅速抬升，后拉張致使水位不斷下降，而后慢慢恢復(fù)。軟件反演得到的張性區(qū)及應(yīng)變?yōu)樨?fù)值的反演量級(jí)、數(shù)值和應(yīng)力方向均與實(shí)際觀測(cè)值一致（圖6、圖7），但擠壓區(qū)的反演結(jié)果小于實(shí)際觀測(cè)值，可能與棲霞井受擠壓時(shí)裂隙井所處位置、深度和固體骨架孔隙連通性有關(guān)。

圖5　 尼泊爾MS 8.1地震商河魯09井水位反演的體應(yīng)變曲線

圖6 　尼泊爾MS 8.1地震棲霞魯07井水位反演的體應(yīng)變曲線

圖7 　尼泊爾MS 8.1地震煙臺(tái)地震監(jiān)測(cè)中心臺(tái)的體應(yīng)變實(shí)際觀測(cè)曲線

5 　結(jié)論與討論

本文在Matlab GUI界面下實(shí)現(xiàn)了groundwater.m和groundwater.fig文件的函數(shù)編輯和操作窗口設(shè)計(jì)，由Command Window和Deployment tool將groundwater.m文件打包封裝為可在Windows操作系統(tǒng)中獨(dú)立運(yùn)行的“GroundWater.exe”程序。經(jīng)實(shí)例驗(yàn)證，該程序的用戶界面操作方便，通過(guò)2種模型反演的體應(yīng)變曲線均較真實(shí)地體現(xiàn)了井-含水層系統(tǒng)的時(shí)空變化規(guī)律，運(yùn)算速度和反演精度均能滿足體應(yīng)變的研究需求。

由該可視化系統(tǒng)反演的體應(yīng)變值反映了含水層中孔隙壓隨時(shí)間的變化情況，通過(guò)反演一口井、得出1套體應(yīng)變分鐘值和小時(shí)值曲線，相當(dāng)于在此處布設(shè)1臺(tái)有分鐘值和小時(shí)值采集功能的體應(yīng)變儀。此外，一般臺(tái)站的體應(yīng)變儀測(cè)得的數(shù)據(jù)僅反映地表附近的體應(yīng)變，并不能準(zhǔn)確反映處于地層深處的應(yīng)力狀態(tài)，而該系統(tǒng)根據(jù)與地下深處相連通的地下流體固體潮效應(yīng)和同震響應(yīng)，反演得到的體應(yīng)變則反映了地層一定深處的應(yīng)力狀態(tài)，對(duì)于理解地下流體的同震變化機(jī)理、了解地震的孕育過(guò)程和預(yù)測(cè)地震具有一定意義。同時(shí)，對(duì)于缺少鉆孔應(yīng)變資料的地區(qū)，反演的數(shù)據(jù)可作為有效補(bǔ)充。

致謝：感謝審稿專(zhuān)家提出寶貴修改意見(jiàn)。

曹井泉，朝倫巴根，劉耀煒，2010. 承壓井水位固體潮M2波海潮負(fù)荷改正. 地震研究，33（1）：75—80.

車(chē)用太，劉五洲，魚(yú)金子等，2000. 板內(nèi)強(qiáng)震的中地殼硬夾層孕震與流體促震假設(shè). 地震學(xué)報(bào)，22（1）：93—101.

車(chē)用太，魚(yú)金子，2006. 地震地下流體學(xué). 北京：氣象出版社，420—424.

耿杰，陳安方，潘雙進(jìn)，2008. 山東地下水動(dòng)態(tài)觀測(cè)井對(duì)2007年印尼8.5級(jí)地震的響應(yīng)特征. 西北地震學(xué)報(bào)，30（2）：173—178.

蔣駿，2000. 地震前兆信息處理與軟件系統(tǒng). 北京：地震出版社.

李艷蕓，李紹武，2006. 風(fēng)暴潮預(yù)報(bào)模式在渤海海域中的應(yīng)用研究. 海洋技術(shù)，25（1）：101—106.

劉序儼，鄭小菁，王林等，2009. 承壓井水位觀測(cè)系統(tǒng)對(duì)體應(yīng)變的響應(yīng)機(jī)制分析. 地球物理學(xué)報(bào)，52（12）：3147—3157.

劉序儼，鄭小菁，陳瑩等，2013. 承壓井與非承壓井水位潮汐效應(yīng)及其定量分析. 大地測(cè)量與地球動(dòng)力學(xué)，33（1）：35—39.

秦雙龍，廖麗霞，陳瑩等，2014. 利用福建井水位同震變化反演井-含水層體應(yīng)變及其意義探討. 內(nèi)陸地震，28（4）：353—359.

史浙明，王廣才，劉春國(guó)，2012. 基于汶川地震同震地下水位變化反演含水層體應(yīng)變. 地震學(xué)報(bào)，34（2）：215—223.

蘇淑娟，孫豪，鄒春紅等，2017. 魯07井水位與水溫同震響應(yīng)特征淺析. 齊魯?shù)卣鹂茖W(xué)專(zhuān)輯（2016合輯），3：97—105.

王學(xué)聚，殷海濤，王鵬，2013. 山東地下流體數(shù)字化井網(wǎng)對(duì)汶川8.0級(jí)地震的響應(yīng)分析. 地震地磁觀測(cè)與研究，34（1—2）：225—231.

王學(xué)聚，殷海濤，王慶林，2017. 山東地下流體數(shù)字化井網(wǎng)對(duì)特大地震的響應(yīng)分析. 國(guó)際地震動(dòng)態(tài)，（10）：32—39.

楊柳，馬建英，曹井泉等，2014. 利用華北地區(qū)承壓井水位資料反演含水層體應(yīng)變. 中國(guó)地震，30（2）： 249—259.

尹京苑，趙利飛，2000. 保山井水位異常的數(shù)值模擬. 西北地震學(xué)報(bào)，22（4）：397—401.

張昭棟，劉慶國(guó)，耿杰，1999. 由承壓井水位動(dòng)態(tài)反演水井含水層的應(yīng)力變化. 華南地震，19（1）：37—42.

趙永紅，謝雨晴，王航等，2017. 地震預(yù)測(cè)方法：地下流體方法. 地球物理學(xué)進(jìn)展，32（4）：1539—1547. Biot M. A., 1941. General theory of three-dimensional consolidation. Journal of Applied Physics, 12(2): 155—164.

Biot M. A., 1955. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 26(2): 182—185.

Biot M. A., 1956. General solutions of the equations of elasticity and consolidation for a porous material. Journal of Applied Mechanics, 78: 91—96.

Bodvarsson G., 1970. Confined fluids as strain meters. Journal of Geophysical Research, 75(14): 2711—2718.

Bredehoeft J. D., 1967. Response of well-aquifer systems to earth tides. Journal of Geophysical Research, 72(12): 3075—3087.

Narasimhan T. N., Kanehiro B. Y., Witherspoon P. A., 1984. Interpretation of earth tide response of three deep, confined aquifers. Journal of Geophysical Research: Solid Earth, 89(B3): 1913—1924.

Rhoads Jr. G. H., Robinson E. S., 1979. Determination of aquifer parameters from well tides. Journal of Geophysical Research: Solid Earth, 84(B11): 6071—6082.

Rice J. R., Cleary M. P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Reviews of Geophysics, 14(2): 227—241.

Roeloffs E., 1996. Poroelastic techniques in the study of earthquake-related hydrologic phenomena. Advances in Geophysics, 37: 135—195.

Terzaghi K., 1923. Die berechnung der durchlassigheitsziffer des tones aus dem verlauf der hydrodynamischen spanningserscheinungen. Sber Akad Wiss Wien, 132: 105—124.

Toll N. J., Rasmussen T. C., 2007. Removal of barometric pressure effects and Earth tides from observed water levels. Groundwater, 45(1): 101—105.

Wangle H. F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton: Princeton University Press, 4-10.

The Development of Visualization System for Inversion of Aquifer Volumetric Strain by the Data of the Well Water Level Through Matlab

Su Shujuan1,2), Zou Chunhong2), Wang Heqian2), Sun Hao2), Zou Zhongyi2)and Wu Yuebo2)

1) Shandong Earthquake Agency, Jinan 250000, China 2) Yantai Earthquake Monitoring Center Station, Yantai 264000, Shandong, China

Based on the mathematical models for the inversion of aquifer volumetric strain, we developed the system to obtain the value and the curve of continuous volumetric strain within the schedule time on the interface of Matlab GUI. After verification, we find that the user interface developed in this paper is easy to operate, and the operation speed and results can meet the needs of calculation and analysis.

Data of the well water level; Aquifer volume strain; Inversion; Software development

10.11899/zzfy20190214

山東省地震局重點(diǎn)研發(fā)項(xiàng)目（YF1703）

2018-11-01

蘇淑娟，女，生于1979年。工程師。研究方向：地震地下流體。E-mail：shujuan_su@163.com

蘇淑娟，鄒春紅，王合乾，孫豪，鄒鐘毅，吳月波，2019．基于Matlab的井水位反演含水層體應(yīng)變可視化系統(tǒng)的研制．震災(zāi)防御技術(shù)，14（2）：411—422．