基于Matlab的井水位反演含水层体应变可视化系统的研制1

苏淑娟　邹春红　王合乾　孙　豪　邹钟毅　吴月波

基于Matlab的井水位反演含水层体应变可视化系统的研制1

苏淑娟1，2）邹春红2）王合乾2）孙豪2）邹钟毅2）吴月波2）

1）山东省地震局，济南 250000 2）烟台地震监测中心台，山东烟台 264000

基于反演含水层体应变的数学模型，在Matlab GUI界面下研发1款高效、便捷的体应变反演软件，实现了规定时间阈内连续体应变值和实时体应变曲线的一键获得功能。软件的用户界面操作简便，经实例验证，其运算速度和结果均可满足计算和分析需求。

井水位数据含水层体应变反演软件研制

引言

目前，利用井水位固体潮效应和同震变化反演含水层参数的计算软件多使用地震前兆信息系统EIS2000（蒋骏，2000），该系统可实现观测井水位的预处理，后经维涅第科夫（Venedikevo）调和分析求得井水位的潮汐因子。由于EIS2000缺少由潮汐因子反演含水层体应变值的计算模块，因此在反演体应变方面，诸多学者（史浙明等，2012；杨柳等，2014；秦双龙等，2014）依赖人工读取、逐个统计井水位的差值，经EIS2000求得潮汐因子，再利用井水位变量与体应变关系式（张昭栋等，1999；刘序俨等，2009，2013）得出体应变。上述计算方法仅显示最大同震振幅时间节点或长周期体应变均值，无法获取连续体应变变化曲线，不能满足井含水层体应变实时、连续的研究需求。

基于上述研究现状，本文整合孔隙弹性固体潮效应和耦合效应数学模型，利用Matlab强大的数学运算能力，在Matlab GUI界面下研发1款高效、便捷的体应变反演软件，实现分钟值与小时值时间阈内体应变值的一键获得功能。承压井水位本身可看作灵敏的含水层体应变仪，若能实现井体应变值的实时可视化导出，则可视为在此处布设了一个体应变仪，进而可密切关注地震构造断裂附近的体应变变化情况，及时获取地震前兆信息，为地震分析预报服务。

1　研究现状及软件开发条件

1.1　研究现状

1.1.1定量研究现状

利用承压井水位固体潮效应和同震变化反演含水层体应变，国内外学者做了大量定量研究。Bredehoeft（1967）基于井-含水层系统水位对潮汐的响应，计算孔隙度和储水系数等含水层参数；Rhoads等（1979）通过理论公式推导建立了潮汐频率和气压对水位的影响；Roeloffs（1996）的研究中考虑了井水位与含水层孔隙压力的差别和滞后；尹京苑等（2000）经拟合实际观测井水位变化得到观测井附近含水层内的平均应力变化率；车用太等（2006）通过井水位固体潮和推导解析式得到含水层参数；Narasimhan等（1984）、张昭栋等（1999）、刘序俨等（2009，2013）基于孔隙弹性介质理论，得到井水位变化量与体应变之间的关系式。19世纪，人们发现在潮汐和重荷载等因素影响下，井水水位会发生波动（Wangle，2000）；20世纪初，固体和流体之间的相互影响的关系被分析研究并提出有效应力的原理（Terzaghi，1923）；此后，小应变情况下各向同性介质的一般控制方程及三维各向异性多孔介质固结方程被提出（Biot，1941，1956），Rice等（1976）在此基础上推导出孔隙弹性耦合方程式，并提出不同数值模型研究水库诱发地震（RIS）。

1.1.2实践研究现状

利用井水位固体潮效应和同震变化反演含水层参数的实践研究成果颇丰，如史浙明等（2012）利用同震水位阶变资料和水位固体潮效应，反演了含水层对汶川S8.0地震产生的体应变响应；杨柳等（2014）基于华北地区承压井水位动态，反演了华北地区含水层体应变场的等值线变化图；秦双龙等（2014）利用福建井水位同震阶变资料，反演了汶川S8.0地震和日本W9.0地震产生的体应变变化；针对山东省的井水位，耿杰等（2008）、王学聚等（2013，2017）、苏淑娟等（2016）曾对部分井孔同震效应进行不同程度的研究。

1.1.3反演软件研究现状

地震前兆信息系统EIS2000可实现原始井水位观测数据的拟合、滤波、去趋势和固体潮改正等，被广泛用于反演含水层参数。EIS2000辅以BETCO程序（Toll等，2007）去气压效应和去海潮效应（李艳芸等，2006；曹井泉等，2010），进行Venedikevo调和分析求得井水位的潮汐因子。然而，EIS2000并未研发由潮汐因子反演含水层体应变值的计算模块，需要基于国内外反演及正演数学模型，先统计水位差值，再利用井水位变化量与体应变之间的关系式计算得出体应变波峰值，且只能演示最大振幅时间节点或者几个月的平均值，无法绘制体应变曲线图，故反演的体应变值不能满足Bodvarsson（1970）提出的“通达深部含水层的井可被视为天然的应力计”的实时观测与研究需求。

1.2　软件开发条件

1.2.1Matlab优势

Matlab在数值计算方面优势明显，可进行矩阵运算、绘制函数和数据、实现算法、创建界面、连接其它程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模等领域。其多层面、多功能集成的平台便于分析预报人员开发所需的算法和程序。

1.2.2观测井数据支撑

山东地震地下流体井网建设始于1979年，截至1984年11月，共有21口井取得国家地震局的验收证书。随着“九五”时期井网的数字化改造，以聊古一井为代表的部分观测井的部分测项于1998年实现数字化观测，2006—2007年“十五”建设时期，井网观测基本实现数字化。原始水位数据观测周期长、连续性好，为软件研发提供了数据支持。

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后台处理程序

后台处理程序可分为数据导入、分析类型设定、分析计算、结果保存及绘图，其执行顺序依次排列。数据导入与处理部分负责从外部导入数据，判断数据类型；设定部分根据用户的研究需求，将分析类型设定为分钟值和小时值；分析计算部分在前2部分的基础上，按式（1）、式（2）计算；数据保存及绘图则用于保存分析结果和完成可视化绘图。

3.1.3数据导入与处理

水位数据文件可通过地震行业内网运行的前兆处理系统、地震前兆信息系统（EIS2000）等软件下载。软件利用uigetfile函数识别所选择的数据文件类型后，将在后台依据回调函数（callback_function）读入并直接绘出水位埋深曲线图，展示在绘图区第一坐标轴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函数实现数据传递，进行后续计算。

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），将执行不同的数学运算方程，在后台对读取的数据进行运算，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坐标轴产出的图片直观易懂。本文开发的软件设计了“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）包括水位数据选择导入区、潮汐因子输入区、反演体应变值模型选择区、绘图展示区、图片导出区及程序退出区。

图4 “GroundWater.exe”程序运行界面

4 　实例验证

以2008年5月12日汶川S8.0地震、2011年3月11日日本本州东海岸W9.0地震和2015年4月25日尼泊尔S8.1地震为时间节点，选取山东省9口（聊古1井、鲁02井、鲁03井、鲁04井、鲁07井、鲁09井、鲁15井、鲁27井、鲁32井）封闭性好、固体潮效应清晰、同震响应明显、观测时间较长、水位资料完整的数字化井，进行实例验证和分析。

4.1 　2种计算模型所得体应变值的相互验证

以商河鲁09井和栖霞鲁07井为例，运行界面分别如图1、图2、图5和图6所示。通过分析可知，“固体潮效应模型”和“耦合效应模型”2种计算方式在同一地震、同一观测井中反演的体应变的数值、量级以及应力变化曲线基本一致，计算过程实现了水位和体应变的协同分析，使得由水位反演的体应变实时、连续、可视化，且井周围岩体所受的应力可形象地呈现。分钟值模块可为解释水震波及含水层应力变化提供实时观测描述，水位埋深曲线中水震波的震荡、体应变的剧烈变化、水位震荡的振幅和体应变量均可直接读取（图1、图5和图6）。小时值模块可获得体应变随朔望月的月相盈亏而产生的周期性变化（图2），便于获取宏观水位与体应变的周期性更迭特征，捕捉临震前兆信息。

4.2　软件反演结果与实际体应变观测值对比

对基于井水位的同震响应数据反演的含水层体应变值与实际体应变观测值进行对比验证。地应力的增加产生压缩区和张性区，压缩区水位将逐渐抬升，张性区水位下降（车用太等，2000）。由图5的水位数据可知，尼泊尔S8.1地震发生时，栖霞井岩体骨架先压缩致使水位迅速抬升，后拉张致使水位不断下降，而后慢慢恢复。软件反演得到的张性区及应变为负值的反演量级、数值和应力方向均与实际观测值一致（图6、图7），但挤压区的反演结果小于实际观测值，可能与栖霞井受挤压时裂隙井所处位置、深度和固体骨架孔隙连通性有关。

图5　 尼泊尔MS 8.1地震商河鲁09井水位反演的体应变曲线

图6 　尼泊尔MS 8.1地震栖霞鲁07井水位反演的体应变曲线

图7 　尼泊尔MS 8.1地震烟台地震监测中心台的体应变实际观测曲线

5 　结论与讨论

本文在Matlab GUI界面下实现了groundwater.m和groundwater.fig文件的函数编辑和操作窗口设计，由Command Window和Deployment tool将groundwater.m文件打包封装为可在Windows操作系统中独立运行的“GroundWater.exe”程序。经实例验证，该程序的用户界面操作方便，通过2种模型反演的体应变曲线均较真实地体现了井-含水层系统的时空变化规律，运算速度和反演精度均能满足体应变的研究需求。

由该可视化系统反演的体应变值反映了含水层中孔隙压随时间的变化情况，通过反演一口井、得出1套体应变分钟值和小时值曲线，相当于在此处布设1台有分钟值和小时值采集功能的体应变仪。此外，一般台站的体应变仪测得的数据仅反映地表附近的体应变，并不能准确反映处于地层深处的应力状态，而该系统根据与地下深处相连通的地下流体固体潮效应和同震响应，反演得到的体应变则反映了地层一定深处的应力状态，对于理解地下流体的同震变化机理、了解地震的孕育过程和预测地震具有一定意义。同时，对于缺少钻孔应变资料的地区，反演的数据可作为有效补充。

致谢：感谢审稿专家提出宝贵修改意见。

曹井泉，朝伦巴根，刘耀炜，2010. 承压井水位固体潮M2波海潮负荷改正. 地震研究，33（1）：75—80.

车用太，刘五洲，鱼金子等，2000. 板内强震的中地壳硬夹层孕震与流体促震假设. 地震学报，22（1）：93—101.

车用太，鱼金子，2006. 地震地下流体学. 北京：气象出版社，420—424.

耿杰，陈安方，潘双进，2008. 山东地下水动态观测井对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. 山东地下流体数字化井网对特大地震的响应分析. 国际地震动态，（10）：32—39.

杨柳，马建英，曹井泉等，2014. 利用华北地区承压井水位资料反演含水层体应变. 中国地震，30（2）： 249—259.

尹京苑，赵利飞，2000. 保山井水位异常的数值模拟. 西北地震学报，22（4）：397—401.

张昭栋，刘庆国，耿杰，1999. 由承压井水位动态反演水井含水层的应力变化. 华南地震，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

山东省地震局重点研发项目（YF1703）

2018-11-01

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

苏淑娟，邹春红，王合乾，孙豪，邹钟毅，吴月波，2019．基于Matlab的井水位反演含水层体应变可视化系统的研制．震灾防御技术，14（2）：411—422．