一种磁源重力视强度反演成像方法
【专利摘要】本发明涉及一种磁源重力视强度反演成像方法，其包括以下步骤：将磁异常转换成空间域的磁源重力异常；选取不同的切割半径，采用插值切割法对空间域的磁源重力异常进行不同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常；采用迭代法将各空间层的剩余磁源重力异常向下延拓到各空间层的磁性体的顶界面上，获取磁性体顶界面的磁源重力异常；对磁性体顶界面的磁源重力异常的视强度参数进行反演，获取各空间层的视强度参数。本发明具有反演精度高且能够降低反演多解性的优点，本发明可以广泛应用于深部矿产勘探和深部地质构造研究中。
【专利说明】一种磁源重力视强度反演成像方法
【技术领域】
[0001] 本发明涉及一种反演成像方法，特别是关于应用地球物理学领域中的一种磁源重 力视强度反演成像方法。
【背景技术】
[0002] 目前，国内外用于磁场物性（磁化强度）反演的方法分为两类：线性反演和非线性 反演。具体的算法大致分为广义最小二乘法（Gauss法）、最速下降法、阻尼最小二乘法、奇 异值分解法、人工神经网络法BP算法、遗传算法、模拟退火算法和基于Radon变换、小波多 尺度边缘的算法等。上述反演方法在理论创新及小面积实际资料反演应用中都取得了一定 的成绩，但在用于大面积的三维物性反演的实际生产当中没有发挥到应有的作用。其主要 存在问题：1)物性参数与几何参数全局反演时，需要解大型线性代数方程组，计算时间长； 2)求解方程组数学病态十分严重，解的稳定性很差；3)磁源偶极正负异常相互叠加效应导 致异常反演解精度差、分辨率低。
【发明内容】
[0003] 针对上述问题，本发明的目的是提供一种反演精度高且能够降低反演多解性的磁 源重力视强度反演成像方法。
[0004] 为实现上述目的，本发明采取以下技术方案：一种磁源重力视强度反演成像方法， 其包括以下步骤：1)将磁异常AT(x，y)转换成空间域的磁源重力异常Ag(x，y) ;2)选 取不同的切割半径，采用插值切割法对空间域的磁源重力异常Λg(x，y)进行不同深度空 间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常Λgj(x，y，〇) ;3)采用迭代 法将各空间层的剩余磁源重力异常Λg^x，y，0)向下延拓到各空间层的磁性体的顶界面 上，获取磁性体顶界面的磁源重力异常AgO^yjj) ;4)对磁性体顶界面的磁源重力异常 Λg(X，y，hp的视强度参数进行反演，获取各空间层的视强度参数。
[0005] 所述步骤1)中，将磁异常ΛT(x，y)转换成空间域的磁源重力异常Ag(x，y)，其包 括以下步骤：（I)将地面观测到的磁异常ΛT(x，y)转换成频率域的磁异常Δ?，即
[0006] ΔΓ(?5v) =F[^r(x,j)],
[0007] 式中，u和v分别表示x方向和y方向的波数；（II)将频率域的磁异常Δ?^,ν)转 换成化极磁异常之（w，v)，即
【权利要求】
1. 一种磁源重力视强度反演成像方法，其包括以下步骤： 1) 将磁异常ΛT(x，y)转换成空间域的磁源重力异常Ag(x，y); 2) 选取不同的切割半径，采用插值切割法对空间域的磁源重力异常Λg(x，y)进行不 同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常Λgj(X，y，〇); 3) 采用迭代法将各空间层的剩余磁源重力异常Λg^x，y，0)向下延拓到各空间层的磁 性体的顶界面上，获取磁性体顶界面的磁源重力异常Λg(X，y，hp; 4) 对磁性体顶界面的磁源重力异常Λg(x，y，hp的视强度参数进行反演，获取各空间 层的视强度参数。
2. 如权利要求1所述的一种磁源重力视强度反演成像方法，其特征在于：所述步骤1) 中，将磁异常ΛT(X，y)转换成空间域的磁源重力异常Λg(X，y)，其包括以下步骤： (I) 将地面观测到的磁异常AT(x，y)转换成频率域的磁异常Af0，v)，即
式中，u和V分别表示X方向和y方向的波数； (II) 将频率域的磁异常转换成化极磁异常^ ，即
式中，i表示虚分量，（Ptl,Qtl,Rtl)表示磁性体所在地层的方向余弦，（P1,Q1,R1)表示磁性 体在地表的方向余弦，（P2,Q2,R2)表示磁性体在磁倾角为90°的磁化方向下的方向余弦； (III) 将化极磁异常(?,0的偶极磁源作类重力的单极源处理，即对乏沿磁性 体南北极方向积分，得到频率域的磁源重力异常Af(W)，即
式中，M表示磁化强度矢量，G表示万有引力常量，P表示磁性体的剩余密度； (IV) 对频率域的磁源重力异常进行反傅里叶变换，得到空间域的磁源重力异 常Λg(x,y)，即 Ag(x,y) = F'[Ag(?,v)] <?
3. 如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步 骤2)中，获取各空间层的剩余磁源重力异常Λgj (X，y，0)，其具体包括以下步骤： (I) 选取不同的切割半径，其中，η表示切割半径的个数，采用插值切割方 法对磁源重力异常Λg(x，y)进行不同深度空间层的磁源重力异常分离，获取不同切割半 径下的浅源异常： L1 (X，y)，L2 (X，y)，L3 (X，y)…Ln (X，y)， (II) 根据不同切割半径下的浅源异常，得到各空间层的剩余磁源重力异常：
式中，Λgj (x，y，0)表示第j层磁性体在地面反映的磁源重力异常，即剩余磁源重力异 堂 巾。
4.如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步 骤3)中，采用迭代法向下延拓，获取磁性体顶界面的磁源重力异常Λg(x，y，hp，其具体包 括以下步骤： (I) 将地面ΓA上的异常值Λg(x，y，0)中第s点的异常值Gas放到平面ΓΒ的垂直投 影点s上，作为平面ΓB上的异常值Gf，即=G4s ,并通过傅里叶变换得到ΓB上异常 值的频谱Δ#(ιι，ν,0); (II) 当地面ΓΑ与平面ΓΒ之间没有场源时，重力异常满足拉普拉斯方程，采用向上延 拓公式：
将平面ΓΒ上的异常值的频谱Gf代入上式求出地面ΓA上s点的异常值G1，即： (
式中，Δ#(?，ν,0)表示地面ΓΑ(ζ= 0)上的磁源重力异常值Ag(x，y，0)的傅里叶变 换；Γ1表示反傅里叶变换，h表示延拓高度； (III) 利用Gas与GT的差对平面ΓΒ上的异常值进行修正，得到新的异常值Gif. 即： Gf： =Gii;:f 式中，1表示步长； (IV) 重复步骤（II)和步骤（III)，得到如下迭代公式： 当时，由式（16)可知|Gf-<^|<1ε表示一个给定的趋近于零的 数，则GiT" ? ;迭代次数η取20?50次； (V) 采用与步骤（I)?（IV)相同的方法，计算得到平面ΓB上全部点的磁源重力异常 值,并利用向上延拓公式： Ag [x,y,h) = F^i v,〇)J? 作反傅里叶变换得到Λg(x,y, 〇)向下延拓hj深度的磁源重力异常Λg(x,y,hj)。
5.如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步 骤4)中，获取各空间层的视强度参数，其具体包括以下步骤： ⑴将引起磁源重力异常Λg(x，y，hj)的非规则的磁性体分割为M*N个大小规则、质心 位置不同的细长垂直棱柱体； 设棱柱体的质心点坐标为（Xtlk,yd，其剩余密度为Pk，棱柱体的底面边长分别为a和b，棱柱体高度为dh，则在频率域中各棱柱体在顶界面（x，y)点产生的重力异常Agk(x，y) 的频谱为：
式中，r-yjir+v2，e表不自然常数，k表不棱柱体序号，k= 1，2, 3,…，M*N; (II) 将各棱柱体在顶界面（x，y)点产生的重力异常的频谱进行叠加，得到整个深度层 在顶界面（X，y)点产生的重力异常Ag(x，y)的频谱Δ#(Μ,ν)为：
令 …， 则磁源重力异常在顶界面的频谱响应为：
则频域中空间层的剩余磁源重力异常顶界面的视密度为：
(III) 对频域中空间层的剩余磁源重力异常顶界面的视密度进行反傅里叶变 换，得到该深度层顶面各点（x，y)的剩余密度值为：
(IV) 将剩余密度值P(X，y)代入式 κ{χ,ν) = p(x,y)~%? 计算得到磁异常的视磁化率κ(x，y);将计算得到的磁异常的视磁化率Κ (x，y)代入 式 M=K (X，y)*T0, 计算得到磁异常的视磁化强度M; 式中，κ(x，y)为磁异常的视磁化率，M为磁异常的视磁化强度，P(x，y)为磁源重力的 视密度，Ttl为正常磁场强度，κ'和p'分别是泊松公式转换时需要的已知磁化率值和密度 值。
【文档编号】G01V3/40GK104316972SQ201410549240
【公开日】2015年1月28日 申请日期:2014年10月16日 优先权日:2014年10月16日
【发明者】刘春成, 李军, 张益明, 吴兴, 叶云飞, 简兴祥, 汪小将, 王绪本, 牛聪, 贾真, 杨小椿 申请人:中国海洋石油总公司, 中海油研究总院, 成都理工大学