CN113109883B - 基于等参变换全球离散网格球坐标下卫星重力场正演方法 - Google Patents
基于等参变换全球离散网格球坐标下卫星重力场正演方法 Download PDFInfo
- Publication number
- CN113109883B CN113109883B CN202110580088.9A CN202110580088A CN113109883B CN 113109883 B CN113109883 B CN 113109883B CN 202110580088 A CN202110580088 A CN 202110580088A CN 113109883 B CN113109883 B CN 113109883B
- Authority
- CN
- China
- Prior art keywords
- integral
- hexagonal prism
- gravity
- transformation
- grid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
- G01V7/02—Details
- G01V7/06—Analysis or interpretation of gravimetric records
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于等参变换全球离散网格球坐标下卫星重力场正演方法,首先选择一定经纬度范围作为观测网范围,并确定观测高度和地球物理模型,利用等面积投影方法,离散地球物理模型,从而获得一系列dggrid网格;此外,基于六边形棱柱体积分公式,计算获得dggrid网格插值形函数、积分点和积分权系数,继而利用等参变换计算所有dggrid网格于各观测点的异常响应。相比于传统经纬度网格的非均匀性和极点奇异性缺点(在极点附近四边形会退化成三角形),dggrid网格具备一致相邻性特点。与任意四面体解析解对比,基于等参变换全球离散网格(dggrid)的球坐标下卫星重力场正演精度较高。
Description
技术领域
本发明属于地球物理勘探技术领域,具体涉及基于等参变换全球离散网格球坐标下卫星重力场正演方法。
背景技术
快速高效正演计算是大规模位场反演的基础。当前,重力及其梯度正演方法的有效性、时效性是制约着大规模重力资料反演解释的重要原因。对于局部地区的重力资料而言,可将重力观测面近似为平面,其数据处理问题可在直角坐标系中进行理论推导和计算。然而,对于区域性的,乃至全球尺度问题,重力观测面已不能近似为平面,而是一个球面问题,需要利用球坐标进行表达与处理。于球坐标系下,计算重力位方法可大致分为两类:基于频率域的球谐函数方法和基于空间域的直接积分方法。对于前者,球谐函数方法在中低阶(N<360)相对于数值积分方法计算效率高,但在(超)高阶(N>~2700)会受到数值不稳定性的影响。且球谐函数方法不能适当地考虑给定质量分布的垂直延伸,实际计算结果幅值稍大,不适用于精确地全局或局部重力场建模。
大尺度条件下,地球表面是一个曲面,其反演模型空间的离散化形式也与笛卡尔坐标系下有差别。在球坐标系下,物性反演沿纬度、经度及径向等三个方向分别按照等间隔方式划分为一个个单元块体,这类单元称之为经纬度网格的Tesseroid单元体(Anderson,1976)。经纬度网格的最大优点是网格编码过程简单,且与三维空间坐标的转换易于实现。但随着纬度的增加,从赤道到两极会产生退化现象,在极点附近四边形会退化成三角形,产生大量数据冗余,因此所占内存比较多。
基于以上描述,亟需一种卫星重力场正演方法,以解决现有技术计算过程中产生大量冗余数据,且计算结果不精确的问题。
发明内容
本发明目的在于提供一种基于等参变换全球离散网格球坐标下卫星重力场正演方法,通过该方法可快速获得更高精度的重力、重力矢量及重力梯度张量正演计算结果。
为解决上述技术问题,本发明的一种基于等参变换全球离散网格球坐标下卫星重力场正演方法的具体技术方案如下:
一种基于等参变换全球离散网格球坐标下卫星重力场正演方法,包括以下步骤:
S1、构造具有n个节点的正六边形棱柱;
S2、确定正六边形棱柱在局部坐标系中的空间位置和积分点个数;
S3、构建待分析形函数N:
S4、将n节点正六边形棱柱的顶点坐标代入形函数N,构建如下矩阵形式方程组:
其中,(xj,yj,zj)为第j个正六边形棱柱的顶点坐标,
S5、获得各积分点的积分权重;
S6、获得dggrid网格(全球离散网格Discrete Global Grid)插值形函数、积分点和积分权系数;
S7、利用等参变换计算所有dggrid网格与各观测点的异常响应。
作为优选,在步骤S1中,正六边形棱柱是由任意n个节点六边形棱柱等参变换而来。
作为优选,在步骤S2中,正六边形棱柱在局部坐标系中的空间位置和积分点个数是根据六边形棱柱的对称性确定的。
作为优选,在步骤S3中,构建待分析形函数N的方法如下:
利用帕斯卡三角形来确定一由坐标(x,y,z)表示待用多项式列表V:
V={1,x,y,z,xy,…,x5…,y5…,z5…},
作为优选,在步骤S5中,各积分点的积分权重是根据等参变换中雅可比矩阵的定义计算得来的,该矩阵的行列式等于六边形棱柱的体积,即利用各形函数Ni基于六边形棱柱的体积分等于六边形棱柱的体积,构建线性方程组获得。
作为优选,在步骤S7中,异常响应计算方法如下:
球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式表达:
其中,G为万有引力常数;Ω地质体所处的积分域;dm与dv为积分质量元与体积元;r为观测点与积分单元的距离;u(x,y,z)为观测点处的引力位;ρ(x',y',z')为地质体的密度或异常密度分布;
根据上述的积分公式对径向r求导,即得到重力或重力异常:
通过u对x和y求偏导得到重力、重力向量或重力梯度张量的各个分量。
有益效果:
本发明充分利用等参变换可将系统坐标系中复杂形体转换为局部坐标系中规则形体的优点,将球壳六边形棱柱的积分区间通过等参变换转换至正六边形棱柱的积分区间,具有以下优点:
1.球壳六边形棱柱具有极强的对称性,这为核函数的重复使用提供了可能,可极大的降低正演核矩阵内存占用。
2.就物性网格离散而言,各球壳六边形棱柱间具有一致相邻性,这避免了传统tesseroid网格因网格位置的不同,而tesseroid网格大小存在显著差异的问题。
3.避免实现传统tesseroid网格积分时需要获得积分根的问题,从而加快正演计算速度。
4.通过本发明提供的计算流程,易于获得更高次积分型函数和积分点,从而获得更高精度的重力、重力矢量及重力梯度张量正演计算结果。
附图说明
图1为本发明的方法流程图;
图2为本发明用于对比验证的四面体剖分立方体示意图;
图3为本发明用于对比验证的立方体正演获得重力矢量和重力梯度张量结果;
图4为本发明用于对比验证的四面体剖分立方体后正演获得重力矢量和重力梯度张量结果;
图5为本发明用于对比验证的立方体及四面体剖分立方体后正演获得重力矢量和重力梯度张量结果的差;
图6为本发明的全球离散网格示意图(线框图);
图7为本发明的全球离散网格示意图(渲染图);
图8为本发明的观测示意图及四面体剖分单一全球离散网格示意图;
图9为本发明的用于对比验证四面体剖分单一全球离散网格重力梯度张量分量结果对比图;
图10为本发明的用于对比验证四面体剖分单一全球离散网格重力梯度矢量分量结果对比图。
具体实施方式
为了更好地了解本发明的方法及步骤,下面结合附图,对本发明一种基于等参变换全球离散网格球坐标下卫星重力场正演方法做进一步详细的描述。
首先,引入任意四面体网格正演计算技术,以验证本发明所提正演方法的有效性和可靠性,在此之前,需要验证任意四面体网格正演计算方法的准确性。为此,将与具有精确解析解的长方体单元作为验证任意四面体网格单元正演计算精度的参考,可以分如下几步实现该目的。第一步,以任意四面体网格单元为例,通过三维Hammer积分,获得其重力矢量和重力梯度张量的解析解;第二步,如图2所示,利用Delaunay三角化技术,将长方体单元离散为多个任意四面体单元;第三步,分别计算上述两种正演方法的重力矢量和重力梯度张量值,通过对比图3-图5,可以发现基于三维Hammer积分重力矢量/重力梯度张量解析解具有很高的计算精度。
其次,在此基础上,利用任意四面体网格正演计算方法验证本发明所提正演方法的有效性和可靠性,涉及如下几个步骤:
步骤一,如图6和图7所示,确定研究区域球壳厚度H,沿着球壳半径方向将其等分为nz层,根据全球离散网格定义,在球面上构建其中的第一级网格(亦可以进一步在球面上构建其中的第n级网格),其中五边形全球离散网格12个,六边形全球离散网格30个,并利用这些球面上的网格,沿着半径方向逐层离散研究区域球壳,即每一层球壳被离散为五边形棱柱的全球离散网格12个,六边形棱柱的全球离散网格30个,这里将其称之为全球离散网格球壳单元。
步骤二,获得本发明基于等参变换全球离散网格球坐标下卫星重力场正演计算公式,具体方法如下:
球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式为:
其中,G为万有引力常数;Ω地质体所处的积分域;dm与dv分别为积分质量元与体积元;r为观测点与积分单元的距离;u(x,y,z)为观测点处的引力位;ρ(x',y',z')为地质体的密度或异常密度分布。
根据上述的积分公式对径向r求导,即得到重力或重力异常:
通过u对x和y求偏导得到重力、重力向量或重力梯度张量的各个分量。
为实现上述公式的积分,首先选择一定经纬度范围作为观测网范围,并确定观测高度和地球物理模型,在等参变换下,构造自n个节点六边形棱柱至n个节点正六边形棱柱的投影。
根据六边形棱柱的对称性确定等参变换的与局部坐标系中的空间位置和积分点个数;
利用帕斯卡三角形确定一由坐标(x,y,z)表示待用多项式列表V:
V={1,x,y,z,xy,…,x5…,y5…,z5…};
将n节点正六边形棱柱的顶点坐标代入形函数,构建如下矩阵形式方程组:
其中,(xj,yj,zj)为第j个正六边形棱柱的顶点坐标,
根据等参变换中雅可比矩阵的定义,该矩阵的行列式等于六边形棱柱的体积,即利用各形函数Ni基于六边形棱柱的体积分等于六边形棱柱的体积,构建线性方程组,以获得各积分点的积分权重。
基于六边形棱柱体积分公式获得dggrid网格(全球离散网格Discrete GlobalGrid)插值形函数、积分点和积分权系数。
步骤三,选择任意一个或多个全球离散网格球壳单元,利用随机数生成函数,对所选全球离散网格球壳单元,赋予随机剩余密度值,以模拟地壳中的复杂地质构造,利用等参变换计算所有dggrid网格于各观测点的重力矢量和重力梯度张量异常响应。
再者,如图8所示,利用Delaunay三角化技术,将上述对所选全球离散网格球壳单元,离散一系列的任意四面体网格单元,进而利用基于三维Hammer积分重力矢量和重力梯度张量解析解,计算上述对所选全球离散网格球壳单元的重力矢量和重力梯度张量异常响应。
可选地,在Delaunay三角化过程中,可限制四面体网格的奇强度不小于1.4,以保证正演计算精度。
最后,通过对比图9和图10,可以发现本发明的重力矢量/重力梯度张量正演结果具有很高的计算精度。
可以理解,本领域技术人员知悉的,在不脱离本发明的精神和范围的情况下,可以对该方法进行各种改变或等效替换。另外,在本发明的教导下,可以对该方法进行修改以适应具体的情况而不会脱离本发明的精神和范围。因此,本发明不受此处所公开的具体实施例的限制,所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。
Claims (6)
1.一种基于等参变换全球离散网格球坐标下卫星重力场正演方法,其特征在于,包括以下步骤:
S1、构造具有n个节点的正六边形棱柱;
S2、确定正六边形棱柱在局部坐标系中的空间位置和积分点个数;
S3、构建待分析形函数N:
S4、将n节点正六边形棱柱的顶点坐标代入形函数N,构建如下矩阵形式方程组:
其中,(xj,yj,zj)为第j个正六边形棱柱的顶点坐标,
S5、获得各积分点的积分权重;
S6、获得dggrid网格(全球离散网格Discrete Global Grid)插值形函数、积分点和积分权系数;
S7、利用等参变换计算所有dggrid网格与各观测点的异常响应。
2.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法,其特征在于,在步骤S1中,正六边形棱柱是由任意n个节点六边形棱柱等参变换而来。
3.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法,其特征在于,在步骤S2中,正六边形棱柱在局部坐标系中的空间位置和积分点个数是根据六边形棱柱的对称性确定的。
5.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法,其特征在于,在步骤S5中,各积分点的积分权重是根据等参变换中雅可比矩阵的定义计算得来的,该矩阵的行列式等于六边形棱柱的体积,即利用各形函数Ni基于六边形棱柱的体积分等于六边形棱柱的体积,构建线性方程组获得。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110580088.9A CN113109883B (zh) | 2021-05-26 | 2021-05-26 | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110580088.9A CN113109883B (zh) | 2021-05-26 | 2021-05-26 | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113109883A CN113109883A (zh) | 2021-07-13 |
CN113109883B true CN113109883B (zh) | 2023-01-03 |
Family
ID=76723243
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110580088.9A Active CN113109883B (zh) | 2021-05-26 | 2021-05-26 | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113109883B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113570724B (zh) * | 2021-07-29 | 2023-12-08 | 辽宁大学 | 基于逆球极投影的网格球面保角参数化方法及其应用 |
CN113960690B (zh) * | 2021-09-03 | 2023-05-05 | 中国人民解放军战略支援部队信息工程大学 | 一种海面重力数据测量精度对海底地形反演结果影响计算方法及装置 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106646645A (zh) * | 2016-12-29 | 2017-05-10 | 中南大学 | 一种新的重力正演加速方法 |
CN107561592A (zh) * | 2017-09-11 | 2018-01-09 | 中南大学 | 一种密度为多项式的任意多面体的重力场正演方法 |
CN109375280A (zh) * | 2018-12-10 | 2019-02-22 | 中南大学 | 一种球坐标系下重力场快速高精度正演方法 |
CN109444973A (zh) * | 2018-11-06 | 2019-03-08 | 中南大学 | 一种球坐标系下重力正演加速方法 |
CN110045432A (zh) * | 2018-12-05 | 2019-07-23 | 中南大学 | 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 |
CN111400654A (zh) * | 2020-03-13 | 2020-07-10 | 中南大学 | 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法 |
-
2021
- 2021-05-26 CN CN202110580088.9A patent/CN113109883B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106646645A (zh) * | 2016-12-29 | 2017-05-10 | 中南大学 | 一种新的重力正演加速方法 |
CN107561592A (zh) * | 2017-09-11 | 2018-01-09 | 中南大学 | 一种密度为多项式的任意多面体的重力场正演方法 |
CN109444973A (zh) * | 2018-11-06 | 2019-03-08 | 中南大学 | 一种球坐标系下重力正演加速方法 |
CN110045432A (zh) * | 2018-12-05 | 2019-07-23 | 中南大学 | 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 |
CN109375280A (zh) * | 2018-12-10 | 2019-02-22 | 中南大学 | 一种球坐标系下重力场快速高精度正演方法 |
CN111400654A (zh) * | 2020-03-13 | 2020-07-10 | 中南大学 | 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法 |
Non-Patent Citations (4)
Title |
---|
3D Inversion of Gravity Data with Lanczos Bidiagonalization and Unstructured Mesh;K. Danaei,et al.;《EAGE 2020 Annual Conference & Exhibition Online》;20201231;全文 * |
基于单元细分H-自适应有限元全张量重力梯度正演;曹书锦等;《地球物理学进展》;20100630;第25卷(第3期);全文 * |
球冠体积分的重力异常正演方法及其与Tesseroid单元体泰勒级数展开方法的比较;杜劲松等;《测绘学报》;20120615(第03期);全文 * |
球坐标系密度界面反演方法及在华南大陆的应用;王祥等;《物探与化探》;20200929(第05期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN113109883A (zh) | 2021-07-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113109883B (zh) | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 | |
CN110045432B (zh) | 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 | |
CN109375280B (zh) | 一种球坐标系下重力场快速高精度正演方法 | |
CN106646645B (zh) | 一种重力正演加速方法 | |
Zhang et al. | Meshless analysis of potential problems in three dimensions with the hybrid boundary node method | |
CN113311494B (zh) | 一种卫星重力场反演方法 | |
Heikes et al. | Optimized icosahedral grids: Performance of finite-difference operators and multigrid solver | |
CN106683185B (zh) | 一种基于大数据的高精度曲面建模方法 | |
CN112528546B (zh) | 一种非结构化网格的重磁数据三维正反演方法 | |
CN108984939A (zh) | 基于3D Gauss-FFT的三维重力场正演方法 | |
CN109636912A (zh) | 应用于三维声呐图像重构的四面体剖分有限元插值方法 | |
CN112526625B (zh) | 航空重力测量点的布格重力异常值的计算装置 | |
CN113671570B (zh) | 一种地震面波走时和重力异常联合反演方法与系统 | |
CN113239573B (zh) | 基于无网格波动建模的封闭空间声场重构方法 | |
EP2442246B1 (en) | Method of solving the non-uniformly discretized Poisson equation | |
CN111158059A (zh) | 基于三次b样条函数的重力反演方法 | |
CN113516754A (zh) | 一种基于磁异常模量数据的三维可视化成像方法 | |
CN117233858A (zh) | 基于聚类约束的隧道多源数据交叉梯度联合反演方法 | |
CN112748471B (zh) | 一种非结构化等效源的重磁数据延拓与转换方法 | |
CN117725354B (zh) | 一种大数据量重力与重力梯度联合的快速正反演方法及系统 | |
Yagawa | Free Mesh Method: fundamental conception, algorithms and accuracy study | |
Katz et al. | Discretization methodology for high aspect ratio prismatic grids | |
Surendranath et al. | Upwind Interpolation near 1-D Hanging Node Interfaces for Compressible Euler Equation Building-Cube Method Simulations | |
Dutt | High order stochastic transport and Lagrangian data assimilation | |
Surendranath et al. | Building-Cube Method Simulations Using Upwind Interpolations for Ghost Cells Near 2-D Hanging Node Interfaces |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |