引用本文
童茂松. 泥质砂岩激发极化弛豫时间谱的正则化反演[J]. 物探与化探, 2015,39(1): 186-191.
TONG Mao-Song. The regularization inversion of induced polarization relaxation time spectrum of agrillaceous sand[J]. Geophysical and Geochemical Exploration, 2015,39(1): 186-191.  
Doi:10.11720/wtyht.2015.1.30
Permissions
Copyright©2015, 《物探与化探》编辑部
《物探与化探》编辑部 所有
泥质砂岩激发极化弛豫时间谱的正则化反演
童茂松
大庆钻探工程公司 测井公司,黑龙江 大庆 163412

作者简介: 童茂松(1971-),男,高级工程师,2001年在吉林大学电子科学与技术学院获理学博士学位,2004年从该校博士后流动站和大庆油田博士后科研工作站出站,研究方向为测井方法研究与测井仪器研发。

收稿日期: 2014-08-12
基金: 中国石油天然气集团公司技术中心专项(05H219)
摘要

激发极化衰减谱是一系列单指数衰减的线性叠加,因此可以对衰减谱进行多指数反演,得到的弛豫时间谱可用于储层孔隙结构和渗透性评价,是一种非常有前途的技术。由于实际测井的信噪比相对较低,常规的奇异值分解方法难以获得合理的弛豫时间谱,为此,采用正则化约束,结合奇异值分解,对合成数据和实际测量的衰减谱进行反演,得到连续的弛豫时间谱;同时,讨论了弛豫时间谱的时间常数分布、正则化因子以及信噪比对反演结果的影响。研究结果表明:该方法适用于信噪比较低的实际测量结果的反演,弛豫时间常数的分布点数为32~64较为合适;随着正则化因子的增加,反演结果逐渐变得平滑,存在一个最优的正则化因子,用于获得最为合理的弛豫时间谱;在双对数图中,最优正则化因子随着信噪比的增加而线性降低,因此通过测量结果的信噪比,可以求取最优正则化因子,进而将该技术应用于实际。

关键词: 激发极化; 弛豫时间谱; 正则化反演; 泥质砂岩; 信噪比; 奇异值分解
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2015)01-0186-06
The regularization inversion of induced polarization relaxation time spectrum of agrillaceous sand
TONG Mao-Song
Well Logging Company, Daqing Drilling and Prospecting Engineering Company, Daqing 163412, China
Abstract

The time-domain induced polarization (IP) decay curve is a linear superposition of IP signal from individual pores. It can be transformed to relaxation time spectrum related to the pore structure and can be used to estimate the permeability of the argillaceous sand. So the IP logging is a useful tool for reservoir characterizing. The singular value decomposition (SVD) method can be used to extract the relaxation time spectra from the decay data with high signal-to-noise ratio (SNR). This method is unstable and invalid for the practical downhole use because of the noisy decay data. In this study, a method of regularization was applied to extracting continuous IP relaxation time spectra from simulated and real decay data. Influence of the regularization factor, SNR and relaxation time arrangement on the spectra was also discussed. Firstly, for the simulated decay data with noise, the spectra become smoother with increasing regularization factor. There is an optimal value which can be used to get the most reasonable spectrum. The optimal factor decreases linearly with increasing SNR on a log-log scale. This relationship can be used to predict the optimal factor from a known SNR. Secondly, the influence of relaxation time arrangement on the spectra shows that the suitable relaxation distribution point ranges from 32 to 64. Finally, the developed algorithm and the prediction of the optimal factor are very suitable for the inversion of the spectra of the real data.

Keyword: induced polarization; relaxation time spectrum; regularization method inversion; argillaceous sand; signal-to-noise ratio; singular value decomposition

激发极化法, 又称激电法, 广泛应用于金属、石油矿产资源以及地下水等方面的探测, 是目前最有效的地球物理勘探方法之一[1, 2]。根据研究对象的不同, 可以分为时间域激发极化法和频率域激发极化法, 前者是研究激发极化效应随时间变化的规律, 后者研究激发极化效应随频率变化规律, 就方法原理而言, 二者是等效的[3]。

何继善院士于1976年提出了“ 双频激电法” 理论, 其核心是同时供含有高、低两种频率的双频电流, 同时测量高、低两种频率的电位差而得到幅频率, 克服了谐波次数越高, 信号越小的缺点, 该方法得到了推广应用, 发现了一大批矿产, 并且在勘查地下水和解决工程问题中得到成功的应用[3, 4, 5, 6]。对于需要提取多个频率地电信息的频谱激电来说, 双频激电法仍然难以一次获得较为完整的激电频谱。在双频法基础上发展起来了an伪随机信号激电法, 进一步克服了变频法和奇次谐波法的缺点, 对双频法进行了提高和完善, 把频率域激电法推进到了一个新的阶段[7, 8, 9, 10, 11]。

傅良魁教授[12, 13, 14, 15, 16]、罗延钟教授[17, 18, 19, 20, 21]、杨进教授[22, 23, 24]、刘崧教授[1, 25, 26, 27, 28]等也对激发极化进行了大量的理论和实验研究, 分析激发极化法的影响因素, 反演激发极化频谱模型的参数, 评价激发极化在找矿、找水、探测油气藏等方面的应用。

在油气勘探开发中, 通过测井曲线进行渗透率解释是储层评价的重要内容之一。目前较先进的是核磁共振测井技术, 但核磁共振测井的探测深度浅, 测井成本高, 且信噪比低[29]。激发极化测井可用于储层的孔隙结构和渗透率评价[30, 31, 32, 33, 34], 与核磁共振测井相比, 激发极化测井具有探测深度深、信噪比高、成本低廉等优点, 是一种具有潜质的井下渗透率测量方法。该方法是在供电结束后, 测量电压衰减谱或者极化率, 从而评价储层的孔渗特性[30, 34, 35, 36]。对于含水的泥质砂岩而言, 激发极化衰减谱由具有不同时间常数的指数衰减谱叠加而成, 通过多指数反演, 可以得到弛豫时间谱, 其包含的信息较极化率更为丰富[37, 38], 可以求取岩石的毛管压力曲线和孔隙结构[30, 34, 38, 39], 通过适当的刻度, 激发极化弛豫时间谱测井可以得到地层的渗透率曲线[40]。

多指数反演方法对于获得合理的弛豫时间谱非常重要, 目前采用的奇异值分解方法, 对于信噪比高的实验室测量衰减谱是适用的[38, 39], 但是在实际测井中, 由于信噪比较低, 该方法并不适合。为此, 笔者研究了正则化与奇异值分解结合的方法, 对实际衰减谱进行反演, 得到了较好的效果。

1 岩石激发极化弛豫基础理论

激发极化理论建立在双电层基础之上, 即处于溶液中的岩石颗粒(尤其是黏土颗粒)因为带有多余的负离子而吸引溶液中的阳离子, 在其周围形成双电层(图1)。当在岩石的两端施加恒定外电场后, 则会使得双电层发生形变或局部浓度的变化, 即产生极化电场。当外电场去掉后, 由于扩散作用, 极化电场逐渐消失并且使离子恢复到原来的状态, 同时形成扩散电位, 即二次场衰减电位, 这就是激发极化电位。

图1 泥质砂岩激发极化原理示意

理论分析表明, 单个孔隙中的激发极化二次场电位的衰减满足单指数衰减规律[29, 30], 其时间常数为T, 满足条件T=R2/D[31], D为离子的扩散系数, R为孔径。因此, 每一个时间常数对应一个孔径。

实际的岩石内部是由一系列大小不等的孔隙群体组成的, 所以, 岩石的激发极化二次场电位V(t)是一系列单个孔隙极化电位的叠加, 二次场V(t)按照t=0时的二次场V(0)归一化:y(t)=V(t)/V(0), 然后进行多指数拟合, 即

yt=∑i=1N1Nfiexp⁡exp-tTi+ε(t)(1)

式中:ε (t)是噪声; N是孔隙的个数, 对于第i个孔隙, 其时间常数为Ti。泥质砂岩激发极化弛豫时间谱反演的目的是从式(1)中求解出系数fi。

将反演得到的f对弛豫时间常数T作图, 可以得到弛豫时间谱(图2)。对于时间常数为Ti的孔隙, 其半径为Ri= DTi, 不同的弛豫时间(Ti)对应着不同的孔隙半径(Ri)。图中纵坐标f代表不同弛豫时间常数的孔隙在总孔隙中所占的份额, 因此激发极化时间谱与泥质砂岩的孔隙结构相关[39]。

图2 典型的激发极化弛豫时间谱

2 激发极化多指数反演

式(1)可以看成非负约束最小二乘问题:

Φ≡‖Af-y‖2→min, Aij≡exp(-ti/Tj); fi≥0（2）

式中:y=(y1, y2, …, ym)T, 为激发极化衰减信号; f=(f1, f2, …, fn)T, 是时间常数为Tj 组分对应的幅度值; Tj(j=1, 2, …, n)为预先指定的弛豫时间系列, 典型的取法为在(Tmin, Tmax)区间对数均匀地选取n个点, 称之为弛豫时间分布点数, Tmin和Tmax为最短和最长的弛豫时间, 对应最小和最大的孔径。

式(2)可以采用奇异值分解方法求解, 但是结果受信噪比的影响非常大, 因此只适用于高信噪比的衰减谱反演。如果信噪比较低, 反演结果不稳定, 尤其是对于实际测井结果而言更不适合。为此, 对方程的解进行约束, 采用正则化反演方法, 求如下优化问题:

min:‖y-Af‖2+α2‖Wf‖2, fi≥0, 1≤i≤m（3）

式中, α 称为正则化因子, W是约束条件矩阵。

采用正则化与奇异值分解结合, 进行激发极化弛豫时间谱反演。为了直接应用奇异值分解对基于模型约束的正则化方程(3)进行求解, 可以将模型约束项用方程组形式表示, 并作为附加限制条件合并到原始数据方程组中, 采用最简单约束矩阵, 即单位阵I, 得到新的方程组:

min:‖y-Af‖2+α2‖If‖2, 　fi≥0, 1≤i≤m（4）

得到解:

f˙=Vdiag1σ1, 1σ2, …, 1σr, 0, …, 0·(UT·y)。（5）

反演的关键是非负约束的实现。这里采用Prammer提出的一种方法, 通过引入投影算子, 借助迭代消去负解[38]。具体做法是:

(1)计算试验解f;

(2)若所有的解分量fi≥ 0, 则接受x为线性方程组的解, 否则, 找出最负的解分量fi= min(fi)1≤i≤m, 删除矩阵A与fi对应的列, 形成新的矩阵A', 同时置fi=0, 并将其从原解向量中删除, 得到长度减少1的新向量f ';

(3)A⇔ A' , 转步骤(1)。通过上述过程能够得到原线性方程组的一组非负解。

3 结果与讨论
3.1 构造数据的反演

为了验证反演算法的有效性, 构造一个双峰特征的光滑弛豫时间谱(图2), 进而通过式(1)求取衰减谱, 在衰减谱中叠加不同的高斯噪声, 使得构造的数据具有不同的信噪比。信噪比定义为:

SNR=V(0)/σ, （6）

其中, σ 为噪声ε (t)的标准偏差。

采用不同正则化因子(α 为10-6, 10-4, 10-2, 100), 对无噪声(SNR=∞ )的衰减谱进行反演, 时间常数分布点数为64点, 分布范围为0.1~10 000 ms, 对数平均分布, 结果如图3所示。可以看出, 正则化因子α 对反演结果有影响, 当α 为10-6和10-4时, 反演结果与标准谱一致; 当该因子进一步增加时, 反演得到的弛豫时间谱与构造谱之间的差异逐渐明晰, 当α =100时, 反演的结果呈现单峰特征。

图3 α 不同时的无噪声数据反演结果

同时还研究了无噪声数据反演过程中, 时间常数T的分布对反演结果的影响。在0.1~10 000 ms范围内, 以对数平均的方式, 得到不同分布点数(8, 16, 32, 64, 128), α =10-4时的反演结果见图4。可以看出, 当分布点数为8和16时, 反演的结果与构造谱差异明显, 随着分布点数的增加, 反演的结果与构造谱逐渐接近, 当分布点数大于32, 反演的结果与构造谱接近。可以预见, 分布点数大于128时, 谱更为接近, 但是分布点数增加, 计算的时间也增加, 因此综合考虑反演效果和计算时间, 分布点数为32~64较为合适。

图4 分布点数对反演结果的影响

采用不同正则化因子(α 分别为0.01, 0.05, 0.2, 10), 对含噪声(SNR=100)的衰减谱进行反演, 时间常数分布点数为64点, 分布范围为0.1~10 000 ms, 对数平均分布, 结果见图5。图中可见, α 对含噪声数据的反演结果的影响比对无噪声数据更加明显, 随着α 的增加, 反演的弛豫时间谱逐渐变平滑; 当α 从0.01增加到10, 峰的个数由4个变化到1个, 当α =0.2, 相应的反演结果与标准谱一致。这些结果表明, 对于含噪声的衰减谱, α 决定了约束矩阵对反演结果的影响程度, 当α 过小时, 反演结果包含噪声信息, α 过大时, 则谱变得过于平滑, 与真实结果差异明显。因此, 存在一个最优正则化因子(α opt), 对于SNR=100的数据而言, α opt=0.2。

图5 α 对含噪声数据(SNR=100)反演结果的影响

同样采用不同正则化因子(α =0.02, 0.2, 2, 20)对SNR=20的数据进行反演, 所得结果(图6)说明正则化因子对此类数据反演结果的影响与图5类似, α opt=0.2。

图6 α 对含噪声数据(SNR=20)反演结果的影响

综合分析图5、图6, 可见对于含噪声的数据, 存在一个最优的正则化因子, 使得反演结果与构造的谱吻合, 而且最优正则化因子与数据的信噪比有关, 信噪比越高, 最优正则化因子越小。

为了将正则化反演方法应用于实际资料, 需要得到最优正则化因子α opt, 以便得到最为合理的弛豫时间谱。前面的分析和文献结果均表明, α opt取决于信噪比SNR, 进而影响反演谱的形态。噪声水平定义为1/SNR, α opt和SNR之间的关系可以写成[42]:α opt=c(1/SNR)b, 其中b和c为系数, 可以通过理论分析得到这两个系数。将噪声叠加到合成衰减谱上, 通过不同的正则化因子反演, 得到不同的弛豫时间谱, 通过将反演结果与标准谱之间的差异最小化, 可以得到对应于不同噪声水平下的最优化正则化因子, 从而可以求出b和c。经过上述过程, 可得到最优正则化因子和信噪比的关系(图7)。

图7 最优正则化因子和信噪比的关系

由图7可见, 在双对数图中, 最优正则化因子随信噪比的增加而线性降低。通过拟合, 可以得到二者的关系:

lgαopt=2.19-1.44lgSNR, r=0.994。（7）

相关系数达到0.994。为了验证该公式的正确性, 重新构造了标准谱, 信噪比SNR=40, 通过上式可求得α opt=0.75, 采用正则化因子分别为0.1α opt、α opt和10α opt进行反演, 结果见图8。可以看出, 从幅度和形态两个方面考虑, 采用最优正则化因子的反演结果与构造的标准谱最为接近, 其他两个明显有差异, 证实了式(6)的有效性。

图8 α 对含噪声数据(SNR=40)反演结果的影响

3.2 实际岩芯测量结果的反演

为了进一步评价正则化反演技术对于激发极化弛豫时间谱的有效性, 将该技术应用于实际岩芯测量结果的反演, 通过四电极法[38, 39], 测量了一块实际岩芯的衰减谱。

该岩芯来自大庆油田, 孔隙度15.6%, 渗透率为5.36× 10-3μ m2。数据采集由NI PCI-6070E 完成, 最大采样速率1.25 Ms/s, 供电和断电时间均为120 s。测量标准偏差为0.6 mV, V(0)约为30 mV, 由式(6)计算得到SNR=50, 由式(7)计算得到α opt=0.55, 采用正则化因子分别为0.1α opt、α opt和10α opt进行反演, 结果如图9所示, 同时图中还给出了通过压汞得到的岩芯孔隙结构。

图9 α 对实测数据(SNR=50)反演结果的影响

激发极化弛豫时间谱和压汞得到的孔隙结构具有对应性, 因此能够从二者的接近程度判断反演结果的有效性[39]。从图9可以看出, 采用α opt反演的结果与孔隙结构最为接近, 从而证明式(7)得到的最优正则化因子可用于实际岩芯测量结果(信噪比较低)的反演。

4 结论

正则化反演适合于含噪声的激发极化弛豫时间谱反演, 分布点数为32~64较为合适。

随着正则化因子的增加, 反演的弛豫时间谱逐渐变得平滑, 存在一个最优正则化因子, 在双对数曲线上, 该因子随信噪比的增加而线性增加。对于实际测量的数据, 采用最优正则化因子反演的弛豫时间谱与实际孔隙结构对应非常好。

The authors have declared that no competing interests exist.

参考文献
[1] 刘崧. 谱激电法[M]. 北京: 中国地质大学出版社, 1998. [本文引用:2]
[2] 孟建国. 浅谈激发极化法[J]. 吉林地质, 2013, 32(2): 82-83, 89. [本文引用:1]
[3] 何继善. 双频激电法[M]. 北京: 高等教育出版社, 2006. [本文引用:2]
[4] 何继善, 李大庆, 汤井田. 频谱激电非线性效应的理论模型[J]. 地球物理学, 1995, 38(5): 662-669. [本文引用:1]
[5] 何继善, 李大庆, 汤井田. 双频道频谱激电非线性效应研究[J]. 地球物理学报, 1995, 38(5): 670-675. [本文引用:1]
[6] 何继善, 熊彬, 鲍力知, 等. 直接消除电磁耦合的斩波去耦方法[J]. 地球物理学, 2006, 49(6): 1843-1850. [本文引用:1]
[7] 何继善. 频率域电法的新进展[J]. 地球物理学进展, 2007, 22(4): 1250-1254. [本文引用:1]
[8] 何继善. 伪随机三频电法研究[J]. 中国有色金属学报, 1994, 4(1): 1-7. [本文引用:1]
[9] 何继善. 电法勘探发展和展望[J]. 地球物理学报, 1997, 40(2) : 308-316. [本文引用:1]
[10] 何继善. 2n系列伪随机信号及应用[C]//中国地球物理学会年刊. 北京: 地震出版社, 1998: 199-199. [本文引用:1]
[11] 何继善, 柳建新. 伪随机多频相位法及其应用简介[J]. 中国有色金属学报, 2002(2): 374-376. [本文引用:1]
[12] 傅良魁. 激发极化法[M]. 北京: 地质出版社, 1982. [本文引用:1]
[13] 傅良魁. 复电阻率异常及空间分布规律[J]. 地质与勘探, 1981, 17(2): 46-53. [本文引用:1]
[14] 傅良魁, 张虎豹. 频谱激电法若干理论研究[J]. 物探与化探, 1985, 9(6): 410-420. [本文引用:1]
[15] 傅良魁, 李金铭, 陈兆洪. 埋藏极化体激电时间谱视参数的实验研究结果[J]. 地质与勘探, 1985, 21(12): 38-43. [本文引用:1]
[16] 傅良魁, 刘若谷. 时间域非线性激电效应的实验研究[J]. 桂林冶金地质学院学报, 1989, 9(1): 81-91. [本文引用:1]
[17] 罗延钟, 张桂青. 频率域激电法原理[M]. 北京: 地质出版社, 1988. [本文引用:1]
[18] 罗延钟, 方胜. 极化椭球体上频谱激电法的异常性态[J]. 桂林冶金地质学院学报, 1985, 1(5): 49-63. [本文引用:1]
[19] 罗延钟, 许妙立, 吴之训. 面极化球体上视复电阻率的频率特性和模拟相似性准则[J]. 桂林冶金地质学院学报, 1986, 6(2): 145-155. [本文引用:1]
[20] 罗延钟, 姚建阳, 何展翔, 等. 球形极化体上点源装置频谱激电异常的高级近似算法[J]. 物化探计算技术, 1987, 6(2): 113-122. [本文引用:1]
[21] 罗延钟, 张胜业, 熊彬. 天然源激电法的可行性[J]. 地球物理学报, 2003, 46(1): 125-130. [本文引用:1]
[22] 杨进, 孟海东, 傅良魁. 激电找水数据处理解释方法与系统[J]. 物探与化探, 1998, 22(3): 204-210. [本文引用:1]
[23] 杨进, 刘兆平. 天然场激发极化法在多金属矿区的野外试验效果[J]. 地学前缘, 2008, 15(4): 217-221. [本文引用:1]
[24] 杨进, 谭捍东, 傅良魁. 被动源激发极化法的野外试验结果[J]. 现代地质, 1998, 12(3): 437-441. [本文引用:1]
[25] 刘崧. 激电法直接探测油气藏的有效性及提高其应用效果的途径谱激电法, 地质科技情报, 1992, 11(4): 81-86. [本文引用:1]
[26] 刘崧. 计算视复电阻率的新的近似公式[J]. 地球物理学报, 1988, 31(6): 687-694. [本文引用:1]
[27] 刘崧, 官善友, 高鹏飞. 求极化椭球体真Cole-Cole 参数的联合谱激电反演[J]. 地球物理学报, 1994, 37(S1): 542-551. [本文引用:1]
[28] 刘崧, 徐建华. 用复电阻率法探测油气藏的研究[J]. 地球科学: 中国地质大学学报, 1997, 22(6): 627-632. [本文引用:1]
[29] 王宏建. 利用激发极化测井确定地层渗透率的研究进展[J]. 测井技术, 2008, 32(6): 509-513. [本文引用:2]
[30] Barker R D. A note on the induced polarization of the Bunter sand stone[J]. Geoexploration, 1975, 13: 227-233. [本文引用:4]
[31] Böner F D, Schopper J R, Weller A. Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements[J]. Geophysical Prospecting, 1996, 44: 583-602. [本文引用:2]
[32] Lesmes D P, Morgan F D. Dielectric spectroscopy of sedimentary rocks[J]. Journal of Geophysical Research, 2001: 3329-13346. [本文引用:1]
[33] Marshall D J, Madden T R. Induced polarization, a study of its causes[J]. Geophysics, 1959, 24: 790-816. [本文引用:1]
[34] Worthington P F, Collar F A. Relevance of induced polarization to quantitative formation evaluation[J]. Marine and Petroleum Geology, 1984(1): 14-26. [本文引用:3]
[35] Titov K, Komarov V, Tarasov V, et al. Theoretical and experimental study of time domain-induced polarization in water-saturated sand s[J]. Journal of Applied Geophysics, 2002, 50: 417-433. [本文引用:1]
[36] 王伟男, 余钦范, 姜亦忠. 激发极化电位衰减谱的连续测井方法[J]. 石油地球物理勘探, 2010, 45(增刊1): 214-216. [本文引用:1]
[37] Vinegar H J, Waxman M H. In-situ method for determining formation permeability[P]. U S Patent : 4743854, 1988. [本文引用:1]
[38] Tong M S, Wang W N, Li L, et al. Estimation of permeability of shaly sand reservoir from induced polarization relaxation time spectra[J]. Journal of Petroleum Science and Engineering, 2004, 45: 1-10. [本文引用:5]
[39] Tong M, Li L, Wang W, et al. Determining capillary pressure curve, pore size distribution and permeability from induced polarization of shaly sand [J]. Geophysics, 2006, 71(3): 33-40. [本文引用:5]
[40] Vinegar H J, Waxman M H. In-situ method for determining pore size distribution, capillary pressure and permeability[P]. U S Patent : 4644283, 1987. [本文引用:1]
[41] Prammer M G. NMR pore size distributions and permeability at the well site[J]. SPE, 1994: 28368. [本文引用:1]
[42] Gallegos D P, Smith D M. A NMR technique for the analysis of pore structure: determination of continuous pore size distributions. Journal of Colloid and Interface Science, 1988, 122(1): 143-153. [本文引用:1]