CN107710023B - 反演问题的有效解 - Google Patents

Abstract

一种求解地球物理学反演问题以用于估计物理参数的方法，所述方法包括：提供代表所述物理参数的模型向量，由第一算子变换所述模型向量以提供第一变换模型向量，针对所述第一变换模型向量求解所述反演问题以提供第一解，由第二算子变换所述模型向量以产生第二变换模型向量，针对所述第二变换模型向量求解所述反演问题以提供第二解，计算所述第一解和所述第二解的加权和。

Description

反演问题的有效解

技术领域

本发明涉及物理数据的线性和非线性反演，并且更具体但不限于受控源电磁(CSEM)数据或声学数据的线性和非线性反演。

背景技术

已经开发了若干种用于探索地球地表下的技术，其基于向地球地表下的区域中发射波或信号。发射的信号与地球相互作用，并且通常信号的一部分传播回到表面，该信号在此被记录并用于基于信号如何与地球相互作用而获得关于地表下结构的信息。CSEM方法例如使用在海底上方牵引的偶极子源来发射电磁信号，并且使用放在海床上的接收器阵列来检测已经行进通过海底下方的地层的信号。然后需要对检测到的信号进行反演以推导物理参数。物理参数可以任选地用于估计烃或水的存在。可以推导的物理参数的示例是地层的电导率。电导率可以在模拟中用作参数，由此模拟能够模拟所记录的数据。电导率的最优值是优化所记录的数据的模拟和数据之间的一致性的那些值。

CSEM数据的非线性反演涉及求解大型线性方程组，以在迭代优化方法的每次迭代时计算电导率的更新，以便使数据和模拟的数据之间的距离最小化。模型向量所依据的空间三维网格的节点的数量通常超过百万个，并且在该网格上求解标准方程变得不可行。可以使用优化算法，如有限存储器Broyden-Fletcher-Goldfarb-Shanno算法，但是这些算法需要非常准确的起始模型向量才能提供好的反演结果。减少反演参数(也称为自由参数)的数量是提高反演算法的效率的重要方式。

可以使用强几何约束来减少独立自由参数的数量。例如，在地质上定义的主体中的电导率值可以被设置为恒定值，并且然后对那些恒定值进行反演。然而，这种方法的缺点是它们需要很多先验信息，并且在电阻或导电体跨越所限定的几何结构伸展时，这种方法是不适合的。

发明内容

根据本发明的第一方面，提供了一种求解地球物理学反演问题以用于估计物理参数的方法，所述方法包括：提供代表所述物理参数的模型向量，由第一算子变换所述模型向量以提供第一变换且压缩的模型向量，针对所述第一变换且压缩的模型向量求解反演问题以提供第一解，由第二算子变换所述模型向量以产生第二变换且压缩的模型向量，针对所述第二变换且压缩的模型向量求解反演问题以提供第二解，以及计算所述第一解和所述第二解的加权和。计算所述第一解和所述第二解的加权和可以提供所述物理参数的估计。

所述第一解和所述第二解的加权和可以使一组测量结果和所述一组测量结果的模拟之间的距离最小化，其中所述模拟取决于所述模型向量。由第一算子变换所述模型向量可以包括根据第一压缩因子减少所述反演问题中的自由参数的数量，并且其中，由第二算子变换所述模型向量包括根据第二压缩因子减少自由参数的数量。替代地，由第一算子变换所述模型向量可以包括向所述模型向量应用第一傅里叶变换，并且其中，由第二算子变换所述模型向量包括向所述模型向量应用第二傅里叶变换，由此可以使用滤波器来选择傅里叶谱的一部分。

减少反演问题中的自由参数的数量可以包括：对在细网格上定义的模型向量进行插值以提供在第一粗网格上定义的第一变换模型向量，其中所述第一粗网格的分辨率低于所述细网格的分辨率；以及对在细网格上定义的所述模型向量进行插值以提供在第二粗网格上定义的第二变换模型向量，其中所述第二粗网格的分辨率低于所述细网格的分辨率。所述粗网格可以是均匀或非均匀的，并且所述粗网格可以取决于地层的物理结构的先验知识。所述第一粗网格相对于所述第二粗网格偏移。

所述一组测量结果可以是一组地震、电、磁或重量测定测量结果之一。模型向量包括迭代反演方法中的模型更新，并且可以利用线搜索来估计最小距离。所述反演问题可以为线性或非线性反演问题。

根据本发明的第二方面，提供了一种计算机系统，其被布置成执行根据本发明的第一方面的方法。

根据本发明的第三方面，提供了一种计算机软件，其在被安装在根据本发明的第二方面的计算机系统上时，被布置成使所述计算机系统执行根据本发明的第一方面的方法。

附图说明

现在将仅通过示例的方式并参考附图来描述本发明的一些实施例，在附图中：

图1是示出了本文公开的方法的流程图，

图2示意性示出了在该方法中使用的数据集，

图3示意性示出了该方法的部分解，

图4示意性示出了该方法的两个部分解，

图5示意性示出了在该方法中使用的解的组合。

具体实施方式

如图1所示，本方法通过如下方式求解地球物理学反演问题以用于估计物理参数：提供代表物理参数的模型向量(S1)，通过第一算子变换所述模型向量以提供第一变换模型向量(S2)，针对第一变换模型向量求解反演问题以提供第一解(S3)。该模型然后通过第二算子变换所述模型向量以产生第二变换模型向量(S4)，并针对第二变换模型向量求解反演问题以提供第二解(S5)。该方法最后计算第一解和第二解的加权和(S6)，这提供了对物理参数的估计。尽管给出了两个解的示例，但该方法在大部分应用中使用很多解的组合。

所述变换可以是傅里叶变换，由此将滤波器应用于变换模型向量以过滤掉一些频率分量。替代地，所述变换可以是对模型向量中的自由参数的数量的减少，以减小问题的大小。所述变换也可以是主成分分析，其中排除了一些较不重要的成分。

本方法的示例涉及求解地球物理学反演问题，其中直接计算模型参数的N维向量(线性反演)或使用优化策略来迭代地更新模型参数的N维向量(非线性反演)。该方法包括在与在其上定义模型参数的向量的网格相比具有减小的维度的网格上减小模型参数的向量。求解减小的模型向量没有求解完整的模型向量那么麻烦。使用相对于第一网格偏移的第二网格重复减小的步骤，并且在第二网格上再次求解该问题。将第一解与第二解组合以提供组合解，该组合解的分辨率比第一解或第二解更高。

该网格可以是单元的阵列，由此在每个单元内，模型向量或解是恒定的。例如，所述单元在真实空间中可以是地层的三维区，并且在那些单元内，可以将诸如电阻率的物理值选择为恒定。

在大部分地球物理学反演问题中，必须要在反演过程中对以下形式的一个或多个方程组求解一次或几次(方程1)：

该方程中的量具有以下含义：

符号m是代表要估计的未知参数的向量。这些参数常常被称为模型参数和m模型向量，因为在很多反演问题中，估计这些参数允许构建地球地表下的某些性质(例如，机械、电气或磁性性质)的3D模型。这种性质的示例是电阻率、压缩或切变波的速度、声阻抗、密度和磁化率。替代地，m代表初始模型向量或方程组(1)的解的先前近似值的可能更新。m的大小为(N,1)，其中N是参数的数量。

符号d代表包含测量数据或包含测量数据和模拟数据之间的差的列向量。这些数据可以具有不同类型：地震、电、磁、重量测定等等。d的大小为(N_d,1)，其中N_d是数据测量结果的数量。C_d通常是数据不确定性的协方差矩阵或其近似值，无论是否为对角矩阵。

是C_d的逆矩阵或其近似值，无论是否为对角矩阵。

符号F是大小(N_d,N)的线性建模算子。在线性反演问题的情况下，F代表正向算子(也称为建模算子)，以使F乘以给定模型向量代表对应于该特定模型向量的模拟数据。替代地，在非线性反演问题的情况下，F代表非线性正向/建模算子F的线性近似。在该情况下，F常常被称为Jacobian。

星号*代表共轭转置，或者在仅涉及实量时，星号*简单地代表转置。

符号Q和Q_P是为了使数值反演问题稳定或使其正则化(即，为了在很多不同模型向量能够给出类似模拟数据的情况下支持特定类型的解)而引入的矩阵。Q的典型选择是Q＝I(例如在Levenberg Marquardt算法中)，其中I是单位矩阵，Q＝S^*S，其中S是一阶或二阶微分算子(例如Occam算法和特定类型的Tickhonov正则化)。

常常等于Q但也可以采用其它值。

符号λ是引导稳定化或正则化的量的标量参数。λ_g也是标量参数，通常等于λ，但也可以采取不同的值，例如零值。这些参数可以是固定的，或者可以在运行若干次迭代时在迭代过程期间发生改变，或者也可以通过优化过程(例如，Levenberg-Marquardt和OCCAM算法)来估计这些参数。

在使用单次迭代估计p时，符号Δm通常为零。在使用若干次迭代时，它典型地代表在前一次迭代中估计的参数向量，或者代表在前一次迭代中估计的参数向量和先验参数向量之间的差。

在迭代算法中，可以通过将通过求解方程组(1)所获得的模型更新向量m加到该估计值来细化模型向量的前一估计值。替代地，可以将αm加到模型向量的该前一估计值，其中可以通过不同技术来优化标量参数α。不同技术的示例为所谓的线搜索技术，其中尝试找到α的最优值以进一步减小数据和模拟数据之间的距离。利用高斯-牛顿算法，例如，可以使用线搜索流程。线搜索需要一些模拟，并且除了更全局性的迭代过程之外，可以使用内部“线搜索迭代”。对于线性反演，线搜索不是必需的，因为α的最优值为1，并且不需要再对其进行搜索。对于非线性反演，m会被解释为“搜索方向”，αm会被解释为“模型更新”。已经选择了覆盖线性和非线性反演的标记。

在某些情况下，上文刚定义的量的部分由更容易计算的量来替代。例如，在迭代过程期间可以近似地建立矩阵对(1)的左手侧的影响的近似，而不需要计算矩阵或对其进行反演(例如，共轭梯度、BFGS和LBFGS算法)。

必须要求解类似于(1)的方程组的所有算法的一个共同点在于，当模型向量的尺寸显著增大时，实际上可能无法求解该方程组：涉及的矩阵可能太大，以致于无法存储在计算机存储器中，或者数值运算的次数可能变得过大。例如，在必须要反演大型3维地理数据集时可能是这种情况。

由于该原因，已经开发了模型压缩技术。大型模型向量m由大小为(N_c,1)的较小(或压缩)的向量m_c替代，其中N_c<N。模型压缩因子可以定义为k＝N/N_c。N维原始模型空间必须与N_c维压缩模型空间相关。这通常是由线性关系(方程2)实现的

m＝Rm_c (2)

如果模型向量中的自由变量的数量减少，算子R可以被解释为插值。换言之，模型m包含N个值，但N个值中仅有N_c个自由度。如果k充分大，则估计向量m_c的N_c个参数比估计向量m的N个参数更简单，因为原始方程组(1)中涉及的矩阵和向量被小得多的矩阵和向量替代。一旦通过求解类似于(1)但小得多且更容易求解的方程组而计算出m_c，就可以通过方程(2)获得m。替代地，算子R可以被解释为傅里叶变换。可以在傅里叶变换之后应用滤波器以过滤掉高频分量或低频分量，由此也减少了自由变量的数量。在该情况下，参数m_c可以被解释为傅里叶系数，并且可以进一步被解释为低频系数。

然而，通过压缩获得的参数向量m_c包含明显比通过求解原始方程组(1)而直接估计的参数向量更少的细节。模型输出是低维参数向量的插值版本。典型地，这看起来是方程组(1)的理想解的低分辨率版本。

本文公开的方法受益于存储器要求和模型压缩技术所允许的数值运算的大幅降低，但具有高得多的分辨率。该方法对线性反演(即使在执行单次迭代时)和非线性反演都有效。对于迭代非线性反演，本发明即使在执行少量迭代时也有效。

发明人认识到，上文所述的尝试计算向量m或模型更新αm的所有反演问题都可以被分成一组L个小得多的问题，其具有解m_i,i＝1,L或α_im_i,i＝1,L。这些解中的每一个代表原始问题的解m(或αm)的低分辨率版本R_im_i(或R_iα_im_i)。然而，可以通过将L个低分辨率解叠加而获得分辨率好得多的解：

或者替代地

其中w_i代表在需要时可以与1不同的权重因子。权重可以被预先定义或例如通过线搜索技术来求解。

在求解较小的反演问题并将解组合时，可以像针对完整问题那样构造N维向量。尽管对于很多真实的数据集不能对完整问题进行反演，但可以对减小的问题进行反演。最终的解具有比个体成分R_im_i更好的分辨率。

L个较小的压缩问题不需要具有完全一样的压缩因子，并且对于不同模型，粗的单元可以具有不同的形状和大小。算子R_i可以各不不同。网格可以随机变化或以系统方式变化。也可以设置网格的单元的形状，以使它们遵循已知的地质结构，例如地震层位。在使用傅里叶变换的情况下，可以选择傅里叶系数m_i，使得每次计算对应于频谱的略微不同的部分，由此解的组合覆盖相关频率的整个频谱。

如果在节点之间进行插值R，这些节点不必放置成矩形和规则图案。相邻节点之间的距离例如可以随着深度而变化，或者可以遵循一些已知结构，或者被约束到相同区域中的另一反演问题的解。例如，反演问题可以是CSEM反演问题，并且可以选择网格以遵循从先验估计已知的地震层位。

所述一组较小的压缩问题不需要在每次迭代时是相同的(如果使用迭代估计方法)，较小问题的数量和类型可以在迭代过程期间变化。例如，可以在迭代进行时使用较小的压缩因子和较少数量L的压缩问题，只要判断其足以求解特定问题。

附图示出了针对所述一组较小问题的不同可能选择。图示被强烈简化以使要点更容易理解。模型被定义在二维矩形网格上。算子R_i是在规则间隔的节点上定义的简单的最近的相邻插值器。这意味着，向量R_im_i在粗矩形单元内部具有恒定值，节点位于粗单元的中心。粗矩形单元中包含的细网格(代表m)的单元的数量给出了压缩因子k的良好近似值。对于所有L个压缩反演问题，除了网格的边缘之外，所有单元都具有相同的大小和形状。因子w_i和α_i应该等于1。本发明当然不限于这些特定简化情况。

该方法是利用简单得多的问题来例示的，而不是特定的地球物理学问题。图2示出了网格11，在网格11上利用异常(anomaly)12定义了地层。我们假设要估计的正确模型向量(例如电阻率)在整个网格11上具有相等的恒定值，除了在一个特定单元12中，在单元12中所述模型向量具有较高值。

图3示出了具有异常12的网格11、连同粗网格上的解。所述解被图示为矩形21，并且在矩形21内，所述解具有比跨粗网格其余部分的周围值更高的恒定值。异常12位于矩形21内，并且所述解已经正确识别了异常的位置，只是分辨率低得多。压缩因子在竖直方向上为3，并且在水平方向上为4，因此总共为12。针对粗网格示出了参考点22。

图4示出了图2的特征，并且另外示出了已经在水平方向和竖直方向上偏移的第二粗网格的第二解。偏移由参考点32例示，使用细网格的坐标，将参考点32相对于参考点22向右偏移一个单元并向上偏移一个单元。第二解31也向右并向上偏移。同样，异常已经被正确识别，但具有较低的分辨率。

图5示出了两个解之和。所述解在异常的外部具有恒定值42，只是在区域21和31中具有较高的值。交叠区域41具有区域21和31的组合值。该交叠区域41的大小小于区域21和31，因此分辨率增大，尽管分辨率仍然低于图1的网格。显然，许多解的组合将开始精确地对图1的异常进行近似，但没有完整模型的计算成本。

本文以示例的方式描述的方法涉及CSEM数据的非线性反演，但所主张的发明也可以应用于其它方法。例如，其它类型的地球物理学数据的非线性反演，如大地电磁数据、地震数据、在钻孔中采集的声学数据和地面穿透雷达数据的非线性反演。对于地震数据，本发明的应用是所谓的完整波形反演。该方法也可以应用于联合反演，其中同时使用若干类型的数据。

可以从地球物理学数据反演的一些常见地球性质包括声波速度、地层和流体密度、声阻抗、泊松比、地层压缩率、抗剪刚度、孔隙度和流体饱和度。

确定性反演方法基于地球模型的输出与观测的现场数据的比较以及连续更新地球模型参数以使函数最小化，该函数通常是模型输出与现场观测值之间的某种形式的差。使目标函数最小化的一组模型参数将产生与收集的现场地震数据相比最好的数值地震图。还在能够在迭代之间改变的网格上执行更新的步骤。还可以使用随机反演方法以使用诸如克里金插值之类的地球统计工具来生成如在储层流动模拟中使用的约束模型。与产生单组模型参数的确定性反演方法相反，随机方法生成全部服从模型约束的一套替代的地球模型参数。

尽管已经依据上文阐述的优选实施例描述了本发明，但应当理解，这些实施例仅仅是例示性的，并且权利要求不限于那些实施例。鉴于本公开，本领域的技术人员将能够做出修改和替换，它们被视为落在所附权利要求的范围内。本说明书中公开或例示的每个特征可以被并入本发明，无论是单独的或与本文公开或例示的任何其它特征进行任何适当组合。

Claims (13)

1.一种求解地球物理学反演问题以用于估计物理参数的方法，所述方法包括：

提供代表所述物理参数的模型向量，

由第一算子变换所述模型向量以提供第一变换且压缩的模型向量，

针对所述第一变换且压缩的模型向量求解所述反演问题以提供第一解，

由第二算子变换所述模型向量以产生第二变换且压缩的模型向量，

针对所述第二变换且压缩的模型向量求解所述反演问题以提供第二解，

计算所述第一解和所述第二解的加权和，

其中，由第一算子变换所述模型向量包括根据第一压缩因子减少所述反演问题中的自由参数的数量，并且其中，由第二算子变换所述模型向量包括根据第二压缩因子减少所述自由参数的数量，

其中，减少所述反演问题中的所述自由参数的数量包括：对在细网格上定义的所述模型向量进行插值以提供在第一粗网格上定义的第一变换模型向量，其中，所述第一粗网格的分辨率低于所述细网格的分辨率；以及对在细网格上定义的所述模型向量进行插值以提供在第二粗网格上定义的第二变换模型向量，其中所述第二粗网格的分辨率低于所述细网格的分辨率，

其中，所述第一粗网格相对于所述第二粗网格偏移。

2.根据权利要求1所述的方法，其中，计算所述第一解和所述第二解的所述加权和提供了所述物理参数的估计。

3.根据权利要求1所述的方法，其中，所述第一解和所述第二解的所述加权和使一组测量结果与所述一组测量结果的模拟之间的距离最小化，其中，所述模拟取决于所述模型向量。

4.根据权利要求1所述的方法，其中，由第一算子变换所述模型向量包括将第一傅里叶变换应用到所述模型向量，并且其中，由第二算子变换所述模型向量包括将第二傅里叶变换应用到所述模型向量。

5.根据权利要求4所述的方法，其中，由所述第一算子或所述第二算子变换所述模型向量还包括应用滤波器。

6.根据权利要求1所述的方法，其中，所述粗网格为均匀或非均匀的。

7.根据权利要求1所述的方法，其中，所述粗网格取决于地层物理结构的先验知识。

8.根据权利要求3所述的方法，其中，所述一组测量结果是一组地震、电、磁或重量测定测量结果之一。

9.根据权利要求1所述的方法，其中，所述模型向量包括迭代反演方法中的模型更新。

10.根据权利要求9所述的方法，其中，最小距离是利用线搜索来估计的。

11.根据权利要求1所述的方法，其中，所述反演问题为线性或非线性反演问题。

12.一种计算机系统，被布置成执行根据权利要求1所述的方法。

13.一种计算机系统，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1中所述方法。

Images