海洋电磁地层各向异性电阻率与发射源姿态联合反演方法
阅读说明：本技术 海洋电磁地层各向异性电阻率与发射源姿态联合反演方法 (Ocean electromagnetic formation anisotropic resistivity and emission source posture joint inversion method ) 是由 罗鸣 李予国 段双敏 李盼 吉芙蓉 吴俊良 于 2020-08-05 设计创作，主要内容包括：本发明提供一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法,包括：读取并转换参与反演的海洋可控源电磁场数据；设置联合反演执行参数；建立观测系统并设计联合反演初始模型；构建各向异性电阻率和发射源位置参数联合反演目标函数；求取电磁场关于各向异性电阻率和发射源位置参数雅各比矩阵和海森矩阵；基于反演参数特性自适应计算正则化因子；求取模型更新量；计算反演迭代模型的目标函数拟合差；判断是否满足反演要求；输出最终反演模型。该方法同时考虑海底介质电阻率各向异性和发射源的位置及姿态参数对反演结果的影响,能够有效解决发射源位置及姿态参数不准确给海洋可控源电磁资料解释带来的影响,适用性更广,容错性更强。(The invention provides a joint inversion method of anisotropic resistivity and emission source attitude of a marine electromagnetic stratum, which comprises the following steps: reading and converting ocean controllable source electromagnetic field data participating in inversion; setting a joint inversion execution parameter; establishing an observation system and designing a joint inversion initial model; constructing an anisotropic resistivity and emission source position parameter joint inversion target function; solving a jacobian matrix and a hessian matrix of parameters of the electromagnetic field, which are related to the anisotropic resistivity and the position of the transmitting source; adaptively calculating a regularization factor based on the inversion parameter characteristics; calculating the model updating quantity; calculating a target function fitting difference of the inversion iteration model; judging whether inversion requirements are met; and outputting the final inversion model. The method considers the influences of the resistivity anisotropy of the seabed medium and the position and attitude parameters of the emission source on the inversion result, can effectively solve the influence of the inaccurate position and attitude parameters of the emission source on the interpretation of the electromagnetic data of the marine controllable source, and has wider applicability and stronger fault tolerance.)
技术领域
本发明涉及多参数联合反演的海洋地球物理技术领域，具体涉及一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法。
背景技术
海洋可控源电磁法(CSEM)是探测海底油气资源和矿产资源的一种海洋地球物理勘探方法。随着资源开发和深地探测的需求不断提高，海洋电磁方法在矿产资源的调查、海洋地壳与地幔结构研究中正发挥越来越大的作用。频率域海洋CSEM方法通常使用拖曳在离海底上方几十米处的水平电偶极源作为发射源，并在被科考船拖曳过程中向电磁采集站发射范围为0.1-10Hz低频电磁信号，通过对布设在海底面采集站接收电磁信号进行反演获得海底介质的电性分布情况。电磁采集站接收电磁信号的强弱和质量高低不仅与海底介质的电阻率有关，同时与发射源和接收站的相对位置和姿态参数精度有直接的关系。在实际海上作业中发现，考察船在拖曳发射源的过程中发射源受海水运动和船速不稳定等因素影响，使得发射源不能按照预设的路径前进，其位置和姿态参数会出现较大范围的浮动，而这种变化量很难被监测和记录下来，致使发射源记录的位置和姿态参数与实际出现偏差，不准确的发射源位置和姿态数据对海洋电磁资料的处理和反演解释带来了极大的困难。因此，需对发射源的位置和姿态参数进行校准以提高海洋电磁信数据的质量。
研究表明，世界上大约30％的油气资源赋存于电性各向异性地层中，并且海底介质的电性各向异性是获得海底介质正确的电性分布的重要影响因素，若在解释海洋电磁资料时忽略电阻率各向异性，则很有可能无法获得合理的海底地电模型。然而，目前应用广泛的反演解释方法仍为传统的地球物理反演方法，该类方法在反演海洋资源勘探获取的海洋电磁资料时，常常假定海底介质为电性各向同性，这样的反演方法可能给后期资料解释提供巨大误差的反演结果，这也可能将导致获得错误的解释结果。为此，解释海洋电磁资料时必须考虑海底介质的电阻率各向异性的性质。
在现有的反演方法中，鲜有考虑电磁发射源位置和姿态参数对反演结果的影响，大多方法均是直接利用存在误差的电磁发射源参数进行反演，其反演结果必定也存在误差。因此，若能提出一种有效改正电磁发射源位置和姿态参数的反演解释方法将有效提高电磁数据解释的准确性。
发明内容
本发明的目的在于提供一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，该方法可应用于海洋电磁探测和反演解释的精确性。
为了实现以上目的，本发明提供如下技术方案：
1.一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，其特征在于，主要包括：
S1、读取并转换参与反演的海洋可控源电磁场数据，所述数据包括电磁场实部与虚部数据、振幅与相位数据、极化椭圆长轴与短轴参数；
S2、设置联合反演执行参数，所述参数包括反演最大迭代次数、目标拟合差、最大迭代步长、步长比例系数、惩罚函数类型、正则化衰减系数；
S3、建立观测系统并设计联合反演初始模型；
S4、构建各向异性电阻率和发射源姿态参数联合反演目标函数；
S5、求取电磁场关于各向异性电阻率和发射源姿态参数雅各比矩阵和海森矩阵；
S6、基于反演参数特性自适应计算正则化因子；
S7、求取各向异性电阻率和发射源姿态参数更新量；
S8、计算反演迭代模型的目标函数拟合差；
S9、判断是否满足反演要求，若满足则转向S10，不满足则转向S5；
S10、输出最终反演模型。
2.如权利要求1所述的海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，其特征在于，构建各向异性电阻率和发射源姿态参数联合反演目标函数为：
式中，φ为反演算法的目标函数；m为模型反演参数向量，其包括海底各向异性电阻率参数mρ、发射源的位置参数mP和发射源的姿态参数mT，即m＝mρ+mP+mT；为模型参数向量的梯度；||·||为标准差算子；d为反演使用的观测数据向量；Wd为数据加权矩阵；Wm为模型加权矩阵；F(m)表示模型m的正演响应算子；μρ、μp和μT分别为反演模型中海底各向异性电阻率参数mρ、发射源位置参数mP和发射源位置参数mT的正则化因子。
3.如权利要求1所述的海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，其特征在于，所述自适应正则化因子按照如下公式确定：
式中，i表示第i次反演迭代；μi为正则化因子；Max|·|为求取矩阵绝对值最大的元素；amj为矩阵乘积[(WdJ)T(WdJ)]的元素；M为矩阵乘积[(WdJ)T(WdJ)]的维度；χ为衰减系数；λ为反演模型海底地层的横向电阻率ρh、垂向电阻率ρv、发射源位置参数(x，y)和发射源姿态参数(Azm，Dip)的加权因子，采用以下公式确定
其中，α、β为权重系数，m为反演模型参数，(x，y)为发射源的横坐标和纵坐标，(Azm，Dip)为发射源的方位角和倾角。
4.如权利要求1所述的海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，其特征在于，基于反演参数计算模型参数更新量的方法为：
其中，Δm为下一次迭代的模型参数更新量，i为第i次反演迭代，Hi为海森矩阵，gi为目标函数的梯度。
5.如权利要求4所述的海底各向异性电阻率与发射源位置联合反演方法，其特征在于，求取目标函数的极小值的方法为：
其中，i为第i次反演迭代，Ji为雅各比矩阵。
6.如权利要求1所述的海底各向异性电阻率与发射源位置联合反演方法，其特征在于，求取雅各比矩阵Ji的方法为：
其中，i为第i次反演迭代；Ji为正演响应F(m)的雅各比矩阵；ρ＝(ρh，ρv)为地层的各向异性电阻率分布，P＝(x，y)为发射源位置参数，T＝(Azm，Dip)为发射源姿态参数。
7.如权利要求1所述的海底各向异性电阻率与发射源位置联合反演方法，其特征在于，求取海森矩阵Hi的方法为：
其中，i表示第i次反演迭代次数。
较现有技术相比，本发明一些实施例中，提供的方法的有益效果在于：
本发明主要针对海洋可控源探测中广泛存在的海底介质各向异性问题，以及海洋可控源电磁发射源记录的位置和姿态参数与实际位置和姿态参数存在偏差导致数据质量不高的问题，提出了一种海洋电磁地层各向异性电阻率与发射源位置和姿态参数联合反演方法，该方法在反演过程中考虑海底介质电阻率各向异性的同时，也考虑发射源的位置和姿态参数对反演结果的影响，能够同时对海底介质各向异性电阻率，以及电磁发射源的位置参数进行反演，最终可以获得与实际数据拟合最佳的海底各向异性电阻率分布情况，以及电磁发射源的准确的位置和姿态参数信息。该反演方法不仅为复杂的海底介质电阻率各向异性电磁资料的解释提供了一种有效的反演方法，同时也为电磁发射源位置和姿态参数存在误差情况下的数据处理问题提供了一种可行的技术手段。较于传统的反演算法，该方法同时考虑海底介质电阻率各向异性和发射源的位置及姿态参数对反演结果的影响，联合二者作为反演参数同时进行反演，该方法能够有效解决发射源位置及姿态参数不准确给海洋可控源电磁资料解释带来的影响；所提出的联合反演方法适用性更广，能够适用于电阻率各向异性资料的反演解释；容错性更强，能够处理电磁发射源位置和姿态参数存在一定误差时的反演问题。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的部分实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明方法的程序流程框图；
图2为一维电阻率各向异性模型结构示意图；
图3为反演迭代模型目标拟合差与反演迭代次数的变化图；
图4为发射源位置参数(x，y)的反演结果示意图；
图5为发射源姿态参数(Azm，Dip)的反演结果示意图；
图6为海底各向异性电阻率的反演示意图。
具体实施方式
为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本发明提供一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，该方法同时考虑海底介质电阻率各向异性和发射源的位置参数对反演结果的影响，联合二者作为反演参数同时进行反演，该方法能够有效解决发射源位置和姿态参数不准确引起海洋可控源电磁资料质量不高的问题。
参见图1，一种海洋电磁地层各向异性电阻率与发射源姿态联合反演方法，主要包括以下步骤：
S1.读取并转换参与反演的海洋可控源电磁场数据，所述数据包括电磁场实部与虚部数据、振幅与相位数据、极化椭圆长轴与短轴参数。
根据用户需求，将参与联合反演的海洋可控源电磁数据文件转换成特定格式的数据文件，并将输入的电磁场数据根据用户要求转换成特定参量(电磁场实部与虚部、振幅与相位、极化椭圆长轴与短轴参数等)，完成反演数据的读入并启动联合反演。
S2.设置联合反演执行参数，所述参数包括反演最大迭代次数、目标拟合差、最大迭代步长、步长比例系数、惩罚函数类型、正则化衰减系数。
联合反演过程是一个复杂且庞大的运算体系，在实现反演迭代直至获得最终的联合反演模型的过程中，需要运用到大量参数。不同的反演执行参数能够获得不同的反演效果，用户可根据参与联合反演数据的特点以及需要用于的特定需求设置反演执行参数。主要涉及的反演参数包括：反演最大迭代次数、目标拟合差、最大迭代步长、步长比例系数、惩罚函数类型、参与反演的数据分量、正则化衰减系数等。
S21.反演最大迭代次数：反演迭代次数的最大值。设置该参数能够使联合反演过程达到该最大迭代次数后退出联合反演程序，避免反演程序持续循环运行，提高程序计算效率。
S22.目标拟合差：反演模型与真实模型拟合程度的目标值。联合反演的最终目标是让反演模型不断逼近于真实模型，数据拟合差即是判断反演模型与真实模型的拟合标准。用户根据数据质量设定目标拟合差，当反演模型的拟合差达到或者低于目标拟合差，联合反演退出；若反演模型的拟合差没有达到目标拟合差，联合反演程序将持续进行。
S23.最大迭代步长：用户设置联合反演过程中模型参数变化量的最大值。在联合反演中，由于海洋可控源电磁场(或其他参量)的数量级跨度较大，所以通常都是取对数后再参与联合反演运算。若迭代的模型变化量过大，其结果变化非常大，例如模型参量(取对数后)变化量为1.0，即其真实值即提高了一个数量级，这样的变化值通常情况是不合理的，因此设计该参数限制模型参数的变化量大小。
S24.步长比例系数：用户设置联合反演过程中模型参数变换量的比例系数。联合反演涉及的反演参数较多，在反演迭代过程中，不同反演参数的变化程度不同，为使得整体的模型变化量在一个较为合理的范围，可将模型参数变化量乘以一个比率系数(即步长比例系数)，以保证反演模型稳健收敛。
S25.惩罚函数类型：模型参量之间的制约方式。在联合反演过程中，所有的模型参量同时参与反演，不同反演参量之间存在一定的联系，通过选择不同的惩罚函数可实现模型参量之间的制约方式，可使得模型更加光滑或曲度增大等效果。
S26.正则化衰减系数：自适应选择正则化因子的衰减系数。在反演迭代过程中，为保证不同参数之间的比率关系，正则化因子需随反演次数的改变而做相应改变，该参数即是设置正则化因子随迭代次数的衰减系数。
S3.设置联合反演初始模型参数，所述参数包括背景层电阻率参数、厚度参数、观测系统参数。
反演初始模型是联合反演模型参数的初始值，通常情况下需根据输入数据的观测系统和已知的水深数据等资料设定。
S4.构建各向异性电阻率和发射源位置、姿态参数联合反演目标函数。
海洋可控源电磁海底各向异性电阻率与电磁发射源位置、姿态参数联合反演涉及多属性的反演参数，参数间关系复杂。鉴于海底介质不同电阻率各向异性参数，以及发射源位置和姿态参数之间差异性，本发明利用正则化约束稳定海洋可控源电磁海底各向异性电阻率与电磁发射源位置参数联合反演过程。所采用的目标函数为：
式中，φ为反演算法的目标函数；m为模型反演参数向量，其包括海底各向异性电阻率参数mρ、发射源的位置参数mP和发射源的姿态参数mT，即m=mρ+mP+mT；为模型参数向量的梯度；||·||为标准差算子；d为反演使用的观测数据向量；Wd为数据加权矩阵；Wm为模型加权矩阵；F(m)表示模型m的正演响应算子；μρ、μp和μT分别为反演模型中海底各向异性电阻率参数mρ、发射源位置参数mP和发射源位置参数mT的正则化因子。
反演模型的海底电阻率反演参数mρ有以下形式：
式中，ρh和ρv分别为海底地层的横向电阻率和垂向电阻率，M为海底地层的层数。
发射源位置参数mp有以下形式：
mP＝[l0g10x1…log10xt，log10y1…log10yt]T；
倾斜电偶极源的姿态反演参数mT有以下形式：
mT＝[Azm1…Azmt，Dip1…Dipt]T；
式中，(x，y)为发射源位置参数(横坐标，纵坐标)，(Azm，Dip)为发射源姿态参数(方位角，倾角)，t为电偶极源的个数。
S5.求取电磁场关于各向异性电阻率和发射源姿态参数雅各比矩阵和海森矩阵。
S51.雅各比矩阵Ji是电磁场关于反演参数的偏导数矩阵，不同的反演参数有不同的形式：
其中，i为第i次反演迭代；Ji为正演响应F(m)的雅各比矩阵；ρ＝(ρh，ρv)为地层的各向异性电阻率分布，P＝(x，y)为发射源位置参数，T＝(Azm，Dip)为发射源姿态参数。
若将电阻率张量ρ为反演模型横向电阻率ρh和垂向电阻率ρv的函数，即ρ＝f(ρh，ρv)，因此，电磁场关于电阻率张量ρ的雅克比矩阵有如下形式：
同理的，可得到电磁场关于倾斜电偶极源的位置反演参数mp和姿态参数mT的雅各比矩阵为：
式中，(x，y)为发射源位置参数(横坐标，纵坐标)，(Azm，Dip)为发射源姿态参数(方位角，倾角)，t为电偶极源的个数。
S52.海森矩阵Hi为目标函数对反演参数的二阶导数，忽略二阶导数项和非对称项，海森矩阵Hi可简化为：
其中，i为第i次反演迭代，Ji为正演响应F(m)的雅各比矩阵；mρ、mP和mT分别为海底电阻率参数、发射源的位置参数和姿态参数；
S6.基于反演参数特性自适应计算正则化因子。
正则化因子μ的选择方式是求解海洋可控源电磁海底各向异性电阻率与电磁发射源参数联合反演问题合理与否的关键。在地球物理反演方法中，出现了许多正则化因子选择方法。在进行海洋可控源电磁海底各向异性电阻率与电磁发射源参数联合反演时，本发明针对不同的反演参数选择不同的正则化因子。首先，选择反演过程中的关于雅各比矩阵的特征参数作为基数，然后，再结合反演参数的特征及参数之间的关系调节不同反演参数的正则化因子大小，从而实现正则化因子的自适应选择。
正则化因子μi可写成如下形式：
式中，i表示第i次反演迭代；μi为正则化因子；Max|·|为求取矩阵绝对值最大的元素；amj为矩阵乘积[(WdJ)T(WdJ)]的元素；M为矩阵乘积[(WdJ)T(WdJ)]的维度；χ为衰减系数；λ为反演模型横向电阻率ρh、垂向电阻率ρv、发射源位置参数(x，y)和发射源姿态参数(Azm，Dip)的加权因子，采用以下公式确定：
其中，α、β为权重系数，m为反演模型参数，权重向量的初始值λ1为单位向量。当反演模型地层中存在电阻率各向异性时，正则化因子可根据电阻率各向异性率自适应地进行调节，从而调节目标函数中各元素之间的权重，从而达到自适应调节反演过程的效果。
S7.求取模型更新量。
通过求取目标函数的极小值可计算模型更新量，即令目标函数梯度gi＝0，可由下式确定
其中，i为第i次反演迭代，Ji为正演响应F(m)的雅各比矩阵；φ为反演算法的目标函数；mρ、mP和mT分别为海底电阻率参数、发射源的位置参数和姿态参数；；为模型参数向量的梯度；||·||为标准差算子；d为反演使用的观测数据向量；Wd为数据加权矩阵；Wm为模型加权矩阵；F(m)表示模型m的正演响应算子；μρ、μp和μT分别为反演模型中海底各向异性电阻率参数mρ、发射源位置参数mP和发射源位置参数mT的正则化因子。
计算得到下一次迭代的模型参量更新量为Δm，可表示为：
其中，Δm为下一次迭代的模型参数更新量，i为第i次反演迭代，Hi为海森矩阵，gi为目标函数的梯度。
S8.计算反演迭代模型的目标函数拟合差。
在反演迭代过程中，反演迭代模型逐步收敛于真实模型，其目标函数拟合差ψ逐渐收敛于1.0。目标函数拟合差ψ有如下形式：
其中，i为反演迭代次数，d为反演使用的观测数据向量；F(m)表示模型m的正演响应算子；N为反演数据的个数；k为反演数据的编号；δk为第k个数据的标准差。
S9.判断是否满足反演要求，若满足则转向S9，若不满足则转向S5。
判断是否退出联合反演的标准有两个：1)是否满足联合反演目标拟合差；2)反演迭代次数是否达到最大迭代次数。若否，则转向S5，继续迭代求解模型更新量；若是，则转向S10。
S10.输出最终反演模型。
参考图2，为一维典型地电模型原理图。为了验证海洋可控源电磁海底各向异性电阻率与电磁发射源参数联合反演的有效性，以图2所示一维电阻率各向异性模型为例，利用合成数据进行海底各向异性电阻率与发射源参数的联合反演。假设26个电磁发射源等间距地布置于测线500m-13000m范围内，且均位于海底正上方50m处，每个发射源的姿态都相同且方位角都为90度，倾角5度；1个接收站布设于海底面(0，0，1000)的位置。设定发射电流为1安培。反演数据由电磁场各分量的实部和虚部组成，各数据加入了2％的随机高斯噪声。在反演算例中，反演初始模型为空气、海水和电阻率为1Ωm的均匀半空间，空气层和海水电阻率和深度固定。设定反演初始模型中电磁发射源测线的位置与真实测线位置交叉(存在一个约45°的夹角)。
参考图3，为反演迭代模型目标拟合差与反演迭代次数的关系示意图。由图可见，在联合反演过程中，目标函数拟合差持续减小并逐渐向真实模式靠近，并由最初较大的拟合差17，经过20次反演迭代后，最终收敛于3.1。
参考图4，为发射源位置参数反演结果。由图4可见，图中纵坐标和横坐标分别为(x，y)的位置，空心三角为电磁发射源的初始位置，实心三角为电磁发射源真实位置，实心圆点为联合反演获得的电磁发射源位置，即使输入的发射源位置与真实位置差异较大(位置差异达100m)，联合反演得仍然能够得到相对准确的位置信息(反演得到的电磁发射源位置与真实位置距离不超过25m)，由此说明，联合反演得到的发射源位置是准确的。
参考图5，为发射源姿态参数反演结果。横坐标为发射源的y坐标，左纵坐标为电磁发射源的方位角，右纵坐标为电磁发射源的倾角，方块和六角形空心线分别为发射源方位角和倾角的初始值，方块和六角形实心线为联合反演得到的电磁发射源方位角和倾角，可见联合反演的发射源方位角和倾角与真实值的偏差分别不超过5度和1度，由此说明，反演得到的电磁发射源的姿态参数(方位角和倾角)均是准确的；
参考图6，为海底各向异性电阻率的反演结果图，浅灰色和深灰色线段为真实的海底介质横向和垂向电阻率分布情况，浅灰色和深灰色阶梯曲线为反演得到的海底介质横向电阻率和垂向电阻率曲线，由图可见，联合反演的围岩的各向异性电阻率与真实情况相符，高阻层的埋深、厚度和垂向电阻率值也得到了较准确的恢复。综上所述，本发明提出的海洋可控源电磁海底各向异性电阻率与电磁发射源参数联合反演方法能够反演出海底介质的各向异性电阻率信息，以及电磁发射源的位置态参数。由此也说明，本发明提出的算法是有效的。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。