一种二维网格数据的补偿方向滤波方法
【专利摘要】本发明涉及地球物理勘探的重力、磁力勘探数据处理技术，是一种二维网格数据的补偿方向滤波方法，将原始数据进行网格化处理，在窗口内划分方位角和每个方位对应的窗口内的节点分布，设置每个方位上背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数，计算窗口内滤波方位的背景场值，循环迭代完成补偿方向滤波。本发明滤波的方向性强，满足了重磁处理提取方向异常信息的需要，同时也能满足压制某方向异常信息的需要。本发明推动了重磁勘探数据处理解释技术的进步，同时也为其它学科的二维网格数据的方向滤波处理提供了新方法。
【专利说明】一种二维网格数据的补偿方向滤波方法
【技术领域】
[0001]本发明涉及地球物理勘探的重力、磁力勘探数据处理领域，属于数字滤波技术，是一种二维网格数据的补偿方向滤波方法。
【背景技术】
[0002]在磁力(包括航磁)资料采集中，由于测线间的一致性差异等原因引起的沿线的线性异常是具有方向的，这些方向性异常是地质研究的干扰(即数据噪音)，因而，在重磁勘探资料处理解释中需要提取具有方向性的异常信息，以满足解释断层、构造的需要，需要在资料处理解释前进行处理，即利用方向滤波方法压制重磁数据中方向性的重磁噪音数据，使得滤波处理后的重磁数据更好地反映地质信息是需要解决的问题。
[0003]目前地球物理勘探的重力、磁力勘探数据处理中常用的频率域椭圆方向滤波是在频率域内实现，它的作用是提取出沿选择方向的方向性异常。这种方法常见的步骤是:
1.将重力(或磁力)数据网格化，2.对重力网格数据进行傅立叶变换，3.选择滤波方向，选取椭圆方向滤波器，4.对傅立叶变换后的数据进行椭圆方向滤波，获得滤波后的频谱数据，
5.对滤波后的频谱数据进行反傅立叶变换，即获得椭圆方向滤波后的重力(或磁力)异常网格数据。这种方法不具有迭代性和局部场补偿作用，而且方向性差。
[0004]地质出版社，2005.6，P194_204公开了一种《重力场与重力勘探》方法，利用空间域滑动窗口加权平均法建立一种圆滑滤波器，对数据具有圆滑滤波作用，可获得区域场，但是这种方法的方向滤波方法性能差，方向性异常往往被圆滑，无法提取重磁勘探资料处理解释中具有方向性的异常信息，不能利用方向滤波方法压制重磁数据中方向性的重磁噪音数据。
【发明内容】
[0005]本发明目的在于为重磁数据处理提供一种具有较好方向滤波功能，可提取方向性重磁异常数据或压制具有方向性的噪音数据的二维网格数据的补偿方向滤波方法。
[0006]本发明通过以下步骤实现:
[0007]I)进行野外重力或磁力勘探，采集得到原始数据并进行网格化处理；
[0008]2)运算窗口采用7行7列窗口，迭代次数取I?10 ;
[0009]步骤2)所述的运算窗口大小为7行7列。
[0010]3)在窗口内划分方位角和每个方位对应的窗口内的节点分布；
[0011]步骤3)所述的窗口内划分方位角是划分为南和北方位、东和西方位、北东和西南方位、北西和东南方位、北北东和南南西方位、北东东和南西西方位、北北西和南南东方位、北西西和南东东方位16个。
[0012]步骤3)所述的窗口内每对方位的节点分布相同，按以下分布:
[0013]南和北这对方位的节点分布是一组，窗口节点的行序为从下到上，窗口节点的列的序为从左到右，坐标组行序，列序为:南和北方位的窗口内节点分布为(I，4)、(2,4)、(3,4)、(4，4)、(5，4)、(6，4)、(7，4)；
[0014]东和西方位的节点分布是一组，东和西方位的窗口内节点分布为(4，I)、(4,2),(4，3)、(4，4)、(4，5)、(4，6)、(4，7)；
[0015]北东和南西方位的节点分布是一组，北东和南西方位的窗口内节点分布为(1，I)、(2，2)、(3，3)、(4，4)、(5，5)、(6，6)、(7，7)；
[0016]北西和南东方位的节点分布是一组，北西和南东方位的窗口内节点分布为(7，I )、(6，2)、(5，3)、(4，4)、(3，5)、(2，6)、(1，7)；
[0017]北北东和南南西方位的节点分布是一组(图1的5)，北北东和南南西方位的窗口内节点分布为(1，3)、(2，3)、(3，3)、(3，4)、(4，4)、(5，4)、(5，5)、(6，5)、(7，5)；
[0018]北东东和南西西方位的节点分布是一组(图1的6)，北东东和南西西方位的窗口内节点分布为(3，I)、(3，2)、(3，3)、(4，3)、(4，4)、(4，5)、(5，5)、(5，6)、(5，7)；
[0019]北北西和南南东方位的节点分布是一组，北北西和南南东方位的窗口内节点分布为(1，5)、(2，5)、(3，5)、(3，4)、(4，4)、(5，4)、(5，3)、(6，3)、(7，3)；
[0020]北西西和南东东方位的节点分布是一组，北西西和南东东方位的窗口内节点分布为(5，I)、(5，2)、(5，3)、(4，3)、(4，4)、(4，5)、(3，5)、(3，6)、(3，7)。
[0021]4)设置每个方位上背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数；
[0022]步骤4)所述的节点系数为:
[0023]东西向滤波时，背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数为:(4，1)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)节点的背景场滤波时节点系数自左向右分别为 3、5、7、12、7、5、3 ;剩余场偿滤波求取时的(1，4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)节点系数分别为1、1、1、1、1、1、1 ；
[0024]南北向滤波时，背景场滤波(1，4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)节点系数分别为 3、5、7、12、7、5、3，剩余场偿滤波(4，1)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)节点的系数分别为1、1、1、1、1、1、1 ；
[0025]其它各个方位的背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数设置采用东西向滤波和南北向滤波设置方法类推。
[0026]5)按照以下公式计算窗口内滤波方位的背景场值gw:
[0027]gbJ = (I/ Σ (k (i, j))) X (Σ (k (i, j) Xg(i, j)))
[0028]式中:
[0029]k (i, j)为节点(i，j)的系数值，g(l, j)为节点(i，j)的重力或磁力值；
[0030]步骤5)所述的滤波方位的背景场值gw，当东西向滤波时，先计算计算点(4，4)的背景场值gw (4，4):
[0031]gbJ (4，4)= (1/42)X (3Xg(4, I)+5Xg(4，2)+7Xg(4，3)+12Xg(4，4)+7Xg(4，5) +5Xg(4, 6)+3Xg(4, 7))；
[0032]式中g(4，I)、g(4，2)、g(4，3)、g(4，4)、g(4，5)、g(4，6)、g(4，7)分别表示计算窗 P内节点(4，I)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)的重力值或者是磁力值；
[0033]依次移动运算窗口，计算工区平面网格的每一个计算点的背景场值；
[0034]其它各方位滤波，替换所选坐标节点和对应节点的系数、按各方位的系数加权平均计算。
[0035]6)按照以下公式计算窗口内滤波方位的剩余场补偿滤波值gbc;h:
[0036]gbch= (1/ Σ (kbch (i, j))) X (Σ (kbch (i, j) Xgbch(i, j)))
[0037]式中:
[0038]kbch (i, j)为节点(i，j)的系数值kbc;h，g(i, j)为节点(i，j)的重力或磁力值；
[0039]步骤6)所述的剩余场补偿滤波值gbc;h,东西向滤波时:
[0040]gbch (4,4) =gbJ (4,4)-( (1/7) X (gbJ(1, 4) +gbJ (2，4) +gbJ (3，4) +gbJ (4，4) +gbJ (5，4) +gbJ(6,4)+gbJ(7,4)))；
[0041]然后,计算计算点的补偿方向滤波输出值gf(4,4),
[0042]gf (4,4) =g (4, 4) -gbJ (4,4) +gbch (4,4)；
[0043]再计算该点的剩余异常值gs(4，4)，
[0044]gs (4，4) =g (4，4) -gf (4,4)；
[0045]依次移动运算窗口，计算每一个计算点的补偿方向滤波值gfl(i，j)、剩余异常值
Ss1 (i，j)。
[0046]7)循环迭代直到完成迭代运算次数，完成补偿方向滤波。
[0047]所述的循环迭代是将第一次计算得到的各个节点的补偿方向滤波值gf计为gfl (i，j)，将第一次计算得到的各个节点的剩余异常值gs计为gs1 (i，j)，对gs1 (i，j)进行下一次的补偿方向滤波运算，即把gs1(i，j)视为步骤5)中的原始数据g(i，j)，进行迭代，计算步骤5)、步骤6)，获得新的g’ f (i，j)、g’ s(i，j)，则第二次迭代的结果为:
[0048]gf2(i, j)=gfl(i，j)+g，f(i，j).[0049]gs2(i, j)=gsl(i，j)_g，f(i，j).[0050]式中(i，j)表示计算点的节点在工区坐标，g’ f、g’ s分别表示本次迭代中计算得到的计算点的补偿方向滤波、剩余异常值。
[0051]本发明对实测数据处理，获得了多次迭代后的补偿方向滤波异常和补偿方向滤波剩余异常，滤波的方向性强，满足了重磁处理提取方向异常信息的需要，同时也能满足压制某方向异常信息的需要。本发明推动了重磁勘探数据处理解释技术的进步，同时也为其它学科的二维网格数据的方向滤波处理提供了一种新方法。
【专利附图】
【附图说明】
[0052]图1是本发明二维网格的窗口选取和16方位节点分布图；
[0053]图2是本发明背景场滤波节点系数和补偿场滤波节点系数分布图；
[0054]图3是原始磁异常平面图；
[0055]图4是本发明南北向滤波磁异常平面图；
[0056]图5是本发明南北向滤波剩余磁异常平面图。
【具体实施方式】
[0057]以下结合实例详细说明本发明。
[0058]本实例通过以下步骤实现:
[0059]1)进行野外航空磁力勘探，采集得到航空磁力原始数据并进行网格化处理，见图3 ；
[0060]2)运算窗口采用7行7列窗口(见图1)，迭代次数取2;
[0061]步骤2)所述的运算窗口大小为7行7列，见图1。
[0062]3)在窗口内划分方位角和每个方位对应的窗口内的节点分布，见图1 ;
[0063]步骤3)所述的窗口内划分方位角是划分为南和北方位、东和西方位、北东和西南方位、北西和东南方位、北北东和南南西方位、北东东和南西西方位、北北西和南南东方位、北西西和南东东方位16个，见图1。
[0064]步骤3)所述的窗口内每对方位的节点分布相同，按以下分布:
[0065]南和北这对方位的节点分布是一组，见图1的1，窗口节点的行序为从下到上，窗口节点的列的序为从左到右，坐标组行序，列序为:南和北方位的窗口内节点分布为(I，
4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)；
[0066]东和西方位的 节点分布是一组，见图1的2，东和西方位的窗口内节点分布为(4，I)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)；
[0067]北东和南西方位的节点分布是一组，见图1的3，北东和南西方位的窗口内节点分布为(1，I)、(2，2)、(3，3)、(4，4)、(5，5)、(6，6)、(7，7)；
[0068]北西和南东方位的 节点分布是一组，见图1的4，北西和南东方位的窗口内节点分布为(7，I)、(6，2)、(5，3)、(4，4)、(3，5)、(2，6)、(1，7)；
[0069]北北东和南南西方位的节点分布是一组，见图1的5，北北东和南南西方位的窗口内节点分布为(1，3)、(2，3)、(3，3)、(3，4)、(4，4)、(5，4)、(5，5)、(6，5)、(7，5)；
[0070]北东东和南西西方位的节点分布是一组，见图1的6，北东东和南西西方位的窗口内节点分布为(3，I)、(3，2)、(3，3)、(4，3)、(4，4)、(4，5)、(5，5)、(5，6)、(5，7)；
[0071]北北西和南南东方位的节点分布是一组，见图1的7，北北西和南南东方位的窗口内节点分布为(1，5)、(2，5)、(3，5)、(3，4)、(4，4)、(5，4)、(5，3)、(6，3)、(7，3)；
[0072]北西西和南东东方位的节点分布是一组，见图1的8，北西西和南东东方位的窗口内节点分布为(5，I)、(5，2)、(5，3)、(4，3)、(4，4)、(4，5)、(3，5)、(3，6)、(3，7)。
[0073]4)设置南和北方位上背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数；
[0074]图2是背景场滤波节点系数和补偿场滤波节点系数分布图。图中I为东-西向，2为南-北向，3为北西-东南向,4为东北-西南向，5为北西西-南南东向,6为北北西-南南东向，7为北东东-西西南向，8为北北东-南南西向。
[0075]步骤4)所述的节点系数为:
[0076]南北向滤波时，背景场滤波(1，4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)节点系数分别为 3、5、7、12、7、5、3，剩余场偿滤波(4，1)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，
7)节点的系数分别为1、1、1、1、1、1、1，节点系数分布见图2的I ;
[0077]5)按照以下公式计算窗口内滤波方位的背景场值gw:
[0078]gbJ = (I/ Σ (k (i, j))) X (Σ (k (i, j) Xg(i, j)))
[0079]式中:k (i，j)为节点(i，j)的系数值，g(i，j)为节点(i，j)的重力或磁力值；
[0080]步骤5)所述的滤波方位的背景场值gw，南北方向滤波时，先计算计算点(4，4)的背景场值gw (4，4):[0081]
【权利要求】
1.一种二维网格数据的补偿方向滤波方法，特点是通过以下步骤实现: 1)进行野外重力或磁力勘探，采集得到原始数据并进行网格化处理； 2)运算窗口采用7行7列窗口，迭代次数取I~10; 3)在窗口内划分方位角和每个方位对应的窗口内的节点分布； 4)设置每个方位上背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数； 5)按照以下公式计算窗口内滤波方位的背景场值gw: gbJ =(1/ Σ (k (i, j))) X (Σ (k (i, j) Xg(i, j))) 式中: k (i, j)为节点(i，j)的系数值，g(i, j)为节点(i，j)的重力或磁力值； 6)按照以下公式计算窗口内滤波方位的剩余场补偿滤波值gbc;h:
gbch= (I/ Σ (kbch (i，j)))X (Σ (kbch (i，j) Xgbch(i，j))) 式中: kbch (i, j)为节点(i，j)的系数值kbc;h，g(i, j)为节点(i，j)的重力或磁力值； 7 )循环迭代直到完成迭代运算次数，完成补偿方向滤波。
2.根据权利要求1所述的方法，特点是步骤2)所述的运算窗口大小为7行7列。`
3.根据权利要求1所述的方法，特点是步骤3)所述的窗口内划分方位角是划分为南和北方位、东和西方位、北东和西南方位、北西和东南方位、北北东和南南西方位、北东东和南西西方位、北北西和南南东方位、北西西和南东东方位16个。
4.根据权利要求1所述的方法，特点是步骤3)所述的窗口内每对方位的节点分布相同，按以下分布: 南和北这对方位的节点分布是一组，窗口节点的行序为从下到上，窗口节点的列的序为从左到右，坐标组行序，列序为:南和北方位的窗口内节点分布为(1，4)、(2,4)、(3,4)、(4，4)、(5，4)、(6，4)、(7，4);东和西方位的节点分布是一组，东和西方位的窗口内节点分布为(4，I)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)； 北东和南西方位的节点分布是一组，北东和南西方位的窗口内节点分布为(1，I)、(2，2)、(3，3)、(4，4)、(5，5)、(6，6)、(7，7)； 北西和南东方位的节点分布是一组，北西和南东方位的窗口内节点分布为(7，I)、(6，2)、(5，3)、(4，4)、(3，5)、(2，6)、(1，7)； 北北东和南南西方位的节点分布是一组(图1的5)，北北东和南南西方位的窗口内节点分布为(1，3)、(2，3)、(3，3)、(3，4)、(4，4)、(5，4)、(5，5)、(6，5)、(7，5)； 北东东和南西西方位的节点分布是一组(图1的6)，北东东和南西西方位的窗口内节点分布为(3，I)、(3，2)、(3，3)、(4，3)、(4，4)、(4，5)、(5，5)、(5，6)、(5，7)； 北北西和南南东方位的节点分布是一组，北北西和南南东方位的窗口内节点分布为(1，5)、(2，5)、(3，5)、(3，4)、(4，4)、(5，4)、(5，3)、(6，3)、(7，3)； 北西西和南东东方位的节点分布是一组，北西西和南东东方位的窗口内节点分布为(5，I)、(5，2)、(5，3)、(4，3)、(4，4)、(4，5)、(3，5)、(3，6)、(3，7)。
5.根据权利要求1所述的方法，特点是步骤4)所述的节点系数为: 东西向滤波时，背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数为:(4，I)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)节点的背景场滤波时节点系数自左向右分别为 3、5、7、12、7、5、3 ;剩余场偿滤波求取时的(1，4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)节点系数分别为1、1、1、1、1、1、1 ； 南北向滤波时，背景场滤波(1，4)、(2，4)、(3，4)、(4，4)、(5，4)、(6，4)、(7，4)节点系数分别为 3、5、7、12、7、5、3，剩余场偿滤波(4，1)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)节点的系数分别为1、1、1、1、1、1、I ； 其它各个方位的背景场滤波时节点系数和剩余场补偿滤波求取时的节点系数设置采用东西向滤波和南北向滤波设置方法类推。
6.根据权利要求1所述的方法，特点是步骤5)所述的滤波方位的背景场值gw，当东西向滤波时，先计算计算点(4，4)的背景场值gw (4，4):
gbJ(4，4)= (1/42)X (3Xg(4, I)+5Xg(4, 2)+7Xg(4, 3)+12Xg(4, 4) +7Xg(4, 5)+5Xg(4，6)+3 Xg (4, 7))； 式中 g(4, I)、g(4, 2)、g(4, 3)、g(4, 4)、g(4, 5)、g(4, 6)、g(4, 7)分别表示计算窗口内节点(4，I)、(4，2)、(4，3)、(4，4)、(4，5)、(4，6)、(4，7)的重力值或者是磁力值； 依次移动运算窗口，计算工区平面网格的每一个计算点的背景场值； 其它各方位滤波，替换所选坐标节点和对应节点的系数、按各方位的系数加权平均计笪
7.根据权利要求1所述的方法，特点是步骤6)所述的剩余场补偿滤波值gbc;h，东西向滤波时:`
gbch (4，4) =gbJ (4,4)-( (1/7) X (gbj (I, 4) +gbJ (2，4) +gbJ (3，4) +gbJ (4，4) +gbJ (5，4) +gbJ (6，4)+gbJ(7,4)))； 然后，计算计算点的补偿方向滤波输出值gf (4，4)， gf (4，4) =g (4，4) -gbJ (4，4) +gbch (4,4)； 再计算该点的剩余异常值gs(4，4)， gs (4,4) =g (4, 4) -gf (4,4)； 依次移动运算窗口，计算每一个计算点的补偿方向滤波值gfl(i，j)、剩余异常值gsl(i，j)。
8.根据权利要求1所述的方法，特点是步骤7)所述的循环迭代是将第一次计算得到的各个节点的补偿方向滤波值gf计为gfl(i，j)，将第一次计算得到的各个节点的剩余异常值gs计为gsi(i，j)，对gsl(i，j)进行下一次的补偿方向滤波运算，即把gsl(i，j)视为步骤5)中的原始数据g(i，j)，进行迭代，计算步骤5)、步骤6)，获得新的g’ f(i，j)、g’ s(i，j)，则第二次迭代的结果为:
gf2(i，j)=gn(i，j)+g，f(i，j);
gs2(i，j)=gsi(i，j)-g，f(i，j); 式中(i，j)表示计算点的节点在工区坐标，g’ f、g’s分别表示本次迭代中计算得到的计算点的补偿方向滤波、剩余异常值。
【文档编号】G01V7/02GK103513288SQ201210211223
【公开日】2014年1月15日 申请日期:2012年6月21日 优先权日:2012年6月21日
【发明者】刘云祥 申请人:中国石油天然气集团公司, 中国石油集团东方地球物理勘探有限责任公司