一种大距离位场向上延拓计算方法与流程

本发明涉及应用于重磁数据处理领域的位场向上延拓计算方法，特别是大距离位场向上延拓快速、高精度数值计算方法。
背景技术：
：重力场和磁场统称位场。位场延拓是重磁数据处理中的一种重要方法，包括向上延拓、向下延拓、曲化平、平化曲和平化平。向上延拓在位场延拓问题中最简单，也是最重要的一种，它是实现其它几类位场延拓问题的基础。传统向上延拓方法是基于快速傅里叶变换算法(fft)的波数域向上延拓方法，通常称为向上延拓fft法。fft法借助快速傅里叶变换算法的高效性，能够快速实现大规模位场数据的向上延拓。但是，由于存在离散傅里叶变换边界效应问题，fft法大距离向上延拓计算精度低，通常采用某种扩边方法来减小边界效应，但扩边也会引入数据误差，降低延拓结果精度。文献(汪炳柱.用样条函数法求重力异常二阶垂向导数和向上延拓计算.石油地球物理勘探,1996,31(3):415-422.)公开了一种重力异常向上延拓的样条函数法，该方法延拓精度高但是计算效率低。文献(陈龙伟,胡小平,吴美平,马涛.空间域位场延拓新方法研究.地球物理学进展,2012,27(4):1509-1518.)和文献(yilezhang,yaushuwong,yuanfanglin.bttb-rrcgmethodfordownwardcontinuationofpotentialfielddata.journalofappliedgeophysics,2016,126:74-86.)分别公开了一种相同的空间域向上延拓方法，根据系数矩阵具有bttb矩阵结构的特点，采用bttb矩阵与向量相乘的快速算法，实现了向上延拓快速计算，并且有效克服了fft法存在的边界效应问题。该方法存在的问题是向上延拓距离小时，系数矩阵接近奇异，延拓结果不稳定。目前已有的向上延拓方法存在的最大问题是不能同时保证计算效率和计算精度，无法满足快速、高精度、大距离向上延拓的要求。因此，寻找一种计算效率高、同时能保证计算精度的位场延拓方法具有重要的现实意义。技术实现要素：本发明针对传统位场基于快速傅里叶变换算法(fft)的波数域向上延拓方法因边界效应而导致向上延拓计算精度低的问题，提出了一种大距离位场向上延拓计算方法，其是一种新的空间域大距离位场向上延拓快速、高精度数值计算方法。为解决上述技术问题，本发明采用以下技术方案：一种大距离位场向上延拓计算方法，其步骤为：步骤一：设定应用场景建立直角坐标系，z轴垂直向下为正。给定高度为z0的规则网格化位场数据图，规则网格化位场数据图其x方向上等间距分布有m个网格节点，其y方向上等间距分布有n个网格节点，这样规则网格化位场数据图上有m×n个网格节点，m×n个网格节点对应m×n个规则网格化位场数据。给定向上延拓高度δz，δz为正数；称高度为z0的水平面为观测面，称高度为z0-δz的水平面为延拓面；高度为z0的规则网格化位场数据图上的m×n个规则网格化位场数据记为u(ξi,ηj,z0)，i＝1,2,…,m，j＝1,2,…,n；(ξi,ηj,z0)表示观测面网格节点坐标，u(ξi,ηj,z0)表示网格节点坐标(ξi,ηj,z0)处观测面位场值；m表示规则网格化位场数据图其x方向上的网格节点数，n表示规则网格化位场数据图其y方向上的网格节点数。高度为z0-δz的延拓面上的规则网格化位场数据图上的m×n个规则网格化位场数据记为u(xm,yn,z0-δz)，(xm,yn,z0-δz)表示延拓面网格节点坐标，u(xm,yn,z0-δz)表示网格节点坐标(xm,yn,z0-δz)处延拓面位场值，m＝1,2,…,m，n＝1,2,…,n。步骤二：根据公式(1)，计算加权系数式(1)中，ω(xm,yn；ξi,ηj)表示加权系数；(xm,yn)表示延拓面网格节点坐标；(ξi,ηj)表示观测面网格节点坐标；arctan()表示反正切函数；δz表示向上延拓高度；μpq，xp，yq，rpq(p＝1,2；q＝1,2)表示辅助计算变量，其计算式为x1＝xm-ξi+0.5δx，x2＝xm-ξi-0.5δxy1＝ym-ηi+0.5δy，y2＝ym-ηi-0.5δy式中，δx表示x方向上相邻网格节点间的间距，δy表示y方向上相邻网格节点间的间距；步骤三：位场向上延拓：根据公式(1)给出的加权系数和一种二维离散卷积快速算法，实现公式(2)给出的位场向上延拓快速、高精度计算式中，(xm,yn,z0-δz)表示延拓面网格坐标，m＝1,2,…,m，n＝1,2,…,n；(ξi,ηj,z0)表示观测面网格坐标，i＝1,2,…,m，j＝1,2,…,n；z0-δz表示延拓面高度；z0表示观测面高度；u(xm,yn,z0-δz)表示网格坐标(xm,yn,z0-δz)处延拓面位场值；u(ξi,ηj,z0)表示网格坐标(ξi,ηj,z0)处观测面位场值。与现有位场延拓方法相比，采用公式(1)给出的方法计算加权系数，能够提高向上延拓的计算精度，并且保证小距离向上延拓时数值计算的稳定性。步骤三中，根据公式(1)给出的加权系数和一种二维离散卷积快速算法，实现公式(2)给出的位场向上延拓快速、高精度计算，其步骤为：(1)将加权系数ω(x1,y1；ξi,ηj)排列成矩阵t，记为式(3)中，矩阵t中的各矩阵元素用ti,j表示；矩阵t中的矩阵元素ti,j与加权系数ω(x1,y1；ξi,ηj)存在关系ti,j＝ω(x1,y1；ξi,ηj)(4)(2)将矩阵t补零扩展成矩阵将矩阵分成四个块矩阵，记为将块矩阵互换位置，得到矩阵cext(3)将观测面位场数据u(ξi,ηj,z0)(i＝1,2,…,m，j＝1,2,…,n)排列成矩阵g，矩阵元素gi,j与位场数据存在关系gi,j＝u(ξi,ηj,z0)(8)将矩阵g补零扩展成矩阵gext式中，表示nx×ny零矩阵；(4)计算式中，fft2()表示二维快速傅里叶变换；(5)计算式中，“.*”表示对应元素相乘运算；(6)计算式中，ifft2()表示二维快速傅里叶反变换；(7)提取矩阵fext的前m行前n列，构成矩阵f，即为二维离散卷积计算结果，也就是位场向上延拓结果。与现有位场向上延拓方法相比，采用上述二维离散卷积快速算法是本专利发明点之一，它保证了在大规模位场数据情况下向上延拓的高效性。本发明是一个有机整体，根据步骤二提供的一种特殊的加权系数计算公式，采用步骤三运用的二维离散卷积快速算法，实现了大距离位场向上延拓在效率和精度上的统一。本发明解决了传统fft法因边界效应而导致大距离向上延拓计算精度低的问题，为满足大规模位场数据处理、位场分离、大深度位场向下延拓等提供了方法支撑。与现有位场延拓方法相比，本发明具有以下优点：(1)能够实现大规模位场数据情况下大距离向上延拓快速、高精度数值计算；(2)能够保证小距离位场向上延拓情况下数值计算的稳定性；(3)大规模位场数据向上延拓时，算法不但计算效率和计算精度高，并且所需计算机内存小。附图说明：图1为大距离位场向上延拓计算流程图；图2为观测面和延拓面网格坐标示意图；图3为棱柱体磁异常组合模型；图4为z＝0m高度面磁异常理论值；图5为z＝20000m高度面磁异常理论值；图6为z＝20000m高度面fft法向上延拓结果；图7为z＝20000m高度面新方法向上延拓结果。图中符号说明如下：si表示国际单位制具体实施方式：下面结合附图对本发明中的方法作进一步详细描述。参照图1，本发明提供一种大距离位场向上延拓计算方法，包括以下步骤：步骤一：设定应用场景建立直角坐标系，z轴垂直向下为正；给定高度为z0的规则网格化位场数据图。图2所示，规则网格化位场数据图是一张规则的网格图，其x方向上等间距分布有m个网格节点，其y方向上等间距分布有n个网格节点，这样规则网格化位场数据图上有m×n个网格节点，m×n个网格节点对应m×n个规则网格化位场数据。给定向上延拓高度δz，δz为正数；称高度为z0的水平面为观测面，称高度为z0-δz的水平面为延拓面。高度为z0的规则网格化位场数据图上的m×n个规则网格化位场数据记为u(ξi,ηj,z0)，i＝1,2,…,m，j＝1,2,…,n；u(ξi,ηj,z0)表示网格节点坐标(ξi,ηj,z0)处的观测面位场值；(ξi,ηj,z0)表示观测面网格节点坐标，u(ξi,ηj,z0)表示网格节点坐标(ξi,ηj,z0)处的观测面位场值；m表示规则网格化位场数据图其x方向网格节点数，n表示规则网格化位场数据图其y方向网格节点数；高度为z0-δz的延拓面上的规则网格化位场数据图上的m×n个规则网格化位场数据记为u(xm,yn,z0-δz)，(xm,yn,z0-δz)表示延拓面网格节点坐标，u(xm,yn,z0-δz)表示网格节点坐标(xm,yn,z0-δz)处延拓面位场值，m＝1,2,…,m，n＝1,2,…,n；步骤二：根据公式(1)，计算加权系数式(1)中，ω(xm,yn；ξi,ηj)表示加权系数；(xm,yn)表示延拓面网格节点坐标；(ξi,ηj)表示观测面网格节点坐标；arctan()表示反正切函数；δz表示向上延拓高度；μpq，xp，yq，rpq表示辅助计算变量，其计算式为x1＝xm-ξi+0.5δx，x2＝xm-ξi-0.5δxy1＝ym-ηi+0.5δy，y2＝ym-ηi-0.5δy其中：p＝1,2；q＝1,2；δx表示x方向上相邻网格节点间的间距，δy表示y方向上相邻网格节点间的间距；根据公式(1)，计算ω(x1,y1；ξi,ηj)，i＝1,2,…,m，j＝1,2,…,n，将加权系数存储于计算机内存。步骤三：计算位场向上延拓结果。根据公式(1)给出的加权系数和一种二维离散卷积快速算法，计算得到公式(2)给出的位场向上延拓结果。位场向上延拓计算公式为式(2)中，(xm,yn,z0-δz)表示延拓面网格节点坐标，m＝1,2,…,m，n＝1,2,…,n；(ξi,ηj,z0)表示观测面网格节点坐标，i＝1,2,…,m，j＝1,2,…,n；z0-δz表示延拓面高度；z0表示观测面高度；u(xm,yn,z0-δz)表示网格节点坐标(xm,yn,z0-δz)处延拓面位场值；u(ξi,ηj,z0)表示网格节点坐标(ξi,ηj,z0)处观测面位场值；m表示规则网格化位场数据图其x方向上的网格节点数，n表示规则网格化位场数据图其y方向上的网格节点数。采用一种二维离散卷积快速计算方法，实现公式(2)的数值计算，得到位场向上延拓结果，其步骤为：(1)将加权系数ω(x1,y1；ξi,ηj)排列成矩阵t，记为式(3)中，矩阵t中的矩阵元素ti,j与加权系数ω(x1,y1；ξi,ηj)存在关系ti,j＝ω(x1,y1；ξi,ηj)(4)(2)将矩阵t补零扩展成矩阵将矩阵分成四个块矩阵，记为将块矩阵互换位置，得到矩阵cext(3)将观测面位场数据u(ξi,ηj,z0)(i＝1,2,…,m，j＝1,2,…,n)排列成矩阵g，矩阵元素gi,j与位场数据存在关系gi,j＝u(ξi,ηj,z0)(8)将矩阵g补零扩展成矩阵gext式中，表示nx×ny零矩阵；(4)计算式中，fft2()表示二维快速傅里叶变换；(5)计算式中，“.*”表示对应元素相乘运算；(6)计算式中，ifft2()表示二维快速傅里叶反变换；(7)提取矩阵fext的前m行前n列，构成矩阵f，即为二维离散卷积计算结果，也就是位场向上延拓结果。下面对本发明方法的效果进行检验。为了说明本发明所提出的方法用于大距离位场向上延拓时的效率和精度，设计了如下棱柱体磁异常组合模型(图3所示)，模型参数由表1给出。背景磁场强度取为50000nt，磁倾角为50度，磁偏角为0度。观测平面范围为：x方向从0m到50000m，y方向从0m到52000m。模拟1:50000比例尺地磁测量，网格间距取为δx＝δy＝500m，观测点数为101×105，向上延拓高度δz＝20000m，即向上延拓40倍点距。大距离位场向上延拓计算方法的程序代码由matlab语言编写，程序运行计算机为个人笔记本电脑，cpu主频为2.8ghz，内存为32gb。由表2可知，新的向上延拓算法运行所需时间约为0.016秒，传统fft法需0.003秒，由此可见本发明提供的新方法保持了传统fft法的计算效率。参见图4至图7，图4,为z＝0m高度面磁异常理论值；图5为z＝20000m高度面磁异常理论值；图6为z＝20000m高度面采用fft法得到的向上延拓结果；图7为z＝20000m高度面采用本发明方法得到的向上延拓结果。对比可知，从形态上看，本发明方法的延拓结果与理论值吻合很好，但fft法延拓结果出现大的边界效应，与理论值在形态上出现较大差别。以均方根误差衡量延拓结果精度，由表2可知，新方法延拓结果均方根误差为0.28nt，而传统fft法延拓结果均方误差为2.98nt，可见新方法延拓精度远优于fft法，具有很高计算精度。表1异常体模型参数表2计算时间和计算精度计算精度(nt)计算时间(s)fft法2.980.003新方法0.280.016以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本
技术领域：
的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。当前第1页12