2017, Vol. 32
2. 吉林大学地球科学学院，长春 130061
3. 核工业二四三大队，内蒙古赤峰 024000
2. College of Earth Sciences, Jilin University, Changchun 130026, China
3. Geologic Party No. 243, CNNC, Chifeng, Inner Mongolia Chifeng 024000, China
频谱激电法 (SIP) 也称复电阻率法 (CR)，是以岩、矿石的电化学性质差异为基础，广泛应用于环境和工程、矿产和土地污染勘查等方面.该方法的特点在于用于解释的参数丰富，Pelton等 (1978)年发表著名论文“多频激电的电磁耦合去除与含矿异常识别”，将Cole-Cole模型引入到地球物理领域中，采用Cole-Cole模型及其组合来拟合激发极化与电磁耦合，从而奠定了频谱激电法的基础，并且他认为频率相关系数c和时间常数τ可以用来区分矿与非矿；罗延钟和方胜 (1986)提出以近似公式为正演模型，用一个测点的视谱求取真谱参数的理论.刘崧等 (1994)认为由于等值性的存在，如果没有附加信息，仅反演一条视激电谱，不可能求得有限极化体的除频率相关系数之外的其他真Cole-Cole参数，于是提出了联合反演两个或多个测点的视谱求真谱参数的方法，并实现了椭球型极化体真谱参数的反演计算，但其并未考虑电磁耦合效应的影响.Ghorbani等 (2009)在考虑电磁耦合的前提下，利用视复电阻率的实部和虚部直接反演Cole-Cole模型参数的算法，他先用低频观测的视复电阻率振幅数据进行一维直流电测深反演，得到电阻率模型作为复电阻率反演中零频电阻率的初始模型，然后再同时反演四个参数；于立波 (2010)在忽略电磁耦合效应的基础上，固定零频电阻率不变，用相位对其他三个参数进行了阻尼最小二乘的反演；周峰 (2015)在忽略电磁耦合的情况下，对Ghorbani等 (2009)的算法进行了改进，他将直流反演得到的电阻率模型作为复电阻率反演中零频电阻率的初始模型后，不再更新零频电阻率这一参数，并用相位反演极化率m、时间常数τ和频率相关系数c；范翠松等 (2012)提出电场振幅和相位联合反演的方案，并同时反演四个Cole-Cole参数.
由于视复电阻率实部和虚部对不同Cole-Cole参数的敏感程度不同，而实部的数量级一般比虚部大很多，这往往导致当拟合差已经达到预先设定的阀值时，虚部拟合的程度依然很低，这时对虚部影响更大的时间常数τ和频率相关系数c往往不能得到理想的结果.针对这个问题，本文在考虑电磁耦合的基础上，提出了利用视复电阻率实、虚部加权反演的方案，并采用先反演零频电阻率ρ0和极化率m再反演频率相关系数c和时间常数τ的方法.本文首先推导了反演算法，接着通过几种理论模型试算，证明了加权反演算法的稳定性和优越性.
1 一维频谱激电实、虚部加权反演 1.1 阻尼最小二乘法反演理论在计算复电率反演问题时，设在某一个排列上，有N个观测数据di，i=1，2，…，N，定义mi为某一Cole-Cole参数，i=1，2，…，M，M为地下模型层数，定义F(m) 为正演数据，Wd为观测数据的协方差矩阵，其结构为
|
(1) |
根据 (尚通晓，2008)，定义目标函数Φ为
|
(2) |
定义模型约束条件为
|
(3) |
式中，α为阻尼因子.
我们将目标函数修改为
|
(4) |
对正演数据F(m) 泰勒展开得：
|
(5) |
式中，mk是第k次迭代的初始模型，J是雅克比矩阵，其元素为
将 (5) 式带入 (4) 式，得到线性化的目标函数为
|
(6) |
从上式可以看出Φ是 (Δm) 的函数，对目标函数Φ取极值条件，即令：
|
(7) |
化简得到阻尼最小二乘线性反演方程组为
|
(8) |
给定初始模型mk，由上式可求出模型修改量Δm，同时得到mk+1=mk+Δm作为新的初始模型进行下一次迭代，直到拟合差小于预先给定的阀值为止，此时的mk+1即为最终的反演结果.
1.2 一维频谱激电实、虚部加权反演前面的推到过程是在复数域下进行的，根据范翠松等 (2012)的理论，这种方法多解性较强，他采取电场振幅和相位联合反演的方法，取得了一定的效果.然而电场振幅和相位对不同Cole-Cole参数的灵敏度不同，电场振幅和相位等权情况下同时反演四个Cole-Cole参数得到的频率相关系数和时间常数的反演结果有时不是特别理想.本文在范翠松等 (2012)联合反演的基础上，提出利用视复电阻率实部和虚部加权反演的策略，采取先反演零频电阻率ρ0和极化率m再反演频率相关系数c和时间常数τ的方法，并在这两个过程中取不同的加权系数，建立目标函数为
|
(9) |
式中，ΦRe和ΦIm分别代表复电阻率实部和虚部的目标函数，λ1和λ2为加权系数.其中，ΦRe和ΦIm具体为
|
(10) |
式中，dRe和dIm分别为观测数据的实部和虚部，FRe(m) 和FIm(m) 分别为正演数据的实部和虚部，WdRe和WdRe分别为观测数据实部和虚部的协方差矩阵.
根据上一章的理论，在目标函数中加入模型约束Φm=α(Δm)T(Δm)，将 (9) 式改写为
|
(11) |
参考前文阻尼最小二乘反演的推到过程：将 (11) 式线性化近似，并令
|
(12) |
式中，JRe和JIm分别为数据实部和虚部的雅克比矩阵，将通过差分法求取.
定义每次迭代反演的拟合差为
|
(13) |
本节将通过对多个理论模型的反演试算，来检验加权反演算法的效果和优越性.在以下两个算例中，反演数据均为视电阻率的实部和虚部，选用偶极-偶极的测量装置.排列数为1个，每个排列道数为30道，发射和接收偶极长度均为100 m，选取7个观测频率，分别为：0.0001 Hz，0.001 Hz，0.01 Hz，0.1 Hz，1 Hz，5 Hz，10 Hz.
反演分别两步：第一步，固定时间常数和频率相关系数的值不变，反演零频电阻率和极化率；第二布，将第一步的反演结果作为初始模型，并固定零频电阻率和极化率的值不变，反演频率相关系数和时间常数.
图 1~图 2为频率相关系数和时间常数初始值不同时，零频电阻率ρ0和极化率m的反演结果.反演时取加权系数：λ1=100000、λ2=1，固定频率相关系数c和时间常数τ不变，只反演零频电阻率ρ0和极化率m.
|
图 1 二层极化模型不同初始模型零频电阻率反演结果 Figure 1 The inversion results of the resistivity of two layers polarization model with different initial models |
|
图 2 二层极化模型不同初始模型极化率反演结果 Figure 2 The inversion results of the chargeability of two layers polarization model with different initial models |
从表 1可知在初始模型1初始模型4中零频电阻率ρ0和极化率m的初始值相同，而时间常数τ和频率相关系数c的初始值不同，在这种情况下固定时间常数τ和频率相关系数c的值不变，反演零频电阻率ρ0和极化率m，令加权系数λ1=100000，λ2=1，反演结果如图 1~图 2所示.从图中可以看出，当频率相关系数c和时间常数τ的初始值各不相同时，零频电阻率ρ0和m极化率可得到较好的反演结果，说明这种分两步反演的方法是可行的.
|
表 1 二层极化模型 Table 1 Polarization model of two layers |
在上述反演算例都是在各参数初始值与真值差别在一定范围内不大的条件下进行的，当模型初始值与真值相差过大时，反演经常陷入局部极小，这是由于阻尼最小二乘法自身的缺陷造成的.根据以上分析，我们按照分两步反演的方法进行下面的反演，即先反演零频电阻率ρ0和极化率m，再反演频率相关系数c和时间常数τ.
首先，固定频率相关系数c和时间常数τ的值不变，对实、虚部目标函数进行四种不同加权，权系数分别为：λ1=100000和λ2=1、λ1=1和λ2=100、λ1=1和λ2=100000，反演零频电阻率ρ0和极化率m，其反演结果如图 3~图 4所示.
|
图 3 二层极化模型初始模型1不同加权系数下零频电阻率反演结果 Figure 3 The inversion results of the resistivity of two layers polarization model with different weighting coefficients in the first initial model |
|
图 4 二层极化模型初始模型1不同加权系数下极化率反演结果 Figure 4 The inversion results of the chargeability of two layers polarization model with different weighting coefficients in the first initial model |
由如图 3~图 4可以看出，当加权系数λ1=100000，λ2=1时，零频电阻率ρ0和极化率m的反演结果与真值最为接近，故将此时的反演结果作为零频电阻率ρ0和极化率m各自的初始值，并固定不变，反演频率相关系数c和时间常数τ.反演时依次取四种不同加权系数：λ1=100000和λ2=1、λ1=1和λ2=1、λ1=1和λ2=100、λ1=1和λ2=100000，频率相关系数c和时间常数τ的反演结果如图 5~图 6所示.
|
图 5 二层极化模型初始模型1不同加权系数下频率相关系数反演结果 Figure 5 The inversion results of the frequency dependence of two layers polarization model with different weighting coefficients in the first initial model |
|
图 6 二层极化模型初始模型1不同加权系数下时间常数反演结果 Figure 6 The inversion results of the time constant of two layers polarization model with different weighting coefficients in the first initial model |
当第一步反演权系数λ1=10000，λ2=1、第二步反演权系数λ1=1，λ2=100000时，反演得到f=0.0001 Hz的观测数据和正演数据最终拟合情况如图 7~图 8所示.
|
图 7 频率f=0.0001 Hz视复电阻率实部反演拟合曲线 Figure 7 The iterative fitting curve of real parts of complex resistivity when f=0.0001 Hz |
|
图 8 频率f=0.0001 Hz视复电阻率虚部反演拟合曲线 Figure 8 The iterative fitting curve of imaginary parts of complex resistivity when f=0.0001 Hz |
图 3~图 6中，横坐标均为深度，纵坐标依次为零频电阻率ρ0、极化率m、频率相关系数c和时间常数τ.图 7~图 8中，横坐标为发射偶极AB长度的一半，纵坐标分别为视复电阻率的实部和虚部.
|
|
表 2 三层极化模型 Table 2 Polarization model of three layers |
同上，取四种不同加权系数：λ1=100000和λ2=1、λ1=1和λ2=1、λ1=1和λ2=100、λ1=1和λ2=100000，固定反演频率相关系数c和时间常数τ的值不变，反演零频电阻率ρ0和极化率m，反演结果如图 9~图 10所示.
|
图 9 三层极化模型不同加权系数下零频电阻率反演结果 Figure 9 The inversion results of the resistivity of three layers polarization model with different weighting coefficients in the first initial model |
|
图 10 三层极化模型不同加权系数下极化率反演结果 Figure 10 The inversion results of the resistivity of three layers polarization model with different weighting coefficients in the first initial model |
从如图 9~图 10中我们可以看出，与两层极化模型时一样，当加权系数λ1=100000，λ2=1时，零频电阻率ρ0和极化率m的反演结果与真值最为接近，故将此时的反演结果分别作为零频电阻率ρ0和极化率m的初始值，并固定不变，然后反演频率相关系数c和时间常数τ，反演时取四种不同加权系数：λ1=100000和λ2=1、λ1=1和λ2=1、λ1=1和λ2=100、λ1=1和λ2=100000，反演结果如图 11~图 12所示.
|
图 11 三层极化模型不同加权系数下频率相关系数反演结果 Figure 11 The inversion results of the frequency dependence of three layers polarization model with different weighting coefficients in the first initial model |
|
图 12 三层极化模型不同加权系数下时间常数反演结果 Figure 12 The inversion results of the time constant of three layers polarization model with different weighting coefficients in the first initial model |
当第一步反演权λ1=100000，λ2=1、第二步反演权λ1=1，λ2=100000时，反演得到f=0.0001 Hz时的观测数据和正演数据最终拟合情况如图 13~图 14所示.图 9~图 12中，与二层模型的情况相同，横坐标均为深度，纵坐标为反演参数的值.图 13~图 14中，横坐标为发射偶极AB长度的一半，纵坐标分别为视复电阻率的实部和虚部.由图 7~图 8和图 13~图 14可以看出，两层和三层模型反演观测数据与正演数据的实部和虚部拟合程度较好，误差较小.
|
图 13 频率f=0.0001 Hz视复电阻率实部反演拟合曲线 Figure 13 The iterative fitting curve of real parts of complex resistivity when f=0.0001 Hz |
|
图 14 频率f=0.0001 Hz视复电阻率虚部反演拟合曲线 Figure 14 The iterative fitting curve of imaginary parts of complex resistivity when f=0.0001 Hz |
从图 3~图 4和图 9~图 10中可以看出，不同加权系数下零频电阻率ρ0和极化率m的反演结果是不同的.对于二层和三层模型，当λ1≥λ2时，由于视复电阻率实部的数量级大于虚部，所以此时实部在反演中起到主要作用，加权系数λ1=100000，λ2=1使零频电阻率ρ0和极化率m的反演结果都要稍好于加权系数λ1=1、λ2=1(等权) 的情况，但差别不大；当λ1 < λ2时，虚部在反演中起到的作用增大，零频电阻率ρ0和极化率m的反演结果与真值相差较大，无法反映真实模型的趋势，极化率m出现负值.
以上分析说明视复电阻率的实部对零频电阻率ρ0和极化率m较为敏感，反演零频电阻率ρ0和极化率m时应令实部的加权系数更大，以得到更好的结果.
从图 5~图 6和图 11~图 12中可以看出，不同加权系数下频率相关系数c和时间常数τ反演的结果是不同的.当λ1≥λ2时，实部在反演中起到主要作用，此时频率相关系数c和时间常数τ的反演结果比加权系数λ1 < λ2时的反演结果要差.对于二层模型，当加权系数λ1=100000，λ2=1时，如图 5a和图 6a所示，频率相关系数c的反演结果比其他三种加权情况下的反演结果都要差 (与真值相差最大)，而时间常数τ的反演没有得到能正确反映模型原始趋势的结果，在深度为300 m处出现假极值，当权λ1=λ2=1时，如图 6b所示，时间常数τ的反演结果趋势不明显；对于三层模型，当权λ1=100000，λ2=1时，如图 11a和图 12a所示，频率相关系数c的反演结果可以反映理论模型的趋势，但相对于其他三种加权情况误差最大，时间常数τ在深度为200 m处出现较为明显的假极值，当权为λ1=λ2=1时，如图 11b和图 12b所示，频率相关系数c的反演结果可以反映理论模型的趋势，但误差比加权系数λ1 < λ2时大，时间常数τ的反演结果不能正确反映模型原始趋势，并在深度为200 m处出现假极值.
当λ1 < λ2时，虚部在反演中起到的作用增大，频率相关系数c和时间常数τ的反演结果都可以很好地反映理论模型原始趋势.分别对比图 6c、6d和12c、12d，可以发现加权系数λ1=1，λ2=100000时的反演结果要稍好于加权系数λ1=1，λ1=100时的反演结果，推测这是由于当加权系数λ1=1，λ2=100000时，虚部在反演中起到的作用更大所引起的.对于频率相关系数c，在两层模型中，对比图 5c与5d，我们发现在权λ1=1，λ2=100时的反演结果要稍好于权λ1=1，λ2=100000时的反演结果，在三层模型中，通过对比图 11c与11d，我们同样可以发现，在权λ1=1，λ2=100的反演结果中，虽然最后一层的值与真值相差稍大，但前两层更接近真值，由此可知并不是λ2相对于λ1越大频率相关系数的反演结果就越好.
以上分析说明视复电阻率的虚部对频率相关系数c和时间常数τ更为敏感.反演频率相关系数c和时间常数τ时应对虚部的加权系数更大，以得到较好的反演结果.
3 结论本文在考虑电磁耦合的前提下，开展了1D频谱激电实、虚部加权反演研究，取得如下结论：
(1) 视复电阻率的实部和虚部对不同Cole-Cole参数的敏感程度不同，视复电阻率的实部对零频电阻率ρ0和极化率m较为敏感，视复电阻率的虚部对频率相关系数c和时间常数τ较为敏感，其中视复电阻率的虚部对时间常数τ的敏感程度比对频率相关系数c的敏感程度更低.
(2) 采用分两步反演的方法，首先固定频率相关系数c和时间常数τ的值不变，反演零频电阻率ρ0和极化率m，接着将反演得到的零频电阻率ρ0和极化率m的值作为各自的初始值，且固定不变，然后再反演频率相关系数c和时间常数τ，并通过理论模型的试算，验证了这种方法是可行的.
(3) 利用视复电阻率的实部和虚部加权反演，通过用加权的方法控制反演时实、虚部所占的比重，通过模型试算，我们发现不同加权系数得到的结果是不同的，且反演零频电阻率ρ0和极化率m时，应对实部加较大权，反演频率相关系数c和时间常数τ时，应对虚部加较大权，该算法可以较好的反演零频电阻率ρ0和极化率m，并且有效提高了频率相关系数c和时间常数τ的反演效果，结果更为可靠.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Brown R J. 1985. EM coupling in multi-frequency IP and a generalization of the Cole-Cole impedance model[J]. Geophysical Prospecting, 33(2): 282–302. DOI:10.1111/gpr.1985.33.issue-2 |
| [] | Fan C S, Li T L, Yan J Y. 2012. Research and application experiment on 2.5D SIP inversion[J]. Chinese Journal of Geophysics (in Chinese), 55(12): 4044–4050. DOI:10.6038/j.issn.0001-5733.2012.12.016 |
| [] | Ghorbani A, Camerlynck C, Florsch N. 2009. CR1Dinv:A Matlab program to invert 1D spectral induced polarization data for the Cole-Cole model including electromagnetic effects[J]. Computers & Geosciences, 35(2): 255–266. |
| [] | Ingeman-Nielsen T, Baumgartner F. 2006. CR1Dmod:A Matlab program to model 1D complex resistivity effects in electrical and electromagnetic surveys[J]. Computers & Geosciences, 32(9): 1411–1419. |
| [] | Liu B, Li S C, Li S C, et al. 2012. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization[J]. Chinese Journal of Geophysics (in Chinese), 55(1): 260–268. DOI:10.6038/j.issn.0001-5733.2012.01.025 |
| [] | Liu S, Guan S Y, Gao P F. 1994. Joint sip inversion for estimation of intrinsic Cole-Cole parameters of a polarizable ellipsoid[J]. Acta Geophysica Sinica (in Chinese), 37(S2): 542–551. |
| [] | Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data[J]. Chinese Journal of Geophysics (in Chinese), 56(12): 4278–4287. DOI:10.6038/cjg20131230 |
| [] | Liu Y L, Li T L, Hu Y C, et al. 2015. Fast quasi-linear approximation and the three-dimensional spectrum of induced polarization inversion study[J]. Chinese Journal of Geophysics (in Chinese), 58(12): 4709–4717. DOI:10.6038/cjg20151231 |
| [] | Luo Y Z, Fang S. 1986. An approximate inversion of the apparent complex resistivity spectrum[J]. Earth Science-Journal of Wuhan College of Geology, 11(1): 93–102. |
| [] | Major J, Silic J. 1981. Restrictions on the use of Cole-Cole dispersion models in complex resistivity interpretation[J]. Geophysics, 46(6): 916–931. DOI:10.1190/1.1441230 |
| [] | Pelton W H, Ward S H, Hallof P G, et al. 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP[J]. Geophysics, 43(3): 588–609. DOI:10.1190/1.1440839 |
| [] | Shang T X. 2008. 1D full-region inversion of CSAMT of bipolar source (in Chinese)[MSc. thesis]. Changchun:Jilin University. |
| [] | Sunde E D. 1968. Earth Conduction Effects in Transmission Systems[M]. New York: Dover Publication. |
| [] | Yu L B. 2010. The study of simulation and inversion of 1D spectrum induced polarization method (in Chinese)[MSc. thesis]. Changchun:Jilin University. |
| [] | Zhang S Z, Li Y X, Zhang S C, et al. 1984. The low frequency electrical phase spectra of mineralized rocks (ores) and some factors which influence them in some sulfide ore deposits in China[J]. Acta Geophysica Sinica (in Chinese), 27(2): 176–189. |
| [] | Zhou F. 2015. 1D and 2.5D modeling and inversion of complex resistivity data (in Chinese)[MSc. thesis]. Nanchang:East China institute of Technology. |
| [] | 范翠松, 李桐林, 严加永. 2012. 2.5维复电阻率反演及其应用试验[J]. 地球物理学报, 55(12): 4044–4050. DOI:10.6038/j.issn.0001-5733.2012.12.016 |
| [] | 刘斌, 李术才, 李树忱, 等. 2012. 基于不等式约束的最小二乘法三维电阻率反演及其算法优化[J]. 地球物理学报, 55(1): 260–268. DOI:10.6038/j.issn.0001-5733.2012.01.025 |
| [] | 刘崧, 官善友, 高鹏飞. 1994. 求极化椭球体真Cole-Cole参数的联合谱激电反演[J]. 地球物理学报, 37(S2): 542–551. |
| [] | 刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究[J]. 地球物理学报, 56(12): 4278–4287. DOI:10.6038/cjg20131230 |
| [] | 刘永亮, 李桐林, 胡英才, 等. 2015. 快速拟线性近似方法及三维频谱激电反演研究[J]. 地球物理学报, 58(12): 4709–4717. DOI:10.6038/cjg20151231 |
| [] | 罗延钟, 方胜. 1986. 视复电阻率频谱的一种近似反演方法[J]. 地球科学-武汉地质学院学报, 11(1): 93–102. |
| [] | 尚通晓. 2008. 双极源CSAMT一维全区反演[硕士论文]. 长春: 吉林大学. |
| [] | 于立波. 2010. 一维频谱激电法正反演研究[硕士论文]. 长春: 吉林大学. |
| [] | 张赛珍, 李英贤, 张树椿, 等. 1984. 我国几个金属矿区岩 (矿) 石的低频电相位频率特性及其影响因素[J]. 地球物理学报, 27(2): 176–189. |
| [] | 周峰. 2015. 复电阻率1D/2. 5D正反演研究[硕士论文]. 南昌: 东华理工大学. |