一种基于Toplitze核矩阵的重力场快速正演方法及反演方法与流程

本发明涉及地球物理及勘探技术领域，具体涉及一种基于toplitze核矩阵的重力场快速正演方法及反演方法。

背景技术：

随着高精度高覆盖率全球重力卫星数据的普及，为重力勘探、深度地质构造、壳幔变形以及月球火星等星体的探索提供了可能。

重力场正反演作为一项重要的技术手段，被广泛运用于地下三维密度成像研究。重力场正反演可以在常规的勘探尺度进行，用于资源勘探(矿产、石油、天然气)和水文环境等研究，大多基于直角坐标系下进行，常见的快速算法包括fft法，gauss-fft法，多级展开法等，使得反演的计算效率和精度得到了保证。然而当我们研究某一区域乃至全球尺度下的重力异常时，就必须考虑到地球曲率的影响，传统的直角坐标系下重力场正反演方法已不再适用。其相应的数据处理、正演算法和反演都将在球坐标系下进行。

由于球面三维积分没有解析解，必须通过数值方法进行离散求解，常见的数值方法包括二维高斯勒让德积分法、三维高斯勒让德积分法、基于泰勒级数展开法等。这些传统方法计算效率低下，很难适应大规模三维反演的要求，正演计算效率将直接制约着大尺度高分辨率重力场反演。在当今计算机计算性能难以有质的飞跃的背景下，想要大幅度提高计算效率，只能寻求新的正演方法。其次，传统三维重力场反演中，大型雅克比矩阵的内存占用是一个严重制约反演分辨率的问题。

另一方面，21世纪以来，卫星重力测量技术的发展不仅革命性地推动了地球重力场的研究，也推动了利用地球重力场信息探测地球内部组构及物质迁移的研究，以及利用grail卫星探测月球内部质量瘤盆地及月幔隆起的研究。日益丰富的卫星重力(包括卫星测高)、海洋重力、陆地重力及gps测量数据使地球重力场的解算精度和分辨率得到极大提高。全球覆盖的高分辨率、高精度重力数据为获取区域及全球壳幔高分辨密度扰动模型创造了条件，也是固体地球科学研究重要的基础数据，同时为全球大尺度大规模三维卫星重力反演提供了新的契机。传统的重力场正反演方法已难以使用海量高分辨率数据的基本要求，亟需开发一套高效高精度正演和反演方法。

技术实现要素：

本发明提供一种基于toplitze核矩阵(托普利兹核矩阵)的重力场快速正演方法，使得重力场核矩阵与密度向量相乘的复杂计算转换到频域，使得计算效率提高了三个数量级，同时减少了大型稀疏矩阵的存储，有利于大尺度高分辨率重力场反演，具体技术方案如下：

一种基于toplitze核矩阵的重力场快速正演方法，包括以下步骤：

步骤一、将地下三维场源在经度方向、纬度方向和深度方向分别剖分成等间隔的nλ、和nr段，共计个tesseroid单元体，剖分间距分别为δλ、和δr；观测面位于三维场源正上方，且观测点位置分布与三维场源剖分tesseroid单元体中心点位置一一对应，观测点个数共计

步骤二、地下三维场源tesseroid单元体在观测点上所产生的重力异常的矩阵-向量乘积采用表达式1)表示：

k·ρ＝g1)；

其中：k是重力异常核矩阵，ρ是地下tesseroid单元体的密度向量，g是观测点上重力异常向量；

步骤三、获得重力异常核矩阵k的子矩阵、密度向量ρ的子向量以及观测重力异常向量g的子向量；

获得重力异常核矩阵k的子矩阵具体是：

将核矩阵k划分成个子矩阵，每个子矩阵的大小是nλ×nλ，如表达式4)：

其中：核矩阵k的每一个子矩阵ki,j都是托普利兹矩阵，如表达式5)所示，

定义ki,j＝toeplitz(t)，表示ki,j是托普利兹矩阵，t＝(t1-n,…,t-1,t0,t1,…,tn-1)是该托普利兹矩阵的特征向量，且有ta＝t-a，a＝0,1,…,n-1，n＝nλ是该特征向量t的元素的个数；

将密度向量ρ和观测重力异常向量g剖分成子向量形式，如表达式6)：

其中：密度向量ρ的每一个子向量ρj和观测重力异常向量g的每一个子向量gi的长度均为nλ，且和

步骤四、找出与核矩阵k的子矩阵ki,j对应的密度向量ρ的子向量ρj和观测重力异常向量g的子向量gi，如表达式7)所示：

则该核矩阵k的子矩阵ki,j与相应的密度向量ρ的子向量ρj的乘积写成表达式8)：

ki,jρj＝gi8)；

步骤五、调用快速傅里叶变换算法进行计算并得到观测重力异常向量g的子向量gi，具体是：

通过调用快速傅里叶变换算法来实现子矩阵ki,j与密度向量ρ的子向量ρj的卷积运算如表达式11)：

其中：是核矩阵k的子矩阵ki,j的扩展矩阵，大小是2n×2n；t^ext是向量t的扩展向量，长度是2n；和是观测重力异常向量g的子向量gi和密度向量ρ的子向量ρj的扩展向量，长度都是2n，其前n项元素与gi和ρj相同，后n项元素为0，0代表零向量，长度为n；s为一中间变量矩阵，用于存储临时计算结果，大小是n×n；和分别表示一维傅里叶逆变换和正变换；

将取其前n项得到观测重力异常向量g的子向量gi；

步骤六、取i＝i+1或j＝j+1，若i小于等于或j小于等于返回步骤四；获得与核矩阵k的每一个子矩阵ki,j对应的gi；

步骤七、将步骤六获得的gi按表达式6)的顺序排列得到观测重力异常向量g。

以上技术方案中优选的，所述步骤三中获得重力异常核矩阵k的子矩阵先要计算得到重力异常核矩阵k，具体是：

采用表达式2)计算重力异常核矩阵k：

其中：r0是编号为p和q的观测点的径向分量，是观测点的纬度分量，λp是观测点的经度分量，p＝1,2,…,nλ，rl′、和λ′m表示编号为l、n和m的tesseroid单元体中心点坐标三分量，l＝1,2,…,nr，m＝1，g是重力常数；r′、和λ′分别是径向、纬度方向和经度方向上的积分变量，积分区间分别为rl′-0.5δr至rl′+0.5δr、至和λ′m-0.5δλ至λ′m+0.5δλ；为观测点到tesseroid单元体内任一点的距离，cosψ为观测点到tesseroid单元体内任一点的中间角度的余弦，采用表达式3)计算得到：

以上技术方案中优选的，所述步骤四中：

由于托普利兹矩阵ki,j的对称性，该矩阵和向量的乘积运算是一个一维的离散卷积运算，其元素gk有表达式9)：

其中：k＝0,1,…,n-1；t是托普利兹矩阵ki,j的特征向量，tk-b是向量t的第k-b个元素；ρj是密度向量ρ的子向量，ρb是密度子向量ρj的第b个元素，j＝0,1,…,n-1，b＝0,1,…,n-1；n是该特征向量t的元素个数；

将托普利兹矩阵ki,j构造成循环矩阵其中n×n循环矩阵c＝circulant(c)具有表达式10)；

其中：c＝(c0,c1,…,cn-1)表示循环矩阵c的特征向量；

将托普利兹矩阵ki,j转换成一个2n×2n的循环矩阵t^ext＝(t0,t1,…,tn-1,0,t1-n,…,t-1)表示循环矩阵的特征向量。

本发明还提供一种基于toplitze核矩阵的重力场快速反演方法，具体包括如下步骤：

获取球坐标系下反演目标函数φ(ρ)如表达式13)：

其中：wd是一个观测数据权重的对角矩阵，每一个对角线元素代表该观测点数据的可信度；wρ是一个大型稀疏的模型权重矩阵；β是正则化系数，用于权衡模型目标函数与数据目标函数之间的比重；ρ^ref表示参考模型；k是重力异常核矩阵，ρ是地下tesseroid单元体的密度向量；

通过最小化反演目标函数φ(ρ)得到如下线性方程组如表达式14)：

其中，k^t是k的转置，wd^t是wd的转置，wρ^t是wρ的转置；

将上述正演方法所得观测重力异常向量g代入表达式14)获得地下tesseroid单元体的密度向量ρ。

应用本发明的技术方案，具有以下有益效果：

1、本发明提出了新的球坐标系下三维重力场快速正演方法，通过将正演核矩阵剖分成多个toplitze子矩阵(托普利兹子矩阵)，然后分块采用频域fft快速计算，最后累加求和，将耗时的矩阵-向量乘积运算转换为求和运算，将正演中核矩阵计算时间减少了2个数量级，将核矩阵与密度向量相乘的重复性计算效率提高了20多倍，正演计算效率整体提高3个数量级。

2、在计算重力场核矩阵时，由于采用特殊剖分方式，无需再计算所有核矩阵元素，仅需计算地下场源第一列tesseroid单元体所产生的核矩阵元素，其它元素的值可以通过等效关系映射出来。这一方法大大降低了核矩阵计算成本，同时降少了内存占用。在反演中，同样无需存储大型雅克比矩阵，仅需存储最终矩阵(βwρ^twρρ和βwρ^twρρ^ref)即可，减少了中间变量的计算及存储。

3、为了证明本发明所提出方法的正确性，本发明给出的一个正演模型算例，模型包含360×720×10个单元体，观测点数位360×720，本发明所提出的正演方法用时1965秒，而传统方法的计算时间是890964秒，计算效率提高432倍。此外，本发明正演方法和反演方法的最大相对误差都是0.024％，证明了本发明所提出方法的有效性。

然后将本正演方法和反演方法运用于月球质量瘤的反演，反演结果和前人得出的莫霍面深度信息基本一致，进一步证明了本发明方法的正确性。

本发明所提出的正演和反演方法既适用于勘探尺度的资源勘查，水文环境监测，又适用于大尺度三维密度成像、研究地球壳幔结构、构造运动及岩石圈变迁等，应用范围广泛。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是球坐标系下tesseroid单元体剖分示意图；

图2是一层3×4个tesseroid单元体与观测点的关系以及核矩阵等效性关系示意图，其中：(a)为核矩阵及其子矩阵；(b)和(c)为核矩阵元素等效性示意图，其中相同线型所表示的核矩阵元素值相等；

图3是本发明所提出方法与传统方法的最大相对误差随剖分间隔的关系图；

图4是本发明所提出方法与传统方法的相对计算时间随剖分间隔的关系图；

图5是本发明优选实施例用于处理月球全球重力数据的过程图片，其中：(a)月球全球地形图；(b)自由空气重力异常；(c)地形改正；(d)布格重力异常，(d)中英文单词为月球上大型质量瘤盆地的名称，且黑色实线所表示的6条断面；

图6是各个深度上的密度分布剖面图，其中：(a)为5km，(b)为15km，(c)为25km，(d)为35km，(e)为45km，(f)为55km，(g)为65km，(h)为75km，(i)为85km，(j)为95km；

图7是图5中(d)上的6条断面的反演密度分布，编号为1-6；图中黑色加粗实线表示前人得出的莫霍面深度。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

本发明提出一种球坐标系下重力场快速正演方法(即核矩阵与密度向量乘积计算方法)，并将该方法运用于重力场三维密度反演，详情如下：

基于频域fft算法的球坐标系下重力场快速正演方法，具体是：

第一步、将地下三维场源剖分成多个tesseroid单元体，详情如下：

对勘探尺度而言，重力场正反演常在直角坐标系下进行，然而当研究区域较大或者研究整个地球岩石圈时，地球曲率的影响不得不考虑，因此，需要在球坐标系下进行。在球坐标系下，场源通常被剖分为如图1所示的tesseroid单元体。将地下三维场源在经度方向、纬度方向和深度方向分别剖分成等间隔的nλ、和nr段，共计个tesseroid单元体，剖分间距分别为δλ、和δr。观测面位于三维场源正上方，且观测点中心点位置与三维场源剖分tesseroid单元体中心点位置一一对应，观测点个数共计个，如图2所示。

第二步、获得地下三维场源tesseroid单元体在观测点上所产生的重力异常的矩阵-向量乘积，具体采用表达式1)表示：

k·ρ＝g1)；

其中，k是重力异常核矩阵，其维度是nobs×ntess；ρ是地下tesseroid单元体的密度向量，其维度是ntess；g是观测点上重力异常向量，其维度是nobs。

采用表达式2)计算重力异常核矩阵k：

其中：r0是编号为p和q的观测点的径向分量，是观测点的纬度分量，λp是观测点的经度分量，p＝1,2,…,nλ，rl′、和λ′m表示编号为l、n和m的tesseroid单元体中心点坐标三分量，l＝1,2,…,nr，m＝1，g是重力常数；r′、和λ′分别是径向、纬度方向和经度方向上的积分变量，积分区间分别为rl′-0.5δr至rl′+0.5δr、至和λ′m-0.5δλ至λ′m+0.5δλ；为观测点到tesseroid单元体内任一点的距离，cosψ为观测点到tesseroid单元体内任一点的中间角度的余弦，采用表达式3)计算得到：

理论公式表达式2)没有解析解，必须通过数值方法求解，可以选择二维高斯勒让德方法、三维高斯勒让德方法，或者泰勒级数展开法，由于泰勒级数展开法在极点处存在解析奇点，二维高斯勒让德方法计算公式复杂且耗时，在本发明选择三个方向上都是2点的三维高斯勒让德数值方法求解上述核矩阵。高斯点数越多，计算精度越高，但相应的计算时间也成倍增加。选择2点高斯系数，可满足高精度正演需求，最大相对误差低于0.1％。

第三步、获得重力异常核矩阵k的子矩阵、密度向量ρ的子向量以及观测重力异常向量g的子向量，详情是：

从图2可以看出，核矩阵k是一个大型分块的toplitze矩阵，也就是说核矩阵里有大量相等的元素存在。因此在核矩阵只计算时，只需对于所有p＝1,2,…,nλ，l＝1,2,…,nr，m＝1，计算。这里m只取1，表示只需计算第一列tesseroid单元体所产生的核矩阵元素，以减少重复计算。通过该策略，减少了nλ倍的计算量以及内存占用量，在大尺度三维正反演中，该策略的优势更加显著。

将计算好的核矩阵k划分成个子矩阵，每个子矩阵的大小是nλ×nλ，采用表达式4)表示：

其中：核矩阵k的每一个子矩阵ki,j都是托普利兹矩阵，如表达式5)，

定义ki,j＝toeplitz(t)，表示ki,j是托普利兹矩阵，t＝(t1-n,…,t-1,t0,t1,…,tn-1)是该托普利兹矩阵的特征向量，且有ta＝t-a，对于a＝0,1,…,n-1，且n＝nλ是该特征向量t的元素的个数；

将密度向量ρ和观测重力异常向量g剖分成子向量如表达式6)：

其中：密度向量ρ的每一个子向量ρj和观测重力异常向量g的每一个子向量gi的长度均为nλ，且和

第四步、找出与核矩阵k的子矩阵ki,j对应的密度向量ρ的子向量ρj和观测重力异常向量g的子向量gi，如表达式7)：

则核矩阵子矩阵ki,j与密度向量子向量ρj的乘积写成表达式8)：

ki,jρj＝gi8)。

由于托普利兹矩阵ki,j的对称性，该矩阵和向量的乘积运算是一个一维的离散卷积运算，其元素gk有表达式9)：

其中：k＝0,1,…,n-1；t是托普利兹矩阵ki,j的特征向量，tk-b是向量t的第k-b个元素；ρj是密度向量ρ的子向量，ρb是密度子向量ρj的第b个元素，j＝0,1,…,n-1，b＝0,1,…,n-1；n是该特征向量t的元素个数；

将托普利兹矩阵ki,j构造成循环矩阵其中n×n循环矩阵c＝circulant(c)为表达式10)：

其中：c＝(c0,c1,…,cn-1)表示循环矩阵c的特征向量；

将托普利兹矩阵ki,j转换成一个2n×2n的循环矩阵t^ext＝(t0,t1,…,tn-1,0,t1-n,…,t-1)表示循环矩阵的特征向量。

第五步、通过调用快速傅里叶变换算法来实现子矩阵ki,j与密度向量ρ的子向量ρj的卷积运算如表达式11)：

其中：是核矩阵k的子矩阵ki,j的扩展矩阵，大小是2n×2n；向量t^ext是向量t的扩展向量，长度是2n；和分别是观测重力异常向量g的子向量gi和密度向量ρ的子向量ρj的扩展向量，长度都是2n，其前n项元素与gi和ρj相同，后n项元素为0；0代表零向量，长度为n；s为一中间变量矩阵，用于存储临时计算结果，大小是n×n；和分别表示一维傅里叶逆变换和正变换；

将取其前n项得到子向量gi。

至此本实施例已经将子矩阵与密度向量的乘积通过频域fft算法快速计算，接下来只需计算每一个子矩阵与相对应的密度子向量的乘积(具体是：完成第四步至第五步后，取i＝i+1或j＝j+1，若i小于等于或j小于等于返回第四步)，然后按表达式6)的顺序排列得到观测重力异常向量g。计算每一个核矩阵子矩阵与相对应密度子向量的乘积，并累加求和，得到重力异常子向量gi，采用下面表达式12)进行运算：

基于频域fft算法的球坐标系下重力场快速反演方法，具体是：

运用吉洪诺夫正则化反演方法，构造球坐标系下反演目标函数φ(ρ)如表达式13)：

其中：wd是一个观测数据权重的对角矩阵，每一个对角线元素代表该观测点数据的可信度；wρ是一个大型稀疏的模型权重矩阵；β是正则化系数，用于权衡模型目标函数与数据目标函数之间的比重；ρ^ref表示参考模型，一般取ρ^ref＝0；k是重力异常核矩阵，ρ是地下tesseroid单元体的密度向量。

通过最小化该反演目标函数可得到如表达式14)的线性方程组：

其中：k^t是k的转置，wd^t是wd的转置，wρ^t是wρ的转置。本发明使用成熟的共轭梯度法求解该方程组。

值得注意的是，上述的基于托普利兹(toplitze)核矩阵频域fft算法的球坐标系下重力场快速正演方法，可以直接用于计算表达式14)中的k·ρ。由于k^t同样也是一个toplitze矩阵，因此可以先计算wd^twdk·ρ和wd^twdg，得到两向量，然后利用上述基于托普利兹(toplitze)核矩阵频域fft算法的球坐标系下重力场快速正演方法，计算k^t这两个向量的乘积。并且无需存储大型矩阵wρ^t和wρ，只需存储最终矩阵βwρ^twρρ和βwρ^twρρ^ref，让于大尺度高分辨率反演成为可能。

在利用上述基于toplitze矩阵的球坐标系下快速反演方法时，具体实施步骤为：首先从网格化grid文件中输入观测数据，程序自动判断观测数据个数，起始范围，并计算出剖分间隔，以及每一个测点的坐标信息等；然后设置观测高度，设置观测数据误差均方差；接着设置反演参数，程序需要输入反演深度的起始范围，该范围可根据该地区的地质、区域构造、前人研究结果得到；最后运行程序，按照前述公式进行运算，得到最终反演结果。

实施例1：

本实施例使用一个月球球壳正演模型来验证本发明所提出方法的正确性与高效性，因为只有球壳模型才存在解析解。该球壳层厚度为100km，起始范围分别为1638到1738km(即月球实际半径)，球壳密度设置为1000kg/m³，观测面度为距离外球壳10km高度的球面，将该球壳层在径向剖分为10段，间隔10km。在经向和纬向上的间隔为0.5°到10°，统计最大相对误差和计算时间随剖分间隔的关系。

最大相对误差随剖分间隔关系如图3所示，相对计算时间随剖分间隔关系如图4所示。为了更加直观感受本发明所提出方法的高效性，表1统计了剖分间隔为0.5°(即剖分单元体个数为360×720×10)时的最大相对误差、核矩阵计算时间、核矩阵与密度向量乘积计算时间以及总时间的对比。

在这个正演测试中，剖分单元体个数为360×720×10时，基于本发明所提出的方法的计算时间为1965秒，约0.5小时，而用传统的正演方法，计算时间为890964秒，约248小时，计算效率提升约500倍。此外本发明所提出的方法同样可以保证高精度，其最大相对误差与传统方法相同，低于0.1％。

表1最大相对误差、核矩阵计算时间、核矩阵与密度向量乘积计算时间以及总时间的对比表

实施例2：

为了进一步证明本发明申请所提出的基于toplitze矩阵频域fft算法的球坐标系重力场快速正演方法的有效性和高效性，接下来本实施例给出一个实际数据反演案例，将所提出方法运用于月球全球尺度三维密度成像研究，详情如下：

月球全球地形是以及最新的月球地形球谐模型lro_ltm05_2050解算得到，如图5中(a)所示。

本发明使用月球最新重力场模型gl1500e来计算自由空气重力异常，如图5中(b)所示，地形改正结果如图5中(c)所示。布格重力异常通过自由空气重力异常(图5中(b))减去地形改正(图5中(c))得到，如图5中(d)所示。在这里，地形校正所用的密度是2560kg/m³，观测高度是10km。用模型的前180阶，对应于空间分辨率为1.0°。

本发明将地下场源剖分成180×360个等间隔的tesseroid单元体，根据前人研究基础，在深度方向上模型范围设置为0km到100km，剖分成10段，间隔为10km，共计180×360×10个tesseroid单元体，观测点数为180×360个。用该布格重力异常做三维重力场反演，各个深度异常密度分布如图6所示。沿着图5中(d)中所示的6条断面的密度分布如图7所示。可以看出，本发明所提出的方法能正确反演出异常体所在位置和深度，高密度异常体主要分布在各个质量瘤盆地的底部，其在深度上的分布范围大概是15-60km。此外，这些高密度异常体的顶界面与前人研究的莫霍面起伏深度一致。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。