Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction

A finite-element method for computing the electric field in a 3-D conductivity model of the Earth for plane wave sources, thus enabling magnetotelluric responses to be calculated, is presented. The method incorporates in the iterative solution of the electric-field system of equations the divergence correction technique introduced for finite-difference solutions by Smith (1996). The correction technique accelerates the development of the discontinuity of the normal component of the approximate electric field across conductivity discontinuities. The convergence rate of the iterative solution is improved significantly, especially for low frequencies. The correction technique involves computing the divergence of the current density for the approximate electric field, computing the static potential whose source is this divergence of the current density, and ‘correcting’ the approximate electric field by subtracting from it the gradient of the potential. This is repeated at regular intervals during the iterative solution of the electric-field system of equations. For the method presented here, the Earth model is discretised using a rectilinear mesh comprising uniform cells. Edge-element basis functions are used to approximate the electric field and nodal basis functions are used to approximate the correction potential. The Galerkin method is used to derive the systems of equations for the approximate electric field and correction potential from the respective differential equations. A bi-conjugate gradient solver was found to be adequate for the system of equations for the correction potential; a generalised minimum residual solver was found to be better for the electric-field system of equations. The method is illustrated using the COMMEMI 3D-1A and 3D-2A models.     

基于高斯牛顿法的频率域可控源电磁三维反演研究

三维反演解释是电磁法勘探发展的重要趋势,而如何提高三维反演的可靠性、稳定性和计算效率是算法开发者们目前的研究重点.本文实现了一种频率域可控源电磁(CSEM)三维反演算法.其中正演基于拟态有限体积法离散化,利用直接矩阵分解技术来求解大型线性系统方程,不仅准确、稳定,而且特别有利于含有大量发射场源位置的CSEM勘探情况;对目标函数的最优化采用高斯牛顿法(GN),具有近似二次的收敛性;使用预条件共轭梯度法(PCG)求解每次GN迭代所得到的法方程,避免了显式求解和存储灵敏度矩阵,减小了计算量.以上这些方法的结合应用,使得本文的三维反演算法准确、稳定且高效.通过陆地和海洋CSEM勘探场景中的典型理论模型的反演测试,验证了本文算法的有效性.     

An airborne tensor VLF system. From concept to realization1

Radio signals from very low frequency (VLF) transmitters distributed world-wide have been used for several decades to study the lateral variations of the electrical conductivity in the upper few hundred metres of the earth's crust. Traditionally, in airborne applications, the total magnetic fields from one or two transmitters are measured to form the basis for construction of maps that primarily show those conductive structures that are parallel or subparallel to the direction to the transmitters. The tensor VLF technique described in this paper makes use of all signals available in a predefined frequency band to construct transfer functions relating the vertical magnetic field and the two horizontal magnetic field components. These transfer functions are uniquely determined for a particular measuring site and contain information about the lateral conductivity variations in all directions. First experiences with real field data, acquired during a test survey in Sweden, show that maps of the so-called peaker, the spatial divergence of the transfer functions, give an image of the conducting structures. Most of the structures can be correlated to small valleys filled with conducting sediments or valleys underlain by conductive fracture zones in the crystalline rocks.     

三维陆地可控源电磁法有限元模型降阶快速正演

三维陆地可控源电磁法有限元快速正演的主要瓶颈在于多频率大型稀疏方程组求解问题.本文引入一种基于模型降阶的Krylov子空间投影算法,推导了有限元刚度矩阵的模型降阶形式,构建了频率域传递函数;采用标准正交向量序列,构建一个远远小于有限元刚度矩阵维度的矩阵,该矩阵与频率无关,通过一次模型降阶即可实现多频点有限元方程快速求解.采用基于电场的变分方程,加入散度校正条件,以消除伪解;引入伪δ函数,消除了源点的奇异性,可适用于复杂背景模型三维有限元数值模拟,并为多源的求解奠定了基础;以层状介质模型解析解为标准,通过和基于Pardiso直接求解器的有限元算法（3DFEM）进行比较,模型降阶法计算时间小于前者的1/10,平均相对误差在1.72%,在满足精度要求下,实现了高效率三维有限元数值求解;分别设计了横向高低阻模型和纵向高低阻模型,分析了从近区到远区电场和卡尼亚视电阻率的变化规律,假极值的表现特征,阴影效应的影响等,从而也验证了该算法的正确性.最后,建立了一个地层陷落柱模型,通过模型降阶有限元正演模拟,发现视电阻率断面图在陷落柱上方出现"凹陷",与模型设计吻合,表明该算法对复杂地层模拟具有同样的适用性.     

An efficient multivariate random field generator using the fast Fourier transform

It is often convenient to use synthetically generated random fields to study the hydrologic effects of spatial heterogeneity. Although there are many ways to produce such fields, spectral techniques are particularly attractive because they are fast and conceptually straightforward. This paper describes a spectral algorithm for generating sets of random fields which are correlated with one another. The algorithm is based on a discrete version of the Fourier-Stieltjes representation for multidimensional random fields. The Fourier increment used in this representation depends on a random phase angle process and a complex-valued spectral factor matrix which can be readily derived from a specified set of cross-spectral densities (or cross-covariances). The inverse Fourier transform of the Fourier increment is a complex random field with real and imaginary parts which each have the desired coveriance structure. Our complex-valued spectral formulation provides an especially convenient way to generate a set of random fields which all depend on a single underlying (independent) field, provided that the fields in question can be related by space-invariant linear transformations. We illustrate this by generating multi-dimensional mass conservative groundwater velocity fields which can be used to simulate solute transport through heterogeneous anisotropic porous media.     

2.5D forward modeling and inversion of frequency-domain airborne electromagnetic data

Frequency-domain airborne electromagnetics is a proven geophysical exploration method. Presently, the interpretation is mainly based on resistivity—depth imaging and one-dimensional layered inversion; nevertheless, it is difficult to obtain satisfactory results for two- or three-dimensional complex earth structures using 1D methods. 3D forward modeling and inversion can be used but are hampered by computational limitations because of the large number of data. Thus, we developed a 2.5D frequency-domain airborne electromagnetic forward modeling and inversion algorithm. To eliminate the source singularities in the numerical simulations, we split the fields into primary and secondary fields. The primary fields are calculated using homogeneous or layered models with analytical solutions, and the secondary (scattered) fields are solved by the finite-element method. The linear system of equations is solved by using the large-scale sparse matrix parallel direct solver, which greatly improves the computational efficiency. The inversion algorithm was based on damping least-squares and singular value decomposition and combined the pseudo forward modeling and reciprocity principle to compute the Jacobian matrix. Synthetic and field data were used to test the effectiveness of the proposed method.     

New induction sounding tested in Central Europe

The hourly data of nine geomagnetic observatories situated in Central Europe have been analyzed using the generalized magnetovariation (GMV) method designed recently for induction soundings of inhomogeneous media. In this method, impedance is one of transfer functions in the differential relation between spectra of the magnetic components and their derivatives. The peculiarity of this impedance is its correspondence to the magnetotelluric one estimated from the linear relations. Three transfer functions have been estimated simultaneously for data of geomagnetic observatories, using three different routines working in the period range from three hours up to two days. Noises in the source field components have been compared with noise in the estimated plane field divergence. The multivariate errors-in-variables method was used to extract spatially and temporally coherent geomagnetic field structure from the partially incoherent geomagnetic variations. This method allows estimating reliably impedances and gradient tippers for each observatory, taking into consideration the Earth’s sphericity. The obtained responses have been used for induction soundings and for detecting a deep inhomogeneity in the region.     

用于波场成像的谱法LU分解

地震波场模拟和偏移成像等有限差分隐格式算法中的重要环节，是实现亥姆霍兹算子表示矩阵H的快速求逆运算. 在螺旋边界条件下，H具有Toeplitz结构的正定厄密矩阵，其快速求逆可由谱法LU分解实现. 本文分析了谱法LU分解对提高计算速度的原理及特点，并着重讨论了在不同类型的介质模型中，采用谱法分解矩阵H时带来的数值误差、误差的分布及其对波场计算的影响. 研究结果表明，对均匀介质而言，矩阵H各列具有相同的非零元素分布，谱法LU分解的误差在吸收边界条件下，不影响波场模拟和成像计算；但对于非均匀介质模型，矩阵H各列具有不同的非零元素分布，谱法LU分解的误差随介质不均匀性程度的增大而增大，势必影响非均匀介质中波场计算. 在波场模拟和成像等有限差分隐格式算法中，采用谱法LU分解完成矩阵求逆时，必须考虑到并尽量减少该方法的误差对波场计算的影响.     

基于非均一场源的球坐标大地电磁模拟方法

大地电磁测深理论与数据处理解释均假定平面电磁波垂直入射地下空间,但随着研究尺度的逐渐扩大,使得因地球弧度产生的影响难以忽略.此时,传统笛卡尔坐标体系及平面波场源不再适用于大尺度的大地电磁数据正反演解释.本文提出并实现了一种基于球坐标系的大地电磁交错网格有限差分三维正演算法,并对电场进行极向-环向分解,结合球谐函数和贝塞尔函数构建了可取代平面波的场源模型.首先利用经纬度信息构建三维地电模型,将场源设置于模型空间正上方,然后通过直接求解球坐标系下麦克斯韦方程来获得大地电磁响应.在此基础上,本文设计了球坐标下具有不同分辨率的多个三维地电模型,阐述了由球体模型到笛卡尔模型的转换方法,详细对比了两种坐标体系在计算效率、所求得的电场和视电阻率方面的差异.结果表明二者差异度主要与电性横向分布和地图投影方法有关,与周期并不存在明显的单调递增关系.     

三维三分量CSAMT法有限元正演模拟研究初探

首先从麦克斯韦方程出发,用伽里金方法推导了三维三分量CSAMT法的有限元方程.在研究过程中,认识到加入散度条件的必要性,在公式中强加了散度条件,提高了解的完备性.其次将成功应用于二维线源频率域电磁法有限元模拟中的两种技术推广到三维中,一是边界条件统一采用一阶吸收边界,使线源产生的电磁波在边界上按波的传播规律被吸收,以降低平面波假设造成的影响;二是总体系数矩阵的存储,用两个二维数组分别记录总体系数矩阵的非零元素及其在总体结点编号中所处的位置,使总体系数矩阵的存储量达到最小的同时,物理意义明确,迭代求解时迅速简便.最后用均匀半空间模型进行了验证.     

Seasonal and interannual behaviour of groundwater catchment boundaries in a Chalk aquifer

Groundwater catchment boundaries and their associated groundwater catchment areas are typically assumed to be fixed on a seasonal basis. We investigated whether this was true for a highly permeable carbonate aquifer in England, the Berkshire and Marlborough Downs Chalk aquifer, using both borehole hydrograph data and a physics‐based distributed regional groundwater model. Borehole hydrograph data time series were used to construct a monthly interpolated water table surface, from which was then derived a monthly groundwater catchment boundary. Results from field data showed that the mean annual variation in groundwater catchment area was about 20% of the mean groundwater catchment area, but interannual variation can be very large, with the largest estimated catchment size being approximately 80% greater than the smallest. The flow in the river was also dependent on the groundwater catchment area. Model results corroborated those based on field data. These findings have significant implications for issues such as definition of source protection zones, recharge estimates based on water balance calculations and integrated conceptual modelling of surface water and groundwater systems. Copyright © 2015 John Wiley & Sons, Ltd.     

利用共轭梯度法的电阻率三维反演研究

利用共轭梯度（ＣＧ）迭代技术,实现了直流电阻率测量数据的三维最小构造反演。 首先,运用共轭梯度迭代算法解反问题的线性方程组,只需求Ｊａｃｏｂｉａｎ矩阵Ｇ与任一向量ｘ的 乘积Ｇｘ及Ｇ^Ｔ与任何一向量ｙ的乘积Ｇ^Ｔｙ,再引入Ｇ的Ｒｏｄｉ算法,则Ｇｘ及Ｇ^Ｔｙ均可在每次反 演迭代中的一次正演计算后一并求得．因此,每次反演迭代仅需一次正演计算,大大加快了 计算速度;而且避免了直接求 Ｇ以及 Ｇ^ＴＧ的逆矩阵,也避免了存储 Ｇ和 Ｇ^ＴＧ所需庞大的存储 量。另外,由于反演参数太多,求模型光滑约束的最小构造反演能够有效地消除多余构造信 息,得到可靠的反演结果。将这３种方法和技术融合于三维反演中,取得了好的反演效果．为 改进传统最小构造反演收敛慢的问题,还提出了一种新的反演迭代技术,仅需１０次左右甚至 更少迭代即可收敛。     

Unified multi‐depth‐level field decomposition

Wavefield decomposition forms an important ingredient of various geophysical methods. An example of wavefield decomposition is the decomposition into upgoing and downgoing wavefields and simultaneous decomposition into different wave/field types. The multi‐component field decomposition scheme makes use of the recordings of different field quantities (such as particle velocity and pressure). In practice, different recordings can be obscured by different sensor characteristics, requiring calibration with an unknown calibration factor. Not all field quantities required for multi‐component field decomposition might be available, or they can suffer from different noise levels. The multi‐depth‐level decomposition approach makes use of field quantities recorded at multiple depth levels, e.g., two horizontal boreholes closely separated from each other, a combination of a single receiver array combined with free‐surface boundary conditions, or acquisition geometries with a high‐density of vertical boreholes. We theoretically describe the multi‐depth‐level decomposition approach in a unified form, showing that it can be applied to different kinds of fields in dissipative, inhomogeneous, anisotropic media, e.g., acoustic, electromagnetic, elastodynamic, poroelastic, and seismoelectric fields. We express the one‐way fields at one depth level in terms of the observed fields at multiple depth levels, using extrapolation operators that are dependent on the medium parameters between the two depth levels. Lateral invariance at the depth level of decomposition allows us to carry out the multi‐depth‐level decomposition in the horizontal wavenumber–frequency domain. We illustrate the multi‐depth‐level decomposition scheme using two synthetic elastodynamic examples. The first example uses particle velocity recordings at two depth levels, whereas the second example combines recordings at one depth level with the Dirichlet free‐surface boundary condition of zero traction. Comparison with multi‐component decomposed fields shows a perfect match in both amplitude and phase for both cases. The multi‐depth‐level decomposition scheme is fully customizable to the desired acquisition geometry. The decomposition problem is in principle an inverse problem. Notches may occur at certain frequencies, causing the multi‐depth‐level composition matrix to become uninvertible, requiring additional notch filters. We can add multi‐depth‐level free‐surface boundary conditions as extra equations to the multi‐component composition matrix, thereby overdetermining this inverse problem. The combined multi‐component–multi‐depth‐level decomposition on a land data set clearly shows improvements in the decomposition results, compared with the performance of the multi‐component decomposition scheme.     

基于二次场电导率分块连续变化的三维可控源电磁有限元数值模拟

从电偶源三维地电断面可控源电磁法的二次电场边值问题及其变分问题出发,采用任意六面体单元对研究区域进行剖分,并且在单元分析中同时对电导率及二次电场进行三线性插值,实现电导率分块连续变化情况下,基于二次场的可控源电磁三维有限元数值模拟.这个新的可控源电磁三维正演方法可以模拟实际勘探中地下任意形状及电性参数连续变化的复杂模型.理论模型的计算结果表明,均匀大地计算的视电阻率误差和相位误差分别为0.002%和0.0005°.分层连续变化模型的有限元计算结果表明,其与对应的分层均匀模型解析结果有明显差异.三维异常体组合模型以及倾斜异常体等复杂模型的有限元计算结果也有效地反映了异常形态.     

Smoothing seismically derived velocities1

It is often very useful to be able to smooth velocity fields estimated from exploration seismic data. For example seismic migration is most successful when accurate but also smooth migration velocity fields are used. Smoothing in one, two and three dimensions is examined using North Sea velocity data. A number of ways for carrying out this smoothing are examined, and the technique of locally weighted regression (LOESS) emerges as most satisfactory. In this method each smoothed value is formed using a local regression on a neighbourhood of points downweighted according to their distance from the point of interest. In addition the method incorporates ‘blending’ which saves computations by using function and derivative information, and ‘weighting and robustness’ which allows the smooth to be biased towards reliable points, or away from unreliable ones. A number of other important factors are also considered: namely, the effect of changing the scales of axes, or of thinning the velocity field, prior to smoothing, as well as the problem of smoothing on to irregular subsurfaces.     

基于有限体积法的二维大地电磁各向异性数值模拟

为了计算带任意地形的各向异性介质中二维大地电磁响应,本文在非结构化网格的基础上,采用有限体积法,开发了二维大地电磁各向异性正演模拟的新算法.首先,从Maxwell方程出发,推导二维各向异性介质中大地电磁场的边值问题;然后,采用三角网格自动生成技术对求解区域进行非结构化网格剖分,进而构建节点中心控制体积单元,利用有限体积方法,得到求解边值问题的大型稀疏线性方程组;最后,利用Pardiso精确地计算了大地电磁响应值.三个各向异性模型的计算结果表明,本文开发的有限体积算法,不仅能够高精度求解带任意地形的大地电磁电导率各向异性问题,而且对于同一模型,该方法的计算消耗和精度都与有限单元法相当.因此,有限体积法是处理电磁法各向异性问题的一种有效方法.     

三维地形频率域人工源电磁场的边界元模拟方法

提出了一种用边界元法计算频率域人工源三维地形电磁场的数值模拟方法.首先用矢量积分理论和电磁场边界条件，将上半空间（空气）和下半空间（地下介质）两个区域电磁场边值问题变为仅对地形界面的两个矢量面积分方程.然后将对地形界面的积分剖分为一系列的三角单元积分.在三角单元积分中，假设单元中电磁场为无限大气空间电磁场与地形影响的叠加,并假设地形影响为常项,这样既保证了计算精度又使得计算方法简便.通过分解和计算，每一个矢量面积分方程分解为对应三个坐标方向的三个常量线性方程，这些线性方程组成了对角占优的线性方程组，可用SSOR方法求解.文中给出了垂直磁偶源的垂直磁场地形影响的例子.     

Analytical configurations of a force-free magnetic cylinder in the solar wind

The problem of a smooth field configuration, which should be an initial configuration in modeling (using the method of coarse particles) the problem of a stationary solar wind flow around a magnetic cloud in the case of a spatially two-dimensional statement (when a magnetic cloud is considered as a force-free magnetic cylinder with a finite radius) is considered. It has been indicated that such a statement is possible only when the magnetic field in the solar wind is parallel to the cylinder axis. The method for finding the magnetic field of a force-free cylinder with a finite radius, when some field component is specified and another component is determined based on this one (which makes it possible to construct fields with preassigned properties), has been proposed. The variant for constructing the initial field configuration in the transition region around a cylinder has been proposed. This variant makes it possible to gradually pass from homogeneous crossed fields in the solar wind to a force-free magnetic and zero electric fields within a cylinder, an electric field being potential and orthogonal to a magnetic field (in the reference system related to a magnetic cloud).     

瞬变电磁2.5维反演中灵敏度矩阵计算方法(Ⅰ)

一种具体的反演方法要分三个重要的组成部分：1)反演算法;2)正演模拟;3)灵敏度矩阵计算．本文给出了快速计算瞬变电磁2．5维反演问题中灵敏度的简单推导．只需求解两个基本相同的边值问题(原始场和辅助场),然后通过在所针对的模型单元上,对上述两场的相应场量点乘,再进行面积积分,就可得到观测量对模型电导率参数的偏导数．模型上的计算结果和计算精度的分析将在本文的续篇(Ⅱ)中介绍。     

Pressure wind adjustment relationships during a multilevel primitive equation prediction process using tropical atmospheric data

Summary In numerical weather forecasting process, with primitive equations, the wind and pressure fields mutually adjust to each other until some form of balance is achieved. The type of balance so achieved by the mass and wind fields during the numerical integration of the primitive equations governing atmospheric motions is not knowna priori. This is particularly so in the case of tropical regions where the pressure wind adjustment laws prevailing in a tropical atmosphere are not well understood.In this study we perform a systematic investigation of the pressure wind adjustment relations during a numerical integration of the primitive equations governing atmospheric motions in a tropical atmosphere. Therefore, a two-day prediction experiment is carried out using the Florida State University Tropical Prediction (FSU) Model (Krishnamurti, 1969;Krishnamurti,et al. 1973;Kanamitsu, 1975). The 200 mb predicted motion (u, v) and height (z) fields are then extracted at 0, 12, 24, 36 and 48 hours of forecast time. Using these motion (u, v) fields three other 200 mb height (z) fields were computed from the inverse nonlinear, linear and quasigeostrophic balance equations. Each of these three diagnostic heights for the 200 mb pressure surface were compared with the respective 200 mb heights obtained from the Florida State University Tropical Preciction Model. The comparison is done by computing the root-mean-square differences between the predicted 200 mb height fields and each of the three 200 mb heights obtained from the inverse non-linear, linear and quasigeostrophic balance equations. The results show that the root-meansquare differences between thez fields from the FSU model and those obtained from the non-linear and linear balance equations lie within the ranges 23 to 44 and 25 to 50 metres respectively. The root-mean-square differences between the predicted heights and the heights computed from the quasigeostrophic balance equation lie in the range 54 to 62 metres. These root-mean-square differences are of significant magnitude since large-scale disturbances in the tropical atmosphere are associated with rather small pressure changes.The variations of these root-mean-square differences as one moves from one forecast time to another exhibit no clear increasing or decreasing trend. In fact the variations appear somewhat random. This rather unsystematic time variation of the root-mean-square differences is a manifestation of the constant changes of the physics in the model as different weather systems evolve in the course of the forecasting process. It seems therefore that the pressure-wind adjustments that take place during a numerical integration of the model equations are of complex nature and cannot simply be approximated by simple diagnostic relations like the ones used in this study.Most of this work was done while the author was at the Florida State University, tallahassee, USA.