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

本发明专利技术提供一种球坐标系下三维重力场快速正演方法，通过将正演核矩阵剖分成多个Toplitze子矩阵，然后分块采用频域FFT快速计算，最后累加求和，将耗时的矩阵‑向量乘积运算转换为求和运算，减少正演中核矩阵的计算时间，提高核矩阵与密度向量相乘的重复性计算效率，提高正演计算的整体效率；本发明专利技术采用特殊剖分方式，无需再计算所有核矩阵元素，仅需计算地下场源第一列tesseroid单元体所产生的核矩阵元素，其它元素的值通过等效关系映射出来，大大降低了核矩阵的计算成本，同时也降少了内存的占用。本发明专利技术还提供一种反演方法，结合上述正演方法，无需存储大型雅克比矩阵，仅需存储最终矩阵，减少了中间变量的计算及存储。

【技术实现步骤摘要】
一种基于Toplitze核矩阵的重力场快速正演方法及反演方法
本专利技术涉及地球物理及勘探
，具体涉及一种基于Toplitze核矩阵的重力场快速正演方法及反演方法。
技术介绍
随着高精度高覆盖率全球重力卫星数据的普及，为重力勘探、深度地质构造、壳幔变形以及月球火星等星体的探索提供了可能。重力场正反演作为一项重要的技术手段，被广泛运用于地下三维密度成像研究。重力场正反演可以在常规的勘探尺度进行，用于资源勘探(矿产、石油、天然气)和水文环境等研究，大多基于直角坐标系下进行，常见的快速算法包括FFT法，Gauss-FFT法，多级展开法等，使得反演的计算效率和精度得到了保证。然而当我们研究某一区域乃至全球尺度下的重力异常时，就必须考虑到地球曲率的影响，传统的直角坐标系下重力场正反演方法已不再适用。其相应的数据处理、正演算法和反演都将在球坐标系下进行。由于球面三维积分没有解析解，必须通过数值方法进行离散求解，常见的数值方法包括二维高斯勒让德积分法、三维高斯勒让德积分法、基于泰勒级数展开法等。这些传统方法计算效率低下，很难适应大规模三本文档来自技高网...

【技术保护点】
1.一种基于Toplitze核矩阵的重力场快速正演方法，其特征在于，包括以下步骤：/n步骤一、将地下三维场源在经度方向、纬度方向和深度方向分别剖分成等间隔的N

【技术特征摘要】
1.一种基于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；text是向量t的扩展向量，长度是2n；giext和ρjext分别是观测重力异常向量g的子向量gi和密度向量ρ的子向量ρj的扩展向量，长度都是2n，其前n项元素与gi和ρj相同，后n项元素为0，0代表零向量，长度为n；S为一中间变量矩阵，用于存储临时计算结果，大小是n×n；和分别表示一维傅里叶逆变换和正变换；
将giext取其前n项得到观测重力异常向量g的子向量gi；
步骤六、取i＝i+1或j＝j+1，若i小于等于或j小于等于返回步骤四；获得与核矩阵K的每一个子矩阵ki,j对应的gi；
步骤七、将步骤六获得的gi按表达...

【专利技术属性】
技术研发人员：赵广东，柳建新，陈波，
申请(专利权)人：中南大学，
类型：发明
国别省市：湖南;43