基于二维NLCG反演的水文地质结构辨识研究

黄磊，侯泽明，韩萱，刘志强，苟青松

(1.内蒙古农业大学水利与土木建筑工程学院，内蒙古 呼和浩特 010018；2.内蒙古自治区水资源保护与利用重点实验室，内蒙古 呼和浩特 010018；3.重庆地质矿产研究院，重庆 401120)

在物探技术不断发展的背景下，各国学者在进行水文地质勘察时采用了多种物探手段。赵国泽等(2004)综合AMT和TEM对煤矿断层进行了富水性探测；而Matthew等(2005)采用了电阻率法、电磁法解译分析了含水层沉积结构；Alastair F等(2011)对地下水流场的运动路径进行测定时，主要运用了电阻率法、地质雷达和地震法。地球物理方法在含水层结构识别中的综合应用已经成为近年来的发展趋势(麻昌英等，2015；李承泽等，2018；刘瑞平等，2019)，但研究的难点和热点始终是含水层结构的精准反演解译(李学兰等，2016)。笔者通过大地电磁法(MT)识别理论含水层模型，再采用二维NLCG反演来解释电学性质与地层之间的对应关系，并尝试建立目标区域的地电模型来反映含水层结构，为核磁共振法获得含水层岩性、弛豫时间、含水率、渗透系数等相关水文地质参数奠定基础，从而形成一套水文地质结构综合辨识方法。其辨识成果可以有效弥补研究区水文地质资料的不足，也可为同类数据反演提供方法借鉴。

1 非线性共轭梯度反演方法

非线性共轭梯度法(Nonlinear Conjugate Gradient，NLCG)在二维大地电磁反演中得到了广泛的应用(张昆等，2011；赵维俊等，2014)。该方法过程简单，收缩速度快，内存占用小，具有高稳定性和高准确性(R.Dehghani et al.，2019)，可表示为：

d=F(m)+e

(1)

其中，d表示数据向量，m表示模型向量，e表示残差向量，F表示正演模拟函数。

WL Rodi等(2001)采用“正则化解”的方式得到了目标函数极小值。定义为：

(2)

笔者选用Polak-Ribiere的非线性共轭梯度算法，在线性方向搜索极小化模型序列：

m0=m

(3)

ψ(mq+αqpq)=minψ(mq+αpq)

(4)

mq+1=mq+αqpq，q=0,1,2……

(5)

然后，沿pq计算搜索步长αq来进行目标函数迭代达到极小化，搜索方向为：

p0=-C0g0

(6)

pq=-Cqgq+βqPq-1，q=1,2……

(7)

对高斯-牛顿线性搜索方法改进后为：

φ(α)=ψ(mq+αpq)

(8)

(9)

则最小值为：

(10)

2 含水层理论模型反演试验

笔者建立了2种不同产状的低阻含水层模型，用来验证二维NLCG反演算法对含水层反演是否准确有效。

图1是层状含水层模型。上覆30 m的覆盖层，电阻率为30 Ω·m；中间为厚30 m，10 Ω·m的低阻含水层；下部隔水层厚140 m，电阻率为120 Ω·m。通过有限元进行TE模式正演(图2)。由正演结果可以得出：电阻率分为3层，顶底两层为高阻层，中间层电阻率较低。尽管正演结果会受到低阻体积效应的影响，不同层的电阻值在分界处会渐进过渡，理论模型和实际结果在垂向分辨上存在一定误差，但可以分辨出不同阻值的分层情况，尤其对低阻含水层有一定的反应(李培熙等，2019；常威等，2019)。总体来说，正演模型可以较好地反映实际情况。

图1 层状含水层模型图Fig.1 Layered aquifer model

从二维NLCG反演结果(图2)可以清晰地分辨出大约30 m处为一、二层的分界线，但第二层底边界下延了100 m的深度。大体上地层模型和电性结果的相关性较好，特别是中间含水层能够被识别，说明二维NLCG算法对层状低阻层反演有一定的效果。

图2 层状含水层模型正反演图Fig.2 Forward and inversion of layered aquifer model

图3是局部含水层理论模型，上覆厚100 m、1 000 Ω·m的高阻层。局部为埋深200 m、长500 m、宽100 m的10 Ω·m低阻含水层；局部含水层所处的围岩电阻率为100 Ω·m，总深度600 m。

图3 局部含水层模型图Fig.3 Local aquifer model

图4为通过TE模式进行正演后得到局部含水层模型。图4中表层为高阻，中间夹低阻含水层，其体积效应的影响导致电阻值向四周不断升高。对正演模型进行二维NLCG反演解译，其反演结果见图4。高阻体下边界反映清晰，以200 Ω·m等值线为分界线。边界以下是低电阻围岩，其电性特征与理论模型对应良好。埋深在200～300 m的中心位置为低电阻层，其边界向外延扩。由上述分析可知二维NLCG反演对局部低阻体有较好的识别效果。

图4 局部含水层模型正反演图Fig.4 Forward and inversion of local aquifer model

3 含水层结构辨识试验

试验数据来源于内蒙古锡林河流域，数据采集选用EH-4电导率仪和地面核磁共振NUMIS-POLY系统(虎维岳等，2017；黄磊等，2018)。根据测区地质资料和电性特征，将视电阻区间划分为三级。Ⅰ级：电阻率5～25 Ω·m，岩性为湖积砂-砾石；Ⅱ级：电阻率25～120 Ω·m，岩性主要为粗砂-砾岩；Ⅲ级：电阻率大于125 Ω·m，岩性常见于凝灰岩、玄武岩。笔者通过对比实测数据的解译结果与钻孔资料，以验证二维NLCG反演解译的准确性。

图5为电阻率剖面图，该测线均匀布设9个测点，剖面深度约为180 m。浅表层为高阻层，厚度不大，产状平缓，推断其为凝灰岩薄层。30 m以下为低阻层，推断岩性为砂土，并赋存孔隙水，该层可能是富水性较好的含水层。

图5 电阻率剖面图Fig.5 Resistivity profile

在大地电磁二维NLCG反演基础上，对剖面中点进行核磁共振探测，其测试结果见图6，解译数据见表1。分析得出：5 m处为表层土壤水。18 m以下为较厚含水层，解译推测含水率为1.20%～5.53%，渗透系数为1.51×10-7～4.93×10-5m/s，岩性为砂砾石层，核磁测深未见隔水底板。

图6 SNMR解译结果图Fig.6 Interpretation result graph of SNMR

表1 SNMR解译数据表Tab.1 Interpretation data sheet for SNMR

通过对实测数据资料的综合分析，将垂直剖面进行了分层(图7)。垂直方向的电阻率在浅层显示出较大的变化趋势，而在深层显示出较小的变化趋势。从凝灰岩高阻层到中粗砂低阻层电阻率变化明显。根据钻孔资料，含水层主要为厚层松散湖积砂砾石与裂隙发育的砾岩组成，反演的视电阻率较低。泥岩和黏性土为隔水层，电阻率较高。

图7 钻孔资料与反演解译结果对比图Fig.7 Comparison chart of borehole data and inversion interpretation results

4 结论

(1)笔者设计了2种常见的含水层模型，并基于非线性共轭梯度法实现了大地电磁二维NLCG反演。通过模型试验和实测分析验证了二维NLCG反演算法适用于低阻含水层地质体的电磁探测解译。借此使得含水层地电结构的精准辨识工作可以顺利开展，同时为水文物探数据解释提供了理论依据和方法借鉴。

(2)基于二维NLCG反演理论，解译了实测数据，得出测区电阻率与岩性特征能较好的对应，并且钻孔资料也佐证了含水介质的岩性与赋存规律。通过SNMR反演解译与钻孔对比，发现含水层在深度上发生了一定程度的错位，产生这种情况可能是由于测区高压线的电磁干扰所致，而将大地电磁法和核磁共振相互参照解译可以有效对含水层的地质结构进行辨识。