Comparative study on two kinds of electrode grounding resistance measurement methods in geo-resistivity observation
-
摘要: 利用传统三极接地电阻率测量方法对井下电极接地电阻进行测量会存在较大的测量误差，表明该方法已不适用于井下电极接地电阻测量。针对此问题，提出了利用台站现有地电阻率仪器进行电极接地电阻测量的新方法—直接测量法，并使用传统三极法和直接测量法对北京平谷地震台井下观测供电电极接地电阻进行了测量。测量结果表明：与传统三极法对比，直接测量法的测量结果更为准确、可靠，而且节省人力、物力和时间。Abstract: When using the traditional three-pole method to measure the underground electrode grounding resistance, there will be a large measurement error, which suggests that the method will not be suitable for the underground electrode grounding resistance measurement. In view of this problem, this paper proposes a new method for measuring the electrode grounding resistance by using the geo-electric resistivity observation system, which is so called direct measurement method. And then, the grounding resistances of the underground observation power supply electrodes of Pinggu seismic station in Beijing were measured by using the traditional three-pole method and the direct measurement method. The measurement results show that, compared with the traditional three-pole method, the direct measurement method can get more reliable measurement results with less error, and can save human and material resources and measuring time.
-
引言
在勘探地震学中，由于近地表（0—500 m）分布的沙漠、沼泽、戈壁和风化层等松软结构的影响，地震波在传播过程中，出现明显的衰减，给地球物理成像带来了很大干扰。传统的全波形反演方法忽略了衰减对波形的影响，但通过品质因子Q则可以判断地下介质对波形的衰减能力，Q值越高，衰减能力越弱。如果已知具体的Q值分布，则可提高速度反演的精度以及地球物理解释的可靠性（Pramanik et al，2000 ；Wang，2002；Zhu et al，2014 ）。除此之外，Q值也可以作为地下天然气存储的指示依据（Odebeatu et al，2006 ）。
由波动方程可知，地震波在地球介质中传播时，速度和衰减是耦合的（Kjartansson，1979；Aki，Richards，1980；Zhu，Carcione，2014），因此，利用波形同时反演速度和衰减难度较大。目前已经发展了较多的利用地震数据反演同时得到速度和衰减结构的方法，例如，Stewart（1983）利用频率域上行波和下行波的波形同时反演出速度和地下介质品质因子Q值。Q值模型的反演需要以较准确的速度数据为基础，才能从多重外在因素所导致的衰减中独立出固有Q值的作用，故Kamei和Pratt （2008）率先提出利用波形首先反演出纵波速度结构，然后进行速度和Q值联合反演的方法。Smithyman等（2009）则利用声波波形反演了二维P波速度及其衰减。Zhu和Harris （2015）利用井间地震走时数据联合反演纵波速度、横波速度和Q值。
地球物理学方面的联合反演始于Vozoff和Jupp （1975），其利用直流电磁测深数据和大地电磁测深数据联合反演一维层状电阻率模型。自此之后，联合反演在地震数据领域得到了广泛的发展。杨文采和焦富光（1987）利用阻尼最小二乘法联合反射波走时和均方根速度反演出速度结构；张树林等（1993）利用井间地震和逆垂直地震剖面（vertical seismic profiling，简写为VSP）联合层析反演三维速度结构；Tryggvason等（2002）利用P波和S波走时信息联合反演纵横波速度，并用其进行地震定位；Gallardo和Meju （2003）提出利用交叉梯度联合反演的方法，并得到了广泛应用（Gallardo，Meju，2004，2007；Tryggvason，Linde，2006；Colombo et al，2008 ）。
通常我们会利用一维弹性波正演模拟数据反演一维速度结构，利用二维弹性波正演模拟数据反演二维结构，或者利用三维弹性波正演模拟数据反演三维结构（Zhou et al，1995 ；Pratt，Shipp，1999）。然而由于一维弹性波波形反演计算效率高，但由于近地表结构复杂，一维模型很难准确地描述实际的地质结构；二维或三维弹性波波形反演能够提供更加精确的地下结构，但计算耗时多，而且对计算机的存储有较高的要求，因此，如何利用一维弹性波模拟反演得到地下二维结构是值得探讨的。Auken和Christiansen （2004）提出伪二维反演电阻率的方法，即利用一维正演模拟，通过增加水平约束条件来反演二维结构。目前地球物理联合反演依然面临着不同数据权重选择的难题（Moorkamp et al，2007 ；Colombo et al，2016 ）。为了提高运算效率，降低计算所用存储空间，本文拟提出一种利用一维弹性波正演，通过二维Tikhonov正则化矩阵（Tikhonov，Arsenin，1977）对模型参数的约束，来反演二维结构的方法。该方法具体步骤为：首先，分析纵波速度、横波速度、QP和QS等4个参量联合反演的可行性，并设计应用交叉梯度联合反演的目标函数；其次，进行理论合成模型测试，讨论不同的权重因子对反演结果的影响；最后，对某地区的地震波数据进行处理，联合反演近地表的速度和Q值。
1. 弹性波正演
地震波是能量的传播，通常包括体波和面波两种。体波是指脉冲短且能够传到地下深部的波，包含P波和S波，遵循由地下不同区域不同的弹性模量和密度引起的反射、折射定律；面波通常指在地表传播的波，包括瑞雷波和勒夫波，其频率越低，能够到达的地下深度越大，而且面波的振幅随深度的增加衰减得很快。本文主要研究地震波中的P波和瑞雷波。
由于地球介质的黏弹性或内在摩擦，地震波能量在传播过程中逐渐衰减（Shearer，1999）。衰减的强度用1/Q来表示，即：
$\frac{1}{Q} {\text{＝}} {\text{－}} \frac{{\Delta E}}{{2\pi E}}{\text{，}}$
(1) 式中，E表示峰值应变能，−ΔE表示一个周期内能量的衰减，因此Q值越大，衰减越小。Knopoff （1964）认为高频地震波的Q值与频率有关，而低频地震波的Q值受频率影响则较小。
在地震波传播过程中，Q值通过
$v\left(\omega \right) {\text{＝}} {v_0}\left({1 {\text{＋}} \frac{1}{{\pi Q}}\ln \left({\frac{\omega }{{2\pi }}} \right) {\text{－}} \frac{\rm i}{{2Q}}} \right)$
(2) 与速度结合，影响弹性波传播（Aki，Richards，1980）。式中，v（ω）表示频率为ω时的速度值，v0表示参考频率为f0时的速度值，故非线性波形正演计算可以简化为
${{d}} {\text{＝}} {G}\left({{m}} \right){\text{，}}$
(3) 式中：d为计算出的波形数据；m为模型参数向量，包含vP，vS，QP和QS；G表示正演计算函数。随后，利用有限差分公式分别求取弹性波对vP，vS，QP和QS的偏微分，并用波形将敏感度体现出来，其表达式为
$\frac{{\partial { d}\left({{ m_i}{\text{,}} t} \right)}}{{\partial {m_i}}} {\text{＝}} \frac{{{ d}\left({{ m_i} {\text{＋}} \delta { m_i}{\text{,}} t} \right) {\text{－}} { d}\left({{m_i}{\text{,}} t} \right)}}{{\delta { m_i}}}{\text{，}}$
(4) 式中，d（mi，t）表示t时刻的地震记录，mi表示第i个模型的参数，δmi表示对mi的扰动。用扰动后模型正演得出的波场减去扰动前模型正演得出的波场，再除以扰动量，便可得到波场对扰动模型的偏微分。
2. 弹性波反演
弹性波的初至波含有vP和QP信息，面波含有vS和QS的信息（Wang，Zhang，2017）。勘探地球物理中，初至波包含详细的地下结构信息，而面波对浅层结构信息更敏感。因此，综合初至波信息和面波信息，可以同时反演vP，vS，QP和QS。尽管目前的研究尚未发现QP与QS之间存在解析关系，但由于在地质学分析中，二者在结构上存在一定的联系，对于干燥或部分饱和的岩石，QP和QS的值的变化趋势相同（Winkler，Nur，1979），所以，我们考虑在联合反演过程中增加结构约束，使QP与QS的模型在结构上保持一致，并与速度模型相近。因此，所使用交叉梯度的目标函数为
$\begin{aligned} \min \Big[ {\varphi \big({m} \big)} \Big] {\text{＝}}& \min \Big( {{\varphi _{{d}}} {\text{＋}} {\varphi _{{L}}} {\text{＋}} {\varphi _t}} \Big) {\text{＝}} \min \Biggl\{ \frac{1}{2}{\Bigg\|\!\! {\begin{array}{*{20}{c}} {{\omega _{\rm{P}}}\big[{{{d}_{\rm{P}}} {\text{－}} {{G}_{\rm{P}}}\big({{{ m}_{{v_{\rm{P}}}}}, {{ m}_{{Q_{\rm{P}}}}}}\big)} \big]} \\ {{\omega _{\rm{S}}}\big[{{{d}_{\rm{S}}} {\text{－}} {{G}_{\rm{S}}}\big({{{ m}_{{v_{\rm{S}}}}}, {{ m}_{{Q_{\rm{S}}}}}} \big)} \big]} \end{array}} \!\!\Bigg\|^2} {\text{＋}}\big. \\& \big.\frac{1}{2}{\Bigg\|\!\! {\begin{array}{*{20}{c}} {{\alpha _{{v_{\rm{P}}}}}{{L}}{{ m}_{{v_{\rm{P}}}}}} \\ {{\alpha _{{v_{\rm S}}}}{{L}}{{{ m}}_{{v_{\rm{S}}}}}} \end{array}} \!\!\Bigg\|^2} {\text{＋}} \frac{1}{2}{\Bigg\|\!\! {\begin{array}{*{20}{c}} {{\alpha _{{Q_{\rm{P}}}}}{{L}}{{ m}_{{Q_{\rm{P}}}}}} \\ {{\alpha _{{Q_{\rm{S}}}}}{{L}}{{ m}_{{Q_{\rm{S}}}}}} \end{array}}\!\! \Bigg\|^2} {\text{＋}} \frac{1}{2}{\alpha _t}{\left\| {t\big({{{ m}_{{v_{\rm{P}}}}}, {{ m}_{{Q_{\rm{P}}}}}, {{ m}_{{v_{\rm{S}}}}}, {{ m}_{{Q_{\rm{S}}}}}} \big)} \right\|^2} \Biggr\} {\text{，}}\end{aligned}$
(5) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
$\nabla \varphi ({m}) {\text{＝}} {{{\textbf{\textit{0}}}}}{\text{，}}$
(6) 对式（6）添加扰动，可以得到
${\nabla ^2}\varphi \left({m} \right) \bullet \delta {m} {\text{＝}} {\text{－}} \nabla \varphi \left({m} \right){\text{，}}$
(7) 根据式（7），可得到反演过程中模型更新量的公式，即
\begin{equation}\begin{array}{l}\left({\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}\!\!\!\!\!\!\!\!{\omega _{\rm{P}}^2{J}_{{v_{\rm{P}}}}^{\rm{T}}{{J}^{}_{{v_{\rm{P}}}}} {\text{＋}} \alpha _{{v_{\rm{P}}}}^2{{L}^{\rm{T}}}{L}}&{\omega _{\rm{P}}^2{J}_{{v_{\rm{P}}}}^{\rm{T}}{{J}_{{Q_{\rm{P}}}}}}\\{\omega _{\rm{P}}^2{J}_{{Q_{\rm{P}}}}^{\rm{T}}{{J}_{{v_{\rm{P}}}}}}&{\omega _{\rm{P}}^2{J}_{{Q_{\rm{P}}}}^{\rm{T}}{{J}_{{Q_{\rm{P}}}}} {\text{＋}} \alpha _{{Q_{\rm{P}}}}^2{{L}^{\rm{T}}}{L}}\end{array}}&{\begin{array}{*{20}{c}}{\textbf{\textit{0}}}&\quad \quad \quad \quad \quad\quad {\textbf{\textit{0}}}\\[4.5pt]{\textbf{\textit{0}}}&\quad \quad \quad \quad \quad\quad {\textbf{\textit{0}}}\end{array}}\\{\begin{array}{*{20}{c}}{\textbf{\textit{0}}}&\quad \quad \quad \quad \quad\quad \quad {\textbf{\textit{0}}} \\{\textbf{\textit{0}}}&\quad \quad \quad \quad \quad \quad\quad {\textbf{\textit{0}}}\end{array}}&{\begin{array}{*{20}{c}}{\omega _{\rm{S}}^2{J}_{{v_{\rm{S}}}}^{\rm{T}}{{J}_{{v_{\rm{S}}}}} {\text{＋}} \alpha _{{v_{\rm{S}}}}^2{{L}^{\rm{T}}}{L}}&{\omega _{\rm{S}}^2{J}_{{v_{\rm{S}}}}^{\rm{T}}{{J}_{{Q_{\rm S}}}}}\\{\omega _{\rm{S}}^2{J}_{{Q_{\rm{S}}}}^{\rm{T}}{{J}_{{v_{\rm{S}}}}}}&{\omega _{\rm{S}}^2{J}_{{Q_{\rm{S}}}}^{\rm{T}}{{J}_{{Q_{\rm{S}}}}} {\text{＋}} \alpha _{{Q_{\rm{S}}}}^2{{L}^{\rm{T}}}{L}}\!\!\!\!\!\!\!\!\!\!\end{array}}\end{array}} \right)\!\left(\!\!\!{\begin{array}{*{20}{c}}{\delta {{ m}_{{v_{\rm{P}}}}}}\\{\delta {{ m}_{{Q_{\rm{P}}}}}}\\{\delta {{ m}_{{v_{\rm{S}}}}}}\\{\delta {{ m}_{{Q_{\rm{S}}}}}}\!\!\!\!\!\end{array}} \!\!\right)\!\! {\text{＋}}\\ \quad\quad\quad\quad\quad {\alpha _t}{{B}^{\rm{T}}}{B}\left({\begin{array}{*{20}{c}}{\delta {{ m}_{{v_{\rm{P}}}}}}\\{\delta {{ m}_{{Q_{\rm{P}}}}}}\\{\delta {{ m}_{{v_{\rm{S}}}}}}\\{\delta {{ m}_{{Q_{\rm{S}}}}}}\end{array}} \!\!\right) {\text{＝}} \left({\begin{array}{*{20}{c}}{\omega _{\rm{P}}^2{J}_{{v_{\rm{P}}}}^{\rm{T}}\left({{{d}_{\rm{P}}}{\text{－}} {{G}_{\rm{P}}}} \right)}\\{\omega _{\rm{P}}^2{J}_{{Q_{\rm{P}}}}^{\rm{T}}\left({{{d}_{\rm{P}}} {\text{－}} {{G}_{\rm{P}}}} \right)}\\{\omega _{\rm{S}}^2{{J}}_{{v_{\rm{S}}}}^{\rm{T}}\left({{{d}_{\rm S}} {\text{－}} {{G}_{\rm S}}} \right)}\\{\omega _{\rm{S}}^2{{J}}_{{Q_{\rm{S}}}}^{\rm{T}}\left({{{d}_{\rm S}} {\text{－}} {{G}_{\rm S}}} \right)}\end{array}} \right) {\text{－}} \left({\begin{array}{*{20}{c}}{\alpha _{{v_{\rm{P}}}}^2{{L}^{\rm{T}}}{L}{{ m}_{{v_{\rm{P}}}}}}\\{\alpha _{{Q_{\rm{P}}}}^2{{L}^{\rm{T}}}{L}{{ m}_{{Q_{\rm{P}}}}}}\\{\alpha _{{v_{\rm{S}}}}^2{{L}^{\rm{T}}}{L}{{ m}_{{v_{\rm S}}}}}\\{\alpha _{{Q_{\rm{S}}}}^2{{L}^{\rm{T}}}{L}{{ m}_{{Q_{\rm{S}}}}}}\end{array}} \right) {\text{－}} {\alpha _t}{{B}^{\rm{T}}}t{\text{，}}\end{array}\end{equation}
(8) 式中，JP和JS分别为GP和GS的雅克比矩阵，可以通过式（4）计算出来。B为交叉梯度对vP，QP，vS和QS的偏微分矩阵，则
${B} {\text{＝}} \left({\begin{array}{*{20}{c}} {\displaystyle\frac{{\partial t}}{{\partial {{ m}_{{v_{\rm P}}}}}}, }&{\displaystyle\frac{{\partial t}}{{\partial {{ m}_{{Q_{\rm P}}}}}}, }&{\displaystyle\frac{{\partial t}}{{\partial {{ m}_{{v_{\rm S}}}}}}, }&{\displaystyle\frac{{\partial t}}{{\partial {{ m}_{{Q_{\rm S}}}}}}} \end{array}} \right).$
(9) 目前已发展了多种算法来反演式（8），例如模拟退火法、遗传算法、奇异值分解法、共轭梯度法等。利用全局最优化算法或半全局最优化算法，可以得到方程的最优解，但计算量大，所需计算机的存储空间大，在处理少量模型参数的情况下有优势，例如一维模型的反演；但在实际勘探地震学中，大量的波形数据以及模型参数需要解决，如果使用模拟退火法、遗传算法等会极大地增加计算时间（Beaty et al，2002 ）；而共轭梯度法在求解方程解时，仅需要动态保存雅克比矩阵和反演的模型参数，在处理勘探地震数据时更加有效实用（Hestenes，Stiefel，1952）。因此本文采用共轭梯度法求解式（8）得到模型的更新量。最后，通过迭代反演的模型正演计算的数据与观测到的数据求取差值，当误差达到预设的条件时，迭代停止；或者，当迭代次数达到预设迭代次数N时，反演结束。从目标函数的表达式式（5）可以看出，波形数据的权重以及交叉梯度的权重直接对反演产生作用，因此，我们将在后面的理论模型测试中详细讨论不同权重因子对反演结果的影响。
3. 理论模型测试
建立一个层状模型，如图1所示。vP变化范围为1 900—2 200 m/s，vS根据vP/vS＝2.0得到，故其变化范围为950—1 100 m/s，QP的取值范围为30—60，QS的取值范围为20—60。模型具有明显的分层：在180 m以下为均匀半空间；第二层在水平方向350—550 m之间和850—1 050 m之间设置一凸起的异常体和一凹陷的异常体。采用震源时间函数为雷克子波，中心频率为8 Hz，震源深度设置在10 m，接收器位于地表，偏移距为800 m，且记录了2.048 s的地震波序列，采样率为0.004 s。为了使联合反演过程中数据对速度和衰减的响应能够合理地体现出来，设置速度和衰减的参考量（vP0，vS0，QP0和QS0）来平衡波形的敏感度。通过测试得到一组比较合理的参考量为vP0＝100 m/s；vS0＝100 m/s；QP0＝1；QS0＝10。初至波和面波对vP/vP0，vS/vS0，QP/QP0以及QS/QS0的敏感度，如图2所示。
![]()
图 2 弹性波对参数vP/vP0 （a），vS/vS0 （b），QP/QP0 （c） 和QS/QS0 （d）的敏感度Figure 2. Sensitivity of elastic waveform to parameters vP/vP0 （a），vS/vS0 （b），QP/QP0 （c） and QS/QS0 （d）![]()
图 1 真实理论测试模型（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 1. True models for synthetic tests（a） vP model；（b） vS model；（c） QP model；（d） QS model建立如图3所示的联合反演初始模型。通常波形反演对初始速度模型的要求较高，因此，我们根据观测数据的初至波到时建立初始P波速度模型，然后根据纵横波速度比为2，得到初始S波速度模型。由于初始模型正演模拟得到的波形与真实波形之间的相位差不能过大，否则波形反演会陷入周期跳跃问题，不可得到准确的反演结果，又因波形对纵横波衰减的敏感度较小，联合反演对QP和QS初始模型的要求则较低，因此我们建立如图3c和图3d所示的各向均匀的衰减模型作为初始QP和QS模型。
![]()
图 3 联合反演的初始输入模型（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 3. Initial models for joint inversion（a） vP model；（b） vS model；（c） QP model；（d） QS model随后测试不同权重因子对反演结果的影响。首先，确定初至波和面波的波形数据权重为1，vP，vS，QP和QS模型的光滑因子均为0.1，测试不同交叉梯度权重αt对反演结果的影响。当αt分别为0，0.1，1，3和10的情况下，波形不匹配度随反演迭代次数的变化如图4所示。可以看出，5组参数下的联合反演均收敛，而且随着权重因子的增大，交叉梯度对模型的约束从反演结果来看越来越强，在联合反演中的作用越来越大，致使数据的不匹配度可收敛至很小。
![]()
图 4 不同交叉梯度权重下初至波（a）和面波（b）的波形不匹配度随迭代次数的变化Figure 4. Misfit variation of initial arrival data （a） and surface wave data （b） with iteration number by different cross-gradient weightingsαt＝0表示联合反演中不加入交叉梯度算子的限制，经过20次迭代反演得到的速度模型及其衰减模型如图5所示。横波速度结构的异常体反演清晰（图5a）；而纵波速度结构异常体的轮廓模糊，出现假构造，在垂直方向的分辨率很低，例如深度约为500 m时第二层的凸起异常区延伸至地表和第三层（图5b）；QP模型与真实模型的相差较大（图5c）；QS模型第二层的起伏变化明显，但第三层与第二层的分界模糊（图5d）。由于未引入交叉梯度算子，横波信息的反演结果比纵波信息的反演结果更接近真实模型，所以我们结合图4中波形不匹配度的变化认为，αt＝10的反演结果要优于αt＝0时的结果。对比αt＝10时反演的速度及衰减模型（图6）与无交叉梯度的反演结果（图5）可以看到：vP在图6a中第二层有起伏变化；vS在图6b中的分层较图5b中更加清晰，第二层的起伏更加明显；QP和QS在图6c和d也显示出第二层具有明显的起伏变化且在深度为100 m处形成分界线。此后在进行真实波形与初始模型的波形对比以及真实波形与经过20次迭代后反演所得波形的对比（图7），可见初至波与面波拟合得较好。
![]()
图 5 交叉梯度权重αt= 0时理论数据的联合反演结果（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 5. The synthetic data inversion results with cross-gradient weighting αt= 0（a） vP model；（b） vS model；（c） QP model；（d） QS model![]()
图 7 理论模型测试的初至波和面波的真实数据（黑线）与模拟数据（红线）的对比图（a）中红线表示基于初始模型正演计算的数据，图（b）中红线表示基于交叉梯度为αt=10反演模型正演计算的数据，下同Figure 7. The comparison of true data （black line） and calculated data （red line） for early arrivals and surface waves in the synthetic testsThe red line in Fig.（a） represents the forward calculation date associ-ated with the initial model，the red line in Fig.（b） represents the forward calaulation data resulted from inversion model with cross-gradient weighting αt=10，the same below![]()
图 6 交叉梯度权重αt=10时理论数据的联合反演结果（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 6. The synthetic data inversion results with cross-gradient weighting αt=10（a） vP model；（b） vS model；（c） QP model；（d） QS model由前文得到的反演结果可知，交叉梯度因子能够为反演模型提供结构上的约束，使反演的4个模型在结构上趋于一致。联合反演的模型中，有关S波信息的横波速度和衰减结构与真实模型的拟合度较高，而有关P波信息的纵波速度和衰减结构与真实模型拟合度较低。因此，下文将简单讨论不同波形数据的权重因子对联合反演结果的影响。
确定交叉梯度的权重为10，S波的权重为1，分别设置ωP为0.2，0.5，1.0和5.0等4种不同权重值进行联合反演，计算反演结果与真实模型之间的误差，结果如图8所示。由图可见：随着P波数据权重的增加，反演的P波速度结构与真实结构的差异逐渐变小，但当ωP＝5.0时，差异在第5次迭代时最小，而随着迭代次数的增加而变大（图8a）；S波速度的反演结果随着ωP的改变无明显变化，与真实结构之间的差异降至较低水平（图8b）；当ωP＝1.0时，反演得到的QP模型与真实模型之间的差异达到最小（图8c）；经过20次迭代，ωP＝0.5时的QS模型与真实模型之间差异最小，但ωP＝1.0的结果在可接受的范围内（图8d）。综合对比vP，vS，QP和QS与真实模型之间的差异，可知ωP＝1.0时的反演结果较好。
![]()
图 8 ωP取不同值时反演的模型与真实模型之间的差异（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 8. Differences between the true model and the inverted model with different ωp（a） vP model；（b） vS model；（c） QP model；（d） QS model由图1—8的理论模型测试结果可知，联合初至波和面波反演速度与衰减参数的方法收敛速度较快，而且伪二维反演算法能够处理含有微弱起伏变化的近似水平层状结构。通过测试不同的交叉梯度权重对反演结果的影响（图4—8）可得：随着αt的增加，反演的4个模型在结构上越来越相似，有效地弥补了波形数据对QP和QS敏感度低的缺陷。尽管得到的衰减模型较为光滑，且QP值和QS值与真实模型存在差异，但对于近地表地区，分辨率较低的QP和QS结构也足够为我们提供分析近地表的地质构造信息。通过对不同数据权重因子的测试（图8）表明：当减小P波数据的权重时，S波数据占主导地位，S波的相关参量vS和QS，均能反演出较为准确的结构，而P波的相关参量vP和QP的模型更新量受交叉梯度的限制作用变大，趋向于与vS和QS结构上的一致，但是忽略P波数据的拟合度，容易导致反演的vP和QP模型不合理；当增大P波数据的权重时，P波数据在反演中占主导地位，而S波数据的拟合被压制，因此，在图8a和8c中，ωP＝0.5时分别在第5次和第7次迭代横型的误差中降至最小值。然而，由于初至波对vP和QP的分辨率较低，反演的结构与真实结构存在差异，因此交叉梯度的作用引起vS和QS结构随着不准确的vP和QS结构变化，最终导致反演的模型拟合度不高。理论模型测试结果表明，根据反演的数据不匹配度曲线和模型的差异曲线选择最优权重因子，有利于得到合理可靠的反演结果。
4. 实际数据应用
将伪二维联合反演近地表的速度和衰减结构的方法应用到实际数据中。首先，根据到时文件拾取数据的初至到时，利用射线追踪走时反演方法得到一个初始纵波速度模型；其次，对地震波数据进行滤波、截取波形等预处理，滤波范围是2—15 Hz。本文测试了几种频段，包括2—10 Hz，2—15 Hz，2—20 Hz以及2—30 Hz，其中利用2—15 Hz频段滤波的波形反演得到的模型较合理；最后，选取偏移距为500 m、间隔为80 m的地震记录，炮点埋深为41 m，检波器置于地表，假定纵横波速度比为1.732，根据走时反演得到的纵波速度，得到一个横波速度模型。由于近地表浅层存在风化层，沉积层土壤松软，浅层的速度值很低，不完全遵循假定的纵横波速比，因此，我们根据面波到时，不断修改横波速度进行正演计算，直到正演模拟的面波到时与观测数据的面波到时基本吻合。由于地震波对QP和QS敏感度较低，我们建立均匀的QP和QS模型作为初始模型，如图9所示。
![]()
图 9 实际测试数据的初始输入模型（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 9. Initial models for real data tests（a） vP model；（b） vS model；（c） QP model；（d） QS model联合反演中不同交叉梯度权重下，数据不匹配度随迭代次数的变化如图10所示。可见：随着αt的增加，数据的不匹配度变大，但基本随着迭代次数的增加而降低；当αt小于10时，数据的不匹配度基本保持在相同水平；当αt为10时，数据的匹配度和交叉梯度的权重基本平衡了；当αt为100时，数据的不匹配度陡然增高。αt为0和10时联合反演的速度及其衰减模型，如图11和图12所示，对比可知，图12反演的QP和QS模型与速度在结构上有很大相似度。基于初始输入模型计算的波形与图12中真实波形的对比以及基于联合反演结果正演计算的波形与真实波形的对比如图13所示。结果显示，初至波拟合得较好，而面波由于实际采集的数据质量较差，最终拟合效果较初至波差。
![]()
图 11 交叉梯度权重αt= 0时实际数据的联合反演结果（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 11. The real data inversion results with cross-gradient weighting αt= 0（a） vP model；（b） vS model；（c） QP model；（d） QS model![]()
图 12 交叉梯度权重αt=10时实际数据的联合反演结果（a） vP模型；（b） vS模型；（c） QP模型；（d） QS模型Figure 12. The real data inversion results with cross-gradient weighting αt=10（a） vP model；（b） vS model；（c） QP model；（d） QS model![]()
图 13 实际数据测试的初至波和面波的真实数据（黑线）与基于初始模型（a）和反演模型（b）模拟数据（红线）的对比Figure 13. The comparison of true data （black line） and calculated data （red line） for early arrivals and surface waves associated with initial models （a） and inversion models （b）![]()
图 10 不同交叉梯度权重下初至波（a）和面波（b）的不匹配度随迭代次数的变化Figure 10. Misfit variation of initial arrival data （a） and surface wave data （b） with iteration number by different cross-gradient weightings研究区域地表风化严重，沉积层地质结构疏松，因此浅层介质对地震波吸收较多，尤其是地表至60 m深度的地层，Q值较小。由于面波的信噪比较低，因此，在反演过程中可利用交叉梯度的作用，得到与纵波速度和QP结构类似的横波速度和QS模型。通过对不同交叉梯度权重因子的讨论（图10），结果显示：如果交叉梯度权重过小，则4个参量之间的结构限制作用变小，联合反演结果与不引入交叉梯度的反演结果类似；如果交叉梯度的权重过大，则联合反演被交叉梯度因子主导从而忽略数据对反演的作用。因此，在联合反演中，要通过对比不同交叉梯度权重作用下的反演结果，综合分析波形的拟合度，以求得到最优的权重取值。
5. 讨论与结论
联合初至波和面波反演近地表的速度及其衰减，研究结果表明，Q值的影响在波形反演中不可被忽视，尤其是地质结构复杂的近地表区域，介质的Q值越小，对波形的衰减作用越大。利用弹性波的初至波和面波联合反演vP，vS，QP和QS等4个参量，不仅可以为地震解释提供丰富的结构信息，而且可以降低地球物理解释的非唯一性。本文在反演的目标函数中，引入交叉梯度算子，利用其对vP，vS，QP和QS这4个参量在结构上的约束，实现了弹性波全波形反演近地表的速度及其衰减结构。
在理论模型测试中，本文讨论了目标函数中的不同数据权重因子和交叉梯度权重因子对反演结果的影响。如果交叉梯度权重因子过大，反演的4个参数虽然在结构上具有较高的一致性，但是不能体现出地震波数据的影响，致使反演的结果可靠性降低；如果交叉梯度权重因子过小，反演的结果对波形数据的质量依赖性较强，而且波形对衰减的敏感度远低于对速度的敏感度，造成在反演迭代过程中只更新速度结构而衰减结构无明显的变化，反演结果也不可靠。因此，适当提高交叉梯度的权重，而又不忽略地震波数据在联合反演中的主导地位，导致联合反演得到的4个参量更接近真实模型。实际数据处理过程的难点在于对数据和交叉梯度权重因子的选择。由于研究区域近地表风化严重，浅层的衰减强，弹性波数据，尤其面波数据，质量较差，因此，在反演中适当降低面波数据的权重，同时增加交叉梯度的权重，使反演的横波速度和QS倾向于与纵波速度和QP的结构一致，有利于得到可靠的结果。实际操作中，需要通过大量的测试，依据反演结果与真实数据之间的拟合度，选择合理的权重值。
本文提出的伪二维反演是基于一维共中心点假设，利用二维正则化矩阵增加反演模型网格点之间的联系，实现二维结构反演。对于第4节中实际数据的处理，采用15个核心并行计算，反演迭代一次约需6分钟，20次迭代需2个小时，而真正的二维弹性波波形反演迭代一次的时间远大于1个小时。这种算法可以处理近似水平层状结构或者有微小起伏变化的结构。当地震波在一维介质中传播时，地表接收到的波形可以用来反演震源和接收器之间的中心点下方的结构，而如果速度或者Q模型在水平方向上有明显的变化，那么这种利用一维正演反演二维结构的伪二维算法可能会失败。然而，这种算法相比于一维反演结果能够提供更加丰富的地下结构水平方向的变化，相比于真正的二维弹性波反演，提高了计算效率，节省了计算时间。因此，这种算法在处理数据量大的勘探地球物理反演问题方面有很大的发展前景。
-
![]()
图 1 半球形电极接地电阻示意图
Figure 1. Schematic diagram of grounding resistance of a hemispherical electrode
![]()
图 2 直线三极法测量接地电阻示意图
Figure 2. Schematic diagram of measuring grounding resistance by linear three-electrode method
![]()
图 3 测量误差系数c随电位极P位置的变化曲线
Figure 3. Variation of error coefficient e with position of the electrode P changing
![]()
图 4 测量井下电极接地电阻示意图
Figure 4. Schematic diagram of measuring grounding resistance of an underground electrode
![]()
图 5 电极接地电阻相对误差随井下电极埋深的变化
Figure 5. Variation of relative error of electrode grounding resistance with depth of electrodes
![]()
图 6 利用直接测量法测量井下电极接地电阻
（a） 测量方法示意图；（b） 等效电路图
Figure 6. Measuring grounding resistance of underground electrodes with direct measurement method
（a） Schematic diagram of measurement method；（b） Equivalent circuit diagram
![]()
图 7 平谷地震台地电阻率井下布极示意图
Figure 7. Schematic diagram for distribution of underground borehole electrodes at Pinggu seismic station
![]()
图 8 三极法测量A电极接地电阻示意图
Figure 8. Schematic diagram for measuring grounding resistance of the electrode A by the three-pole method
表 1 不同土壤电阻率ρs、不同电极埋深h和不同电极尺寸条件下的接地电阻测量误差e和相对误差δ
Table 1 Measurement error e and relative error δ of grounding resistance with different soil resistivity ρs，electrode buried depth h and electrode size
电极半径r0
/mh＝100 m h＝200 m ρs＝100 Ω•m ρs＝200 Ω•m ρs＝100 Ω•m ρs＝200 Ω•m e/Ω δ e/Ω δ e/Ω δ e/Ω δ 0.2 0.492 1.24% 0.984 1.24% 0.638 1.60% 1.276 1.60% 0.5 0.491 3.08% 0.982 3.08% 0.638 4.01% 1.276 4.01% 1.0 0.492 6.17% 0.984 6.17% 0.638 8.02% 1.277 8.02% 
下载: 导出CSV
表 2 利用不同方法测量得到的井下供电极A和B的接地电阻
Table 2 Grounding resistance of the underground power supply electrodes A and B by different methods
测量次数 被测电极A的接地电阻/Ω 被测电极B的接地电阻/Ω 直接测量法 直线三极法 直接测量法 直线三极法 第一次 59.36 54.7 54.46 51.0 第二次 59.40 54.8 54.39 50.8 第三次 59.43 54.8 54.35 50.9 第四次 59.37 55.4 54.38 51.0 第五次 59.44 54.7 54.41 50.9 均值 59.40 54.9 54.40 50.9 均方差 0.03 0.28 0.04 0.08 实测差值 −4.52 −3.48 实测相对误差 −7.62% −6.39%
下载: 导出CSV
-
杜学彬. 2010. 在地震预报中的两类视电阻率变化[J]. 中国科学:地球科学,40(10):1321–1330. Du X B. 2010. Two types of changes in apparent resistivity in earthquake prediction[J]. Science China Earth Science,54(1):145–156.
冯慈璋. 1983. 电磁场[M]. 第二版. 北京: 高等教育出版社: 108–110. Feng C Z. 1983. Electromagnetic Field[M]. Second edition. Beijing: Higher Education Press: 108–110 （in Chinese）.
冯志伟,马金福. 2009. 电位降法测量接地电阻时电压极补偿点位置分析[J]. 电气应用,28(8):18–19. Feng Z W,Ma J F. 2009. Analysis of the position of the voltage pole compensation point when measuring grounding resistance by fall-of-potential method[J]. Electrotechnical Application,28(8):18–19 (in Chinese).
冯志伟. 2011. 影响接地电阻测量的因素分析[D]. 南京: 南京信息工程大学: 28–37. Feng Z W. 2011. Analysis of the Effect Factors of the Measurement of Grounding Resistance[D]. Nanjing: Nanjing University of Information Science and Technology: 28–37 （in Chinese）.
高曙德. 2016. 深井地电观测技术在地震监测中的应用探讨[J]. 地球物理学进展,31(5):2078–2088. doi: 10.6038/pg20160526 Gao S D. 2016. Discussion on the deep well geoelectric observation technique applied in earthquake monitoring[J]. Progress in Geophysics,31(5):2078–2088 (in Chinese).
钱家栋. 2010. 地震电磁学理论基础与观测技术[M]. 北京: 地震出版社: 293. Qian J D. 2010. Underpinnings and Observation Technology of Seismic Electro-Magnetism[M]. Beijing: Seismological Press: 293 （in Chinese）
钱家栋,马钦忠,李劭秾. 2013. 汶川MS8.0地震前成都台NE测线地电阻率异常的进一步研究[J]. 地震学报,35(1):4–17. doi: 10.3969/j.issn.0253-3782.2013.01.002 Qian J D,Ma Q Z,Li S N. 2013. Further study on the anomalies in apparent resistivity in the NE configuration at Chengdu station associated with Wenchuan MS8.0 earthquake[J]. Acta Seismologica Sinica,35(1):4–17 (in Chinese).
王洪泽, 杨丹, 王梦云. 2007. 电力系统接地技术手册[M]. 北京: 中国电力出版社: 15–17. Wang H Z, Yang D, Wang M Y. 2007. Manual of Grounding Technology of Power System[M]. Beijing: China Electric Power Press: 15–17 （in Chinese）.
王兰炜,张宇,张世中,颜蕊,王子影,张兴国,胡哲. 2015. 我国井下地电阻率观测技术现状分析[J]. 地震地磁观测与研究,36(2):95–102. doi: 10.3969/j.issn.1003-3246.2015.02.018 Wang L W,Zhang Y,Zhang S Z,Yan R,Wang Z Y,Zhang X G,Hu Z. 2015. The status of deep-well geo-electrical resistivity observation in China[J]. Seismological and Geomagnetic Observation and Research,36(2):95–102 (in Chinese).
解滔,于晨,卢军. 2019. 开展小极距井下地电阻率观测的可行性分析[J]. 中国地震,35(1):14–24. doi: 10.3969/j.issn.1001-4683.2019.01.002 Xie T,Yu C,Lu J. 2019. Feasibility analysis on short-electrode spacing well apparent resistivity observation[J]. Earthquake Research in China,35(1):14–24 (in Chinese).
徐程,庞亮. 2010. 直线三极法接地电阻测试原理与方法[J]. 科技创业月刊,(10):177–179. doi: 10.3969/j.issn.1672-2272.2010.10.080 Xu C,Pang L. 2010. Principle and method of grounding resistance measurement by three pole method[J]. Pioneering With Science &Technology Monthly,(10):177–179 (in Chinese).
杨德荣,梁丹. 2007. 接地电阻测量中应注意的几个问题[J]. 气象研究与应用,28(2):72–74. doi: 10.3969/j.issn.1673-8411.2007.02.019 Yang D R,Liang D. 2007. Several problems paying attention to in measuring grounding resistance[J]. Journal of Meteorological Research and Application,28(2):72–74 (in Chinese).
应顺潮. 1992. 接地电阻测量中测试极位置的确定[J]. 高电压技术,(1):80–83. Ying S C. 1992. Spacing of test electrodes for measurement of ground resistance[J]. High Voltage Engineering,(1):80–83 (in Chinese).
曾嵘,何金良,高延庆,孙为民,屠幼萍. 2001. 电极位置对垂直三层土壤结构中接地系统接地电阻测量结果的影响[J]. 清华大学学报(自然科学版),41(3):28–31. Zeng R,He J L,Gao Y Q,Sun W M,Tu Y P. 2001. Effects of electrode placement on grounding resistance of grounding system measured in vertical three-layer soil[J]. Journal of Tsinghua University (Science and Technology)
,41(3):28–31 (in Chinese). 郑兆苾,吴培稚. 1986. 地电阻率法的一个重要干扰因素:线间漏电[J]. 地震,(3):26–30. Zheng Z B,Wu P Z. 1986. An important interference factor of geo-electric resistivity method:Leakage between lines[J]. Earthquake,(3):26–30 (in Chinese).
中国地震局. 2009. DB/T 33.1—2009 地震地电观测方法地电阻率观测第1部分: 单极距观测[S]. 北京: 中国标准出版社: 3. China Earthquake Administration. 2009. DB/T 33.1−2009 The Method of Earthquake-Related Geoelectrical Monitoring: Geoelectrical Resistivity Observation: Part 1: Single Separation Configuration[S]. Beijing: Standards Press of China: 3 （in Chinese）.
中华人民共和国国家质量监督检验检疫总局. 2000. GB/T 17949.1—2000 接地系统的土壤电阻率、接地阻抗和地面电位测量导则第1部分: 常规测量[S]. 北京: 中国标准出版社: 20–24. General Administration of Quality Supervision, Inspection and Quarantine of the People’s Republic of China. 2000. GB/T 17949.1−2000 Guide for Measuring Earth Resistivity, Ground Impedance and Earth Surface Potentials of a Ground System: Part 1: Normal Measurements[S]. Beijing: Standards Press of China: 20–24 （in Chinese）.
IEEE Power and Energy Society. 2012. IEEE Std. 81−2012 Guide for Measuring Earth Resistivity, Ground Impedance, and Earth Surface Potentials of a Grounding System[S]. New York: Institute of Electrical and Electronics Engineers, Inc.: 22
Lu J,Xie T,Li M,Wang Y L,Ren Y X,Gao S D,Wang L W,Zhao J L. 2016. Monitoring shallow resistivity changes prior to the 12 May 2008 M8.0 Wenchuan earthquake on the Longmen Shan tectonic zone,China[J]. Tectonophysics,675:244–257. doi: 10.1016/j.tecto.2016.03.006
-
期刊类型引用(7)
1. 孙一男，王亮，程国亮. 华北克拉通中西部一维地壳结构特征. 华北地震科学. 2022(03): 40-49 . 
百度学术
2. 曾宪伟，冯建刚，龙锋，莘海亮. 鄂尔多斯西缘中上地壳Pg波速度成像研究. 地震研究. 2017(02): 176-185+333 . 百度学术
3. 魏建民，韩晓明，张帆，陈立峰，李娟，杨红樱. 2015年阿拉善左旗M_S5.8地震序列特征及发震破裂面讨论. 地震工程学报. 2017(05): 919-924 . 百度学术
4. 谢辉，马禾青，马小军，李青梅，张楠，任家琪. 基于背景噪声的宁夏及邻区瑞利波层析成像. 地震. 2016(02): 26-37 . 百度学术
5. 金春华，田小惠，蔡新华，何秋菊. MSDP软件定位方法应用对比. 防灾减灾学报. 2015(01): 65-70 . 百度学术
6. 金春华，何秋菊，蔡新华，田小惠. 宁夏地区地壳新速度模型在地震定位中的应用. 地震. 2015(02): 51-60 . 百度学术
7. 李启雷，李玉丽. 利用地震走时计算波速与震源深度. 地震研究. 2014(S1): 94-97 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 837
- HTML全文浏览量: 212
- PDF下载量: 109
- 被引次数: 9
