3D inversion of gravity data by separation of sources and the method of local corrections: Kolarovo gravity high case study

We present a novel methodology for 3D gravity/magnetic data inversion. It combines two algorithms for preliminary separation of sources and an original approach to 3D inverse problem solution. The first algorithm is designed to separate sources in depth and to remove the shallow ones. It is based on subsequent upward and downward data continuation. For separation in the lateral sense, we approximate the given observed data by the field of several 3D line segments. For potential field data inversion we apply a new method of local corrections. The method is efficient and does not require trial-and-error forward modeling. It allows retrieving unknown 3D geometry of anomalous objects in terms of restricted bodies of arbitrary shape and contact surfaces. For restricted objects, we apply new integral equations of gravity and magnetic inverse problems. All steps of our methodology are demonstrated on the Kolarovo gravity anomaly in the Danube Basin of Slovakia.     

海域航空重力快速构建区域大地水准面

我国在海域开展了大规模的航空重力勘探,这些资料对构建高精度大地水准面具有重要价值.基于此,本文提出一种利用海域航空重力测量数据快速构建大地水准面的方法.该方法基于移去-恢复法思想,利用位场最小曲率方法对航空重力数据进行高精度向下延拓并获取相应的扰动位,实现航空重力测量快速构建海域大地水准面.与斯托克斯积分计算相比,采用了处理效率更高的频率域位场转换,解决了向下延拓及垂向积分时航空重力异常数据空白及扩边问题,具有较高的位场转换精度.本文应用EGM2008模拟航空重力数据进行模型验证,计算结果与其给出的水准面的精度相当;同时,也选取GRAV-D计划的航空重力数据进行实际验证,计算结果与xGEOID18B水准面模型精度基本一致.模型验证和实际应用验证了本方法的实用性.     

位场向下延拓的波数域广义逆算法

位场向下延拓是位场数据处理和反演中的重要运算，但是它的不稳定性影响了它在许多处理和反演方法技术中的应用.本文通过把位场向下延拓视为向上延拓的反问题，得到向下延拓的褶积型线性积分方程，再利用Fourier变换矩阵的正交对称特性，并结合矩阵的奇异值分解和广义逆原理，提出了一种稳定的不需要进行求逆运算的位场向下延拓广义逆方法——波数域广义逆算法，解决了位场大深度向下延拓的不稳定性问题.把这种方法用于三维理论模型数据和实际磁场数据的向下延拓获得了理想的结果.     

A geoid solution for airborne gravity data

Airborne gravity data is usually attached with satellite positioning of data points, which allow for the direct determination of the gravity disturbance at flight level. Assuming a suitable gridding of such data, Hotine’s modified integral formula can be combined with an Earth Gravity Model for the computation of the disturbing potential (T) at flight level. Based on T and the gravity disturbance data, we directly downward continue T to the geoid, and we present the final solution for the geoid height, including topographic corrections. It can be proved that the Taylor expansion of T converges if the flight level is at least twice the height of the topography, and the terrain potential will not contribute to the topographic correction. Hence, the simple topographic bias of the Bouguer shell yields the only topographic correction. Some numerical results demonstrate the technique used for downward continuation and topographic correction.     

Optimal Model for Geoid Determination from Airborne Gravity

Two different approaches for transformation of airborne gravity disturbances, derived from gravity observations at low-elevation flying platforms, into geoidal undulations are formulated, tested and discussed in this contribution. Their mathematical models are based on Green's integral equations. They are in these two approaches defined at two different levels and also applied in a mutually reversed order. While one of these approaches corresponds to the classical method commonly applied in processing of ground gravity data, the other approach represents a new method for processing of gravity data in geoid determination that is unique to airborne gravimetry. Although theoretically equivalent in the continuous sense, both approaches are tested numerically for possible numerical advantages, especially due to the inverse of discretized Fredholm's integral equation of the first kind applied on different data. High-frequency synthetic gravity data burdened by the 2-mGal random noise, that are expected from current airborne gravity systems, are used for numerical testing. The results show that both approaches can deliver for the given data a comparable cm-level accuracy of the geoidal undulations. The new approach has, however, significantly higher computational efficiency. It would be thus recommended for real life geoid computations. Additional errors related to regularization of gravity data and the geoid, and to accuracy of the reference field, that would further deteriorate the quality of estimated geoidal undulations, are not considered in this study.     

Atmospheric effects on satellite gravity gradiometry data

Atmospheric masses play an important role in precise downward continuation and validation of satellite gravity gradiometry data. In this paper we present two alternative ways to formulate the atmospheric potential. Two density models for the atmosphere are proposed and used to formulate the external and internal atmospheric potentials in spherical harmonics. Based on the derived harmonic coefficients, the direct atmospheric effects on the satellite gravity gradiometry data are investigated and presented in the orbital frame over Fennoscandia. The formulas of the indirect atmospheric effects on gravity anomaly and geoid (downward continued quantities) are also derived using the proposed density models. The numerical results show that the atmospheric effect can only be significant for precise validation or inversion of the GOCE gradiometric data at the mE level.     

Spectral Combination of Spherical Gradiometric Boundary-Value Problems: A Theoretical Study

The Earth’s gravity potential can be determined from its second-order partial derivatives using the spherical gradiometric boundary-value problems which have three integral solutions. The problem of merging these solutions by spectral combination is the main subject of this paper. Integral estimators of biased- and unbiased-types are presented for recovering the disturbing gravity potential from gravity gradients. It is shown that only kernels of the biased-type integral estimators are suitable for simultaneous downward continuation and combination of gravity gradients. Numerical results show insignificant practical difference between the biased and unbiased estimators at sea level and the contribution of far-zone gravity gradients remains significant for integration. These contributions depend on the noise level of the gravity gradients at higher levels than sea. In the cases of combining the gravity gradients, contaminated with Gaussian noise, at sea and 250?km levels the errors of the estimated geoid heights are about 10 and 3 times smaller than those obtained by each integral.     

位场向下延拓的最小曲率方法

针对位场向下延拓的不适定性,我们将位场向下延拓视为向上延拓的反问题,提出以位场最小曲率作为约束条件来求解稳定的下延位场.我们将剖面位场向上延拓表达式用傅里叶矩阵的形式表示,以矩阵乘法形式给出延拓的表达式,同时向待反演的下延位场引入最小曲率约束,得到向下延拓的最小曲率解,并利用正交变换给出了更为简洁的频率域解.随后,利用Kronecker积将上述全部结果拓展至三维位场,给出了三维位场向下延拓的最小曲率解.此外,我们将位场数据的填充、扩充问题与向下延拓问题统筹考虑,提出一种新的向下延拓迭代格式,该算法面向实际资料处理需求、无须预扩充或填补数据.下延迭代时,对原始数据直接向下延拓,而空白部分利用上一次下延位场估计的上延值替代其空白值并对其向下延拓,直至获得最小曲率约束下稳定的向下延拓结果.同时,我们也讨论了利用改进L曲线和广义交叉验证(GCV)计算正则参数最优估计的问题.对理论模型和实际航空重力资料进行了向下延拓检验,处理结果表明位场向下延拓的最小曲率方法解能满足实际位场资料对向下延拓处理的需求,具有较高的下延精度.     

Does Poisson’s downward continuation give physically meaningful results?

The downward continuation (DWC) of the gravity anomalies from the Earth’s surface to the geoid is still probably the most problematic step in the precise geoid determination. It is this step that motivates the quasi-geoid users to opt for Molodenskij’s rather than Stokes’s theory. The reason for this is that the DWC is perceived as suffering from two major flaws: first, a physically meaningful DWC technique requires the knowledge of the irregular topographical density; second, the Poisson DWC, which is the only physically meaningful technique we know, presents itself mathematically in the form of Fredholm integral equation of the 1st kind. As Fredholm integral equations are often numerically ill-conditioned, this makes some people believe that the DWC problem is physically ill-posed. According to a revered French mathematician Hadamard, the DWC problem is physically well-posed and as such gives always a finite and unique solution. The necessity of knowing the topographical density is, of course, a real problem but one that is being solved with an ever increasing accuracy; so sooner or later it will allow us to determine the geoid with the centimetre accuracy.     

Numerical Studies on the Harmonic Downward Continuation of Band-Limited Airborne Gravity

In this paper, the numerical stability and efficiency of methods of harmonic downward continuation from flying altitudes are treated for sampled gravity field data. The problem is first formulated in its continuous form, i.e. as the inverse solution of the spherical Dirichlet problem, and is then approximated by Gaussian quadrature to yield a finite system of linear equations. The numerical stability of this system is investigated for both error-free gravity data and for the noisy and band-limited gravity measurements usually obtained from airborne gravity surveys. It can be shown that the system becomes ill-conditioned, once the ratio between flying altitude and data sampling rate exceeds a certain limit. It can also be shown that noisy measurements tend to generate a solution that is practically useless, long before the system becomes ill-conditioned. Therefore, instead of treating the general solution of the discrete downward continuation problem, the more modest question is studied, for which range of flying altitudes and sampling rates, the numerical solution of the discrete linear system can be considered as practically useful. Practically useful will be defined heuristically as of sufficient accuracy and stability to satisfy the requirements of the user. The question will be investigated for the specific application of geoid computation from gravity data sampled at flying altitudes. In this case, a stable solution with a standard deviation of a few centimeters is required. Typical flight parameters are heights of 2–6 km, a minimum half-wavelength resolution of 2 km, and data noise between 0.5 and 1.5 mGal. Different methods of geoid determination, different solution techniques for the resulting systems of linear equations, and different minimization principles will be compared. As a result operational parameters will be defined which, for a given noise level, will result in a geoid accuracy of a few centimeters for the estimated band-limited gravity field spectrum.     

Regional gravity field modeling based on rectangular harmonic analysis

Regional gravity field modeling with high-precision and high-resolution is one of the most important scientific objectives in geodesy,and can provide fundamental information for geophysics,geodynamics,seismology,and mineral exploration.Rectangular harmonic analysis(RHA)is proposed for regional gravity field modeling in this paper.By solving the Laplace’s equation of gravitational potential in local Cartesian coordinate system,the rectangular harmonic expansions of disturbing potential,gravity anomaly,gravity disturbance,geoid undulation and deflection of the vertical are derived,and so are the formula for signal degree variance and error degree variance of the rectangular harmonic coefficients(RHC).We also present the mathematical model and detailed algorithm for the solution of RHC using RHA from gravity observations.In order to reduce the edge effects caused by periodic continuation in RHA,we propose the strategy of extending the size of computation domain.The RHA-based modeling method is validated by conducting numerical experiments based on simulated ground and airborne gravity data that are generated from geopotential model EGM2008 and contaminated by Gauss white noise with standard deviation of 2 mGal.The accuracy of the 2.5′×2.5′geoid undulations computed from ground and airborne gravity data is 1 and 1.4cm,respectively.The standard error of the gravity disturbances that downward continued from the flight height of 4 km to the geoid is only 3.1 mGal.Numerical results confirm that RHA is able to provide a reliable and accurate regional gravity field model,which may be a new option for the representation of the fine structure of regional gravity field.     

Numerical behaviour of the downward continuation of gravity anomalies

The numerical results of downward continuation (DWC) of point and mean gravity anomalies by the Poisson integral using point, single mean, and doubly averaged kernel are examined. Correct evaluation of the integral in its innermost zone is a challenging task. To avoid instabilities, an analytical planar approximation is used in the innermost integration zone. In addition it is shown that the single mean mode has the minimum discretization error. Downward continuation of point and mean anomalies by singly and doubly averaged kernel are the same mean anomalies on the geoid.     

位场数据曲化平的迭代法

位场数据曲化平是位场数据处理解释中的重要运算,但是它的计算量和计算的复杂性影响了它在许多处理和解释方法技术中的应用.本文提出一种位场数据曲化平的迭代方法,即通过把位场数据曲化平视为平面位场数据向上延拓的反问题,得到曲化平的线性积分方程,再把曲面上位场数据视为曲面平均高程面上的位场数据,利用向下延拓的波数域广义逆算法把平均高程面上的位场数据向下延拓到设定平面上,再根据曲面和其平均高程面的相对起伏对设定平面上的向下延拓数据进行起伏校正,最后再把所得平面上的位场数据向上延拓得到曲面上的位场数据,并进行迭代.把这种方法用于三维理论模型数据和实际磁场数据的曲化平处理均获得了理想的结果.     

基于重力卫星数据监测地表质量变化的三维点质量模型法

利用卫星重力测量手段监测全球质量变化取得了巨大成功,本文基于牛顿万有引力定律在三维空间直角坐标系中导出利用重力卫星观测数据监测全球质量变化的三维点质量模型法,该方法可直接利用重力卫星的轨道和星间观测数据或时变重力场模型计算全球质量变化,由于利用卫星观测数据计算地表质量变化的向下延拓过程以及观测数据噪声的影响,需要采用合适的空间约束方程或正则化技术对解算结果进行约束或平滑处理.利用合成全球质量变化模型模拟一个月的GRACE双星轨道和星间距离变率数据计算全球质量变化,对三维点质量模型法进行分析验证,采用零阶Tikhonov正则化技术处理病态问题.结果表明,三维点质量模型法可有效用于重力卫星观测数据监测全球质量变化,为利用重力卫星观测数据监测全球质量变化提供一种可选的途径.     

积分迭代法的正则性分析及其最优步长的选择

位场积分迭代法是一种解决位场大数据量、大深度向下延拓的有效方法.本文基于Kirsch正则化子理论,推导了积分迭代法对应的正则化滤子函数,并证明积分迭代法为一种求解位场向下延拓不适定反问题的正则化方法.针对积分迭代法迭代步长固定、迭代次数较多影响收敛速度的问题,提出该迭代法最优迭代步长的选择原理.理论模型和实测数据对比分...     

补偿向下延拓方法研究及应用

位场向下延拓是重、磁处理和解释的常用方法,但其不稳定性限制了其在资料处理及反演中的应用.本文基于补偿圆滑滤波思想以及空间域向下延拓迭代法,通过逐次补偿的办法实现位场的稳定向下延拓.同时,在频率域空间给出了该下延方法的频率域响应因子,并讨论了其低通滤波特性,理论模型和实际位场资料试验表明该方法向下延拓稳定性具有较高的延拓精度.将其应用于重力密度界面反演中,改进反演的稳定性,实际莫霍界面反演表明下延因子具备实用性.     

WAVE SCATTERING DECONVOLUTION BY SEISMIC INVERSION*

We propose a wave scattering approach to the problem of deconvolution by the inversion of the reflection seismogram. Rather than using the least-squares approach, we study the full wave solution of the one-dimensional wave equation for deconvolution. Randomness of the reflectivity is not a necessary assumption in this method. Both the reflectivity and the section multiple train can be predicted from the boundary data (the reflection seismogram). This is in contrast to the usual statistical approach in which reflectivity is unpredictable and random, and the section multiple train is the only predictable component of the seismogram. The proposed scattering approach also differs from Claerbout's method based on the Kunetz equation. The coupled first-order hyperbolic wave equations have been obtained from the equation of motion and the law of elasticity. These equations have been transformed in terms of characteristics. A finite-difference numerical scheme for the downward continuation of the free-surface reflection seismogram has been developed. The discrete causal solutions for forward and inverse problems have been obtained. The computer algorithm recursively solves for the pressure and particle velocity response and the impedance log. The method accomplishes deconvolution and impedance log reconstruction. We have tested the method by computer model experiments and obtained satisfactory results using noise-free synthetic data. Further study is recommended for the method's application to real data.     

Geoid determination by spectral combination of an Earth gravitational model with airborne and terrestrial gravimetry — a theoretical study

Today air-gravimetry is a versatile technique to quickly collect gravity data over large regions, where terrestrial gravity data are sparse and/or of poor quality. The method requires the data to be downward continued to sea level for use in geoid determination, an inverse problem operation that calls for smoothing of the data and/or the kernel function involved (in either spectral or space domain). In this purely theoretical study we avoid this separate computational step by performing direct geoid estimation by so-called spectral combination/filtering of the data, which includes terrestrial gravimetry, airgravimetry, an Earth Gravitational Model (EGM) as well as their signal and error degree variances. Each derived geoid estimator is presented as the sum of one or two integral formulas and the harmonic series of the EGM together with the expected mean square error of the estimator. The article is limited to a theoretical study, leaving its practical tests for future investigation.     

重力场向下延拓Milne法

重力场向下延拓能够突出局部和浅部的异常信息,分离叠加的异常特征.但是向下延拓通常具有过程不稳定、下延深度小、结果不准确等问题.针对向下延拓所存在的不足,本文利用重力场及其垂向一阶导数,基于辛普森(Simpson)求积公式,推导出重力场向下延拓米尔尼(Milne)公式.将本文向下延拓方法应用于模型数据,向下延拓模型结果及误差曲线表明,相对于向下延拓快速傅里叶变换(FFT)法和积分迭代法,向下延拓Milne法的深度更大,相对误差更小;相对模型值,向下延拓Milne法能够获得稳定且准确的结果.对加拿大乃查科(Nechako)盆地地区实测航空重力数据进行本文方法向下延拓验证,处理结果表明,相对于实测异常,本文方法向下延拓结果能够很好还原实测数据,并且在进一步向下延拓中反映原始异常的趋势,增强局部和细小异常信息.     

位场向下延拓的波数域迭代法及其收敛性

提出了位场向下延拓的波数域迭代法. 对水平面上的位场观测值进行Fourier变换,得到其波谱. 根据第一类Fredholm积分方程的空间域迭代解法,推导出计算向下延拓水平面上位场波谱的波数域迭代公式. 在波数域中进行迭代,一直进行到相继两次迭代近似解的差值最大绝对值小于给定的精度,或迭代达到给定的最大迭代次数. 对这种迭代近似解进行Fourier逆变换,得到向下延拓的位场. 数值计算结果表明：与空间域迭代法比较,这种波数域迭代法简单、快速,并有同样好的向下延拓效果. 本文还证明了这种迭代法是收敛的,并给出了它的收敛特性和滤波特性.