一种井中磁力异常估计有效磁化倾角的方法与流程

本方法属于勘探地球物理领域，尤其涉及一种井中磁力异常估计有效磁化倾角的方法。
背景技术：
：地下异常体的磁化方向以及大小不仅受地磁场的感应磁化的影响，还会受到剩磁的影响。在剩余磁化存在的情况下，异常体的磁化方向会偏离地磁场的方向，使得感应磁化矢量化，并可能扭曲磁异常的形态，从而导致产生的异常场变得复杂，磁法勘探得到的数据难于解释。因此，估计异常体的实际磁化方向具有重要的意义。在估计磁化方向的问题上，gerovska等提出基于磁总场模量异常和磁异常化极数据互相关估计的方法估计磁化方向，于是基于该方法的原理，提出利用井中一维数据进行剖面有效磁化方向估计，并得到与利用地面二维数据估计磁化方向的方法的不同点。技术实现要素：一种井中磁力异常估计有效磁化倾角的方法，包括以下步骤：s1、基于井中磁异常三分量包括垂直分量za、水平分量hax和水平分量hay数据计算得到磁异常模量ta；s2、对za、hax和hay三个分量进行傅里叶正变换分别得到相应的频谱；s3、设置初始有效磁化倾角度数；s4、基于s3设置的有效磁化倾角度数以及井与异常体的相对位置关系计算出垂直磁化方向转换因子以及水平磁化方向转换因子；s5、za频谱乘以垂直磁化方向转换因子，hax和hay分别乘以水平磁化方向转换因子；s6、将乘以转换因子的频谱分别进行傅里叶逆变换得到磁化方向转到垂直方向的za数据以及磁化方向转到水平方向的hax和hay数据；s7、将s6中转换后的za、hax和hay分别与s1中磁异常模量ta进行相关性计算，记录得到的相关性数值；s8、将有效磁化倾角增加n°，0°<n°<360°，循环s4至s8，当有效磁化倾角增加到指定的度数，停止循环；s9、基于记录的相关性数值绘制出za、hax和hay的三条相关性曲线，寻找za相关性曲线的最小值点以及hax和hay相关性曲线的最大值点，得到的最值点对应的度数即为估计的有效磁化倾角。进一步的，s1中磁异常模量ta的计算公式如下：其中，ta表示磁异常模量，za表示垂直分量，hax表示水平分量，hay表示水平分量。进一步的，s3中初始有效磁化倾角度数设置为x°，其中0°≤x°≤360°。进一步的，s4中频率域转换因子与井和异常体的相对位置有关，具体关系为：在建立的坐标系下，所研究剖面中，当井位于异常体正方向即在对应横坐标轴正方向时，磁化方向转换因子为：α2＝cosis2γ2＝sinis2α1＝cosis1γ1＝sinis1其中，w表示磁化方向转化因子，is2为转换后磁化倾角，is1为原磁化倾角；当井位于异常体负方向即横坐标轴负方向时，转换因子为：α2＝cosis2γ2＝sinis2α1＝cosis1γ1＝sinis1其中，is2为转换后磁化倾角，is1为原磁化倾角，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；当is2为90°时，转换因子为垂直磁化方向转换因子；当is2为0°时，转换因子为水平磁化方向转换因子。进一步的，s4中垂直磁化方向转换因子以及水平磁化方向转换因子的具体推导过程如下：首先，根据地质体三个方向的延伸度，把磁性体分为二度体和三度体，对于地面来说，数据采集面和线永远处在异常体的上半空间，即在所建立的空间直角坐标系下，数据采集在z＝0的平面进行，异常体中心在z＝300的平面，z轴竖直向下为正，对于井来说，假设异常体为二度体，令建立的坐标系x轴垂直于异常体走向，y轴平行于异常体走向，数据采集沿着井轴即z轴方向进行，由于井轴在异常体的正方向即正x方向和异常体的负方向即负x方向，因此，井中数据与井和异常体的相对位置有关，所以，在推导二度体的频率域转换因子时需要分两种情况，即在所建立的空间直角坐标系下，井在异常体的正方向和负方向两种情况；对于二度体，水平分量hay为零，其位场转换因子是由二维拉普拉斯方程狄里希莱问题推导而来：其中，x和z表示x轴和z轴坐标，uxx为函数u对x的二次求导，uzz为函数u对z的二次求导，为边界条件即测量数据；情况1:井位于异常体负方向：对于这类情况，二维拉普拉斯方程狄里希莱问题解为：其中，u为狄里希莱问题的解，δ表示边界上任意一点的z轴坐标，t(x，z)满足上式，t为磁异常，x和z表示x轴和z轴坐标，即：对t(0，z)与做关于z的一维傅里叶变换：得到磁异常频谱表达式为：其中，st为磁异常频谱；f为频率，i为虚数单位；利用微分定理对式(1)进行求微分：可以得到对于情况1：x方向导数因子为2πf，z方向导数因子为2πif；给出二度体泊松公式：za磁异常垂直分量，ha为磁异常水平分量，m为磁化强度的模，α，γ为磁化强度的方向余弦，δ为地质体密度，g为万有引力常量，v为引力位；vxz和vzx为对v的x方向和z方向的二次偏导数，vzz为对v的z方向的二次导数,vxx为对v的x方向的二次导数；令引力位频谱为sv，对引力位频谱进行x和z方向求导：对进行二阶求导：将上式代入(4)与式(5)可得：其中，sz为磁异常垂直分量频谱，sh为磁异常水平分量频谱；根据式(6)与式(7)可以得到，sz，sh相互转换的分量转换因子为：sz→sh：-i(8)sh→sz：i(9)二度体新磁化方向磁位即α2，γ2为新磁化方向余弦为：u2是新磁化方向的磁位，t1表示原磁化方向，即新磁位是沿着原磁化方向t1反方向的曲线积分令sh为ha的频谱，ha频谱表达式为：令sz为za的频谱，za频谱表达式为：利用式(8)、式(9)得到的分量转换因子可得：由式(10)式(11)式(13)可得：其中，α1，γ1为原磁化方向余弦；根据磁位与分量之间的关系可以推出新分量(ha2，za2)的频谱表达式：因此对于井位于异常体负方向:任意磁化方向转换转换因子为：其中，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；情况2:井位于异常体正方向:对于这类情况，二维拉普拉斯方程狄里希莱问题解为：t(x，z)满足上式，t为磁异常，即：对t(0，z)与做关于z的一维傅里叶变换：得到磁场频谱表达式为：st为磁异常频谱，f为频率，i为虚数单位；利用微分定理对式(17)进行求微分：可以得到对于情况2：x方向导数因子为-2πf，z方向导数因子为2πif给出二度体泊松公式：za磁异常垂直分量，ha为磁异常水平分量，m为磁化强度的模，α，γ为磁化强度的方向余弦，δ为地质体密度，g为万有引力常量，v为引力位。令引力位频谱为sv，对引力位频谱进行x和z方向求导：对进行二阶求导：将上式代入(20)与式(21)可得：根据式(22)与式(23)可以得到，sz，sh相互转换的分量转换因子为：sz→sh：i(24)sh→sz：-i(25)二度体新磁化方向磁位为：α2，γ2为新磁化方向余弦，t1表示原磁化方向，即新磁位是沿着原磁化方向t1反方向的曲线积分ha频谱表达式为：za频谱表达式为：利用式(24)、式(25)得到的分量转换因子可得：由式(26)式(28)式(29)可得：其中，α1，γ1为原磁化方向余弦；根据磁位与分量之间的关系可以推出新分量的频谱表达式：因此对于井位于异常体正方向：任意磁化方向转换转换因子为：其中，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；对于三度体来说，由于三度体hay不为零，故其位场转换因子是由三维拉普拉斯方程狄里希莱问题推导而来：井中一维数据无法作为边界条件来解上述方程，是三维拉普拉斯方程狄里希莱问题的边界条件，因此，直接利用二度体的转换因子来处理三度体数据，相当于把三度体剖面默认为了二度体剖面，结果会存在误差，降低误差的方法包括选取多个井估算结果的平均值作为最终估计结果和选取hay接近于零的一组三分量进行估算。进一步的，s7中相关性计算具体如下：cov表示协方差，d表示方差，za表示垂直分量，ta表示磁异常模量，n表示za和ta的相关性值，hax、hay与ta计算相关性同理。进一步的，s9中所述有效磁化倾角指定度数设置为0°到360°之间。本发明提供的技术方案带来的有益效果是：(1)快速高效的利用测井数据估计出异常体的剖面有效磁化倾角，尤其对于二度体估计结果准确；(2)对于三度体估算结果存在误差，通过多井数据估计结果求平均值的方式减少误差，实用性强。附图说明图1是本发明一种井中磁力异常估计有效磁化倾角的方法流程图；图2是圆柱体模型井位剖面图；图3是球体模型井位俯视图；图4是圆柱体模型有效磁化倾角估计结果图；图5是球体模型有效磁化倾角估计结果图。具体实施方式一种井中磁力异常估计有效磁化倾角的方法，如图1所示，包括以下步骤：s1、基于井中磁异常三分量包括垂直分量za、水平分量hax和水平分量hay数据计算得到磁异常模量ta；公式如下：其中，ta表示磁异常模量，za表示垂直分量，hax表示x方向上的水平分量，hay表示y方向上的水平分量。本实施例在建立的圆柱体模型和球体模型的基础上进行，模型实际磁化方向，剖面有效磁化倾角以及井和模型的相对位置关系分别如图2以及图3所示，图中n方向为剖面横坐标轴的正方向；s2、对za、hax和hay三个分量进行傅里叶正变换分别得到相应的频谱；s3、设置初始有效磁化倾角度数，在本实施例中，x°＝1°；s4、基于s3设置的有效磁化倾角度数以及井与异常体的相对位置关系计算出垂直磁化方向转换因子以及水平磁化方向转换因子；频率域转换因子与井和异常体的相对位置有关，具体关系为：在建立的坐标系下，所研究剖面中，当井位于异常体正方向即在对应横坐标轴正方向时，磁化方向转换因子为：α2＝cosis2γ2＝sinis2α1＝cosis1γ1＝sinis1其中，w表示磁化方向转化因子，is2为转换后磁化倾角，is1为原磁化倾角；当井位于异常体负方向即横坐标轴负方向时，转换因子为：α2＝cosis2γ2＝sinis2α1＝cosis1γ1＝sinis1其中，is2为转换后磁化倾角，is1为原磁化倾角，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；当is2为90°时，转换因子为垂直磁化方向转换因子；当is2为0°时，转换因子为水平磁化方向转换因子；如下表所示为井中二度体频率域位场转换因子，表1井中二度体频率域位场转换因子垂直磁化方向转换因子以及水平磁化方向转换因子的具体推导过程如下：首先，根据地质体三个方向的延伸度，把磁性体分为二度体和三度体，对于地面来说，数据采集面和线永远处在异常体的上半空间，即在所建立的空间直角坐标系下，数据采集在z＝0的平面进行，异常体中心在z＝300的平面，z轴竖直向下为正，对于井来说，假设异常体为二度体，令建立的坐标系x轴垂直于异常体走向，y轴平行于异常体走向，数据采集沿着井轴即z轴方向进行，由于井轴可以在异常体的正方向即正x方向也可以在异常体的负方向即负x方向，因此，井中数据与井和异常体的相对位置有关，所以，在推导二度体的频率域转换因子时需要分两种情况，即在所建立的空间直角坐标系下，井在异常体的正方向和负方向两种情况；对于二度体，水平分量hay为零，其位场转换因子是由二维拉普拉斯方程狄里希莱问题推导而来：其中，x和z表示x轴和z轴坐标，uxx为函数u对x的二次求导，uzz为函数u对z的二次求导，为边界条件即测量数据；情况1:井位于异常体负方向：对于这类情况，二维拉普拉斯方程狄里希莱问题解为：其中，u为狄里希莱问题的解，δ表示边界上任意一点的z轴坐标，t(x，z)满足上式，t为磁异常，x和z表示x轴和z轴坐标，即：对t(0，z)与做关于z的一维傅里叶变换：得到磁异常频谱表达式为：其中，st为磁异常频谱；f为频率，i为虚数单位；利用微分定理对式(1)进行求微分：可以得到对于情况1：x方向导数因子为2πf，z方向导数因子为2πif；给出二度体泊松公式：za磁异常垂直分量，ha为磁异常水平分量，m为磁化强度的模，α，γ为磁化强度的方向余弦，δ为地质体密度，g为万有引力常量，v为引力位；vxz和vzx为对v的x方向和z方向的二次偏导数，vzz为对v的z方向的二次导数,vxx为对v的x方向的二次导数；令引力位频谱为sv，对引力位频谱进行x和z方向求导：svz＝2πifsv对进行二阶求导：将上式代入(4)与式(5)可得：其中，sz为磁异常垂直分量频谱，sh为磁异常水平分量频谱；根据式(6)与式(7)可以得到，sz，sh相互转换的分量转换因子为：sz→sh：-i(8)sh→sz：i(9)二度体新磁化方向磁位即α2，γ2为新磁化方向余弦为：u2是新磁化方向的磁位，t1表示原磁化方向，即新磁位是沿着原磁化方向t1反方向的曲线积分令sh为ha的频谱，ha频谱表达式为：令sz为za的频谱，za频谱表达式为：利用式(8)、式(9)得到的分量转换因子可得：由式(10)式(11)式(13)可得：其中，α1，γ1为原磁化方向余弦；根据磁位与分量之间的关系可以推出新分量(ha2，za2)的频谱表达式：因此对于井位于异常体负方向:任意磁化方向转换因子为：其中，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；情况2:井位于异常体正方向:对于这类情况，二维拉普拉斯方程狄里希莱问题解为：t(x，z)满足上式，t为磁异常，即：对t(0，z)与做关于z的一维傅里叶变换：得到磁场频谱表达式为：st为磁异常频谱，f为频率，i为虚数单位；利用微分定理对式(17)进行求微分：可以得到对于情况2：x方向导数因子为-2πf，z方向导数因子为2πif给出二度体泊松公式：za磁异常垂直分量，ha为磁异常水平分量，m为磁化强度的模，α，γ为磁化强度的方向余弦，δ为地质体密度，g为万有引力常量，v为引力位。令引力位频谱为sv，对引力位频谱进行x和z方向求导：对进行二阶求导：将上式代入(20)与式(21)可得：根据式(22)与式(23)可以得到，sz，sh相互转换的分量转换因子为：sz→sh：i(24)sh→sz：-i(25)二度体新磁化方向磁位为：α2，γ2为新磁化方向余弦，t1表示原磁化方向，即新磁位是沿着原磁化方向t1反方向的曲线积分ha频谱表达式为：za频谱表达式为：利用式(24)、式(25)得到的分量转换因子可得：由式(26)式(28)式(29)可得：其中，α1，γ1为原磁化方向余弦；根据磁位与分量之间的关系可以推出新分量的频谱表达式：因此对于井位于异常体正方向：任意磁化方向转换因子为：其中，α1，γ1为原磁化方向余弦，α2，γ2为新磁化方向余弦；对于三度体来说，由于三度体hay不为零，故其位场转换因子是由三维拉普拉斯方程狄里希莱问题推导而来：井中一维数据无法作为边界条件来解上述方程，是三维拉普拉斯方程狄里希莱问题的边界条件，因此，直接利用二度体的转换因子来处理三度体数据，相当于把三度体剖面默认为了二度体剖面，结果会存在误差，降低误差的方法包括选取多个井估算结果的平均值作为最终估计结果和选取hay接近于零的一组三分量进行估算。s5、za频谱乘以垂直磁化方向转换因子，hax和hay分别乘以水平磁化方向转换因子；s6、将乘以转换因子的频谱分别进行傅里叶逆变换得到磁化方向转到垂直方向的za数据以及磁化方向转到水平方向的hax和hay数据；s7、将s6中转换后的za、hax和hay分别与s1中磁异常模量ta进行相关性计算，记录得到的相关性数值；相关性计算具体如下：cov表示协方差，d表示方差，za表示垂直分量，ta表示磁异常模量，n表示za和ta的相关性值，hax、hay与ta计算相关性同理。s8、将有效磁化倾角增加n°，在本实施例中，n°＝1°，循环s4至s8。s9、当有效磁化倾角增加到指定的度数m°，停止循环，在本实施例中，m°＝360°；基于记录的相关性数值绘制出za、hax和hay的三条相关性曲线，寻找za相关性曲线的最小值点以及hax和hay相关性曲线的最大值点，得到的最值点对应的度数即为估计的有效磁化倾角。在本实施例中，圆柱体模型有效磁化倾角估计结果如图4所示，表2为圆柱体模型有效磁化倾角估计误差表。球体模型有效磁化倾角估计结果如图5所示，表3为球体模型有效磁化倾角估计误差表。表2圆柱体模型有效磁化倾角估计误差表(实际值：54.73°)井号za估计值误差ha估计值误差n155°0.27°54°0.73°n256°1.27°53°1.73°n358°3.27°52°2.73°s155°0.27°54°0.73°s256°1.27°53°1.73°s358°3.27°51°2.73°平均1.6°1.7°表3球体模型有效磁化倾角估计误差表(实际值：63°)井号za估计值误差hax(hay)估计值误差128°35°33°30°224°39°33°30°347°16°162°99°(舍去)465°2°84°21°527°36°33°30°624°39°32°31°745°18°16299°(舍去)866°3°85°22°平均40.8°22.2°50°13°当前第1页12