基于Matlab的井水位反演含水層體應變可視化系統的研制1

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

基于Matlab的井水位反演含水層體應變可視化系統的研制1

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

1）山東省地震局，濟南 250000 2）煙臺地震監(jiān)測中心臺，山東煙臺 264000

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

井水位數據含水層體應變反演軟件研制

引言

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

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

1　研究現狀及軟件開發(fā)條件

1.1　研究現狀

1.1.1定量研究現狀

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

1.1.2實踐研究現狀

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

1.1.3反演軟件研究現狀

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

1.2　軟件開發(fā)條件

1.2.1Matlab優(yōu)勢

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

1.2.2觀測井數據支撐

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

2　反演含水層體應變數學模型的選取

2.1　固體潮效應模型

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

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

2.2 　彈性耦合效應模型

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

其中，為應變張量，、為位移張量。

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

3　軟件設計及功能實現

3.1　軟件設計

3.1.1操作界面

操作界面主要由groundwater.fig和groundwater.m文件完成。groundwater.fig文件負責軟件中各控件的布局、屬性的設定、工具欄及菜單欄的設計；groundwater.m文件由Matlab自動生成，在其基礎上完成用戶所定義的所有函數功能。

3.1.2后臺處理程序

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

3.1.3數據導入與處理

水位數據文件可通過地震行業(yè)內網運行的前兆處理系統、地震前兆信息系統（EIS2000）等軟件下載。軟件利用uigetfile函數識別所選擇的數據文件類型后，將在后臺依據回調函數（callback_function）讀入并直接繪出水位埋深曲線圖，展示在繪圖區(qū)第一坐標軸Axes1，將句柄文件保存在Axes1。數據導入功能分“導入水位分鐘值”和“導入水位小時值”按鈕?；卣{函數（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（'小時'）；

ylabel（'水位/m'）；

3.1.4潮汐因子的輸入

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

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

handles.input=input；

guidata（hObject， handles）；

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

num=str2num（num）

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

3.1.5反演體應變值

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

m=num；

A=1/m；

與去年相比，十大石油公司名次沒有明顯變動。中國三大石油公司排名變化不大。中國石油的綜合排名依然保持在第3 位，中國石化名列第20 位，中國海油名列第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地震時商河魯09井“反演體應變值”界面（分鐘值）

3.1.6體應變值的導出

計算結果隨反演的結束而自動產出groundwater.txt文本文件，打開后可瀏覽計算結果，也可用于數據分析。儲存函數的回調函數（callback_function）為：

handles.y=y；

guidata（hObject， handles）；

save groundwater.txt y -ascii；

3.1.7圖片的導出

Axes坐標軸產出的圖片直觀易懂。本文開發(fā)的軟件設計了“FIG.1”、“FIG.2”和“FIG.3”3個按鈕，分別代表Axes1、Axes2和Axes3產出的圖件，使用時可根據需求下載圖片并自由選擇下載路徑。圖片以figure type類型保存在指定路徑，下載函數的回調函數（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設為“保存固體潮效應模型所得曲線”，axes3設為“保存耦合效應模型所得曲線”

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地震時棲霞魯07井“反演體應變值”界面（小時值）

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

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

3.2　分析功能的實現

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

圖4 “GroundWater.exe”程序運行界面

4 　實例驗證

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

4.1 　2種計算模型所得體應變值的相互驗證

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

4.2　軟件反演結果與實際體應變觀測值對比

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

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

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

圖7 　尼泊爾MS 8.1地震煙臺地震監(jiān)測中心臺的體應變實際觀測曲線

5 　結論與討論

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

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

致謝：感謝審稿專家提出寶貴修改意見。

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

車用太，劉五洲，魚金子等，2000. 板內強震的中地殼硬夾層孕震與流體促震假設. 地震學報，22（1）：93—101.

車用太，魚金子，2006. 地震地下流體學. 北京：氣象出版社，420—424.

耿杰，陳安方，潘雙進，2008. 山東地下水動態(tài)觀測井對2007年印尼8.5級地震的響應特征. 西北地震學報，30（2）：173—178.

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

李艷蕓，李紹武，2006. 風暴潮預報模式在渤海海域中的應用研究. 海洋技術，25（1）：101—106.

劉序儼，鄭小菁，王林等，2009. 承壓井水位觀測系統對體應變的響應機制分析. 地球物理學報，52（12）：3147—3157.

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

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

史浙明，王廣才，劉春國，2012. 基于汶川地震同震地下水位變化反演含水層體應變. 地震學報，34（2）：215—223.

蘇淑娟，孫豪，鄒春紅等，2017. 魯07井水位與水溫同震響應特征淺析. 齊魯地震科學專輯（2016合輯），3：97—105.

王學聚，殷海濤，王鵬，2013. 山東地下流體數字化井網對汶川8.0級地震的響應分析. 地震地磁觀測與研究，34（1—2）：225—231.

王學聚，殷海濤，王慶林，2017. 山東地下流體數字化井網對特大地震的響應分析. 國際地震動態(tài)，（10）：32—39.

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

尹京苑，趙利飛，2000. 保山井水位異常的數值模擬. 西北地震學報，22（4）：397—401.

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

趙永紅，謝雨晴，王航等，2017. 地震預測方法：地下流體方法. 地球物理學進展，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

山東省地震局重點研發(fā)項目（YF1703）

2018-11-01

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

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