Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models

We discuss a new method for approximately decorrelating and non-isotropically filtering the monthly gravity fields provided by the gravity recovery and climate experiment (GRACE) twin-satellite mission. The procedure is more efficient than conventional Gaussian-type isotropic filters in reducing stripes and spurious patterns, while retaining the signal magnitudes. One of the problems that users of GRACE level 2 monthly gravity field solutions fight is the effect of increasing noise in higher frequencies. Simply truncating the spherical harmonic solution at low degrees causes the loss of a significant portion of signal, which is not an option if one is interested in geophysical phenomena on a scale of few hundred to few thousand km. The common approach is to filter the published solutions, that is to convolve them with an isotropic kernel that allows an interpretation as smoothed averaging. The downside of this approach is an amplitude bias and the fact that it neither accounts for the variable data density that increases towards the poles where the orbits converge nor for the anisotropic error correlation structure that the solutions exhibit. Here a relatively simple regularization procedure will be outlined, which allows one to take the latter two effects into account, on the basis of published level 2 products. This leads to a series of approximate decorrelation transformations applied to the monthly solutions, which enable a successive smoothing to reduce the noise in the higher frequencies. This smoothing effect may be used to generate solutions that behave, on average over all possible directions, very close to Gaussian-type filtered ones. The localizing and smoothing properties of our non-isotropic kernels are compared with Gaussian kernels in terms of the kernel variance and the resulting amplitude bias for a standard signal. Examples involving real GRACE level 2 fields as well as geophysical models are used to demonstrate the techniques. With the new method, we find that the characteristic striping pattern in the GRACE solutions are much more reduced than Gaussian-filtered solutions of comparable signal amplitude and root mean square.     

Transformation between surface spherical harmonic expansion of arbitrary high degree and order and double Fourier series on sphere

In order to accelerate the spherical harmonic synthesis and/or analysis of arbitrary function on the unit sphere, we developed a pair of procedures to transform between a truncated spherical harmonic expansion and the corresponding two-dimensional Fourier series. First, we obtained an analytic expression of the sine/cosine series coefficient of the \(4 \pi \) fully normalized associated Legendre function in terms of the rectangle values of the Wigner d function. Then, we elaborated the existing method to transform the coefficients of the surface spherical harmonic expansion to those of the double Fourier series so as to be capable with arbitrary high degree and order. Next, we created a new method to transform inversely a given double Fourier series to the corresponding surface spherical harmonic expansion. The key of the new method is a couple of new recurrence formulas to compute the inverse transformation coefficients: a decreasing-order, fixed-degree, and fixed-wavenumber three-term formula for general terms, and an increasing-degree-and-order and fixed-wavenumber two-term formula for diagonal terms. Meanwhile, the two seed values are analytically prepared. Both of the forward and inverse transformation procedures are confirmed to be sufficiently accurate and applicable to an extremely high degree/order/wavenumber as \(2^{30}\,{\approx }\,10^9\). The developed procedures will be useful not only in the synthesis and analysis of the spherical harmonic expansion of arbitrary high degree and order, but also in the evaluation of the derivatives and integrals of the spherical harmonic expansion.     

GRACE时变重力场各向异性高斯组合滤波方法

受测量误差等因素影响,直接使用GRACE时变重力场模型的地表质量变化反演结果呈现严重条带噪声,必须采用滤波消除。本文对不同滤波方法进行了试验分析,以信噪比最大为准则,确定了不同滤波方法的最优滤波参数,并在此基础上提出了一种各向异性组合滤波方法。该方法根据时变重力场模型球谐系数误差特性,结合各向异性高斯滤波和均方根滤波特点,对精度较高的低次项系数采用较大权重以保留更多有效信号,而对精度较差的高次项系数采用较小权重以压制噪声。不同于传统的两步法组合滤波,该方法仅需进行一步滤波处理。试验结果表明,本文提出的各向异性组合滤波方法计算步骤简单,能够有效消除条带噪声;与单一滤波和传统两步法组合滤波方法相比,提高了反演结果信噪比,保留了更多真实信号。     

Spherical harmonic analysis and synthesis for global multiresolution applications

 On the Earth and in its neighborhood, spherical harmonic analysis and synthesis are standard mathematical procedures for scalar, vector and tensor fields. However, with the advent of multiresolution applications, additional considerations about convolution filtering with decimation and dilation are required. As global applications often imply discrete observations on regular grids, computational challenges arise and conflicting claims about spherical harmonic transforms have recently appeared in the literature. Following an overview of general multiresolution analysis and synthesis, spherical harmonic transforms are discussed for discrete global computations. For the necessary multi-rate filtering operations, spherical convolutions along with decimations and dilations are discussed, with practical examples of applications. Concluding remarks are then included for general applications, with some discussion of the computational complexity involved and the ongoing investigations in research centers. Received: 13 November 2000 / Accepted: 12 June 2001     

Spherical cap harmonic model for mapping and predicting regional TEC

An approach to modeling the regional ionospheric total electron content (TEC) based on spherical cap harmonic analysis is presented. This approach not only provides a better regional TEC mapping accuracy, but also the capability for ionospheric model prediction based on spectrum analysis and least squares collocation. Unlike conventional approaches, which predict the immediate TEC with models using current observations, the spherical cap harmonic approach utilizes models using past observations to predict a model which will provide future TEC values. A significant advantage in comparison with conventional approaches is that the spherical cap harmonic approach can be used to predict the long-term TEC with reasonable accuracy. This study processes a set of GPS data with an observation time span of 1 year from two GPS networks in China. The TEC mapping accuracy of the spherical cap harmonic model is compared with the polynomial model and the global ionosphere model from IGS. The results show that the spherical cap harmonic model has a better TEC mapping accuracy with smoother residual distributions in both temporal and spatial domains. The TEC prediction with the spherical cap harmonic model has been investigated for both short- and long-term intervals. For the short-term interval, the prediction accuracies for the latencies of 1-day, 2-days, and 3-days are 2.5 TECU, 3.5 TECU, and 4.5 TECU, respectively. For the long-term interval, the prediction accuracy is 4.5 TECU for a 2-month latency.     

结合Gram-Schmidt变换的高光谱影像谐波分析融合算法

针对高光谱影像谐波分析融合(HAF)算法在影像融合时不顾及地物光谱曲线整体反射率这一缺陷,提出了结合Gram-Schmidt变换的高光谱影像谐波分析融合(GSHAF)改进算法。GSHAF算法可在完全保留融合前后像元光谱曲线波形形态的基础上,将高光谱影像融合简化为各像元光谱曲线的谐波余相组成的二维影像与高空间分辨率影像之间的融合。它是在原始高光谱影像光谱曲线被谐波分解为谐波余项、振幅和相位后,首先将其谐波余项与高空间分辨率影像进行GS变换融合,这样便可有效地修正融合后像元光谱曲线的反射率特征,随后再利用该融合影像与谐波振幅、相位进行谐波逆变换,完成高光谱影像谐波融合。本文最后利用Hyperion高光谱遥感影像与ALI高空间分辨率影像对GSHAF算法进行可行性分析,再以HJ-1A等卫星数据对其进行普适性验证,试验结果表明,GSHAF算法不仅可以完全地保留光谱曲线波形形态,而且融合后影像的地物光谱曲线反射率更接近真实地物。     

Signal Representation and the Nonlinear Property of Phase Images

This letter proposes a representation of pixel value in a phase image as a complex number, which implies that the phase image is a nonlinear field. Consequently, conventional linear filtering will not properly work. Local statistics filters (LSFs) are chosen as a special case of the signal-processing algorithm to show the significance of signal representation to filtering performance. First, relationships among various LSFs, namely, ordinary LSF, speckle LSF, and phasor LSF, are discussed. We found that they are related to each other, whose general form is a complex-valued LSF (CV-LSF). It is also realized that directional windowing for the CV-LSF is actually not required, which indicates that the complex-valued formulation (representation and algorithm) is the most suitable one for phase image filtering. We demonstrate that the linear filtering to a simulated complex-valued image and an actual interferometric synthetic aperture radar image will start to fail when the window size is subsequently enlarged due to the nonlinear nature of the phase image.     

Global spherical harmonic computation by two-dimensional Fourier methods

A method is presented for performing global spherical harmonic computation by two-dimensional Fourier transformations. The method goes back to old literature (Schuster 1902) and tackles the problem of non-orthogonality of Legendre-functions, when discretized on an equi-angular grid. Both analysis and synthesis relations are presented, which link the spherical harmonic spectrum to a two-dimensional Fourier spectrum. As an alternative, certain functions of co-latitude are introduced, which are orthogonal to discretized Legendre functions. Several independent Fourier approaches for spherical harmonic computation fit into our general scheme.     

Spherical harmonic analysis,aliasing, and filtering

The currently practiced methods of harmonic analysis on the sphere are studied with respect to aliasing and filtering. It is assumed that a function is sampled on a regular grid of latitudes and longitudes. Then, transformations to and from the Cartesian plane yield formulations of the aliasing error in terms of spherical harmonic coefficients. The following results are obtained: 1) The simple quadratures method and related methods are biased even with band-limited functions. 2) A new method that eliminates this bias is superior to Colombo's method of least squares in terms of reducing aliasing. 3) But, a simple modification of the least-squares model makes it identical to the new method as one is the dual of the other. 4) The essential elimination of aliasing can only be effected with spherical cap averages, not with the often used constant angular block averages.     

利用零相位Kaiser窗滤波器改善MEX多普勒数据定轨精度

精密轨道确定在深空探测中至关重要,而定轨数据中的白噪声会影响定轨性能。基于零相位分析,比较了FRR（forward-filter reverse-filter reverse-output）、RRF（reverse-filter reverse-filter forward-output）和Matlab中的filtfilt这3种滤波器的优劣,设计了一种零相位Kaiser窗低通滤波器。利用火星快车号（Mars Express,MEX）的仿真数据和实测数据验证了零相位Kaiser窗低通滤波器的性能,结果发现滤除白噪声后MEX数据的定轨精度有了显著改善。双程测速数据残差均方根（root mean square,RMS）减小为原来的1/3左右,达到了0.031 mm/s;轨道位置和速度与欧空局（European Space Agency,ESA）精密轨道的差异明显变小。该滤波算法作为定轨前的数据预处理可以提高定轨精度,从而为中国火星探测器的轨道数据处理提供一定的参考。     

Shaping filters and the upward continuation problem

An attempt is made to bridge the gap between closed-form harmonic upward continuation (HUC) of analytic covariance functions of the disturbing potential of the anomalous local gravity field and the numerical shaping filter construction when the local gravity vector is modelled in the framework of Kalman filtering. Some fundamental concepts of the local gravity field, interpreted as a stochastic process that is stationary in the plane and harmonic in the upper half space, are reviewed. The shaping-filter modelling technique for the local gravity vector is introduced. To determine the relation between the disturbing potential covariance function and the gravity vector covariance matrix, the role of the so-called admissible pair is established. It is shown that rescaling an admissible pair leads to an analogue rescaling of the shaping filter matrices derived hereof; no cumbersome numerical recalculations are necessary. The class of covariance functions whose corresponding shaping filters possess a closed-form HUC are identified as models whose HUC can be interpreted as a rescaling. Received: 17 December 1997 / Accepted: 7 September 1998     

The Meissl scheme for the geodetic ellipsoid

We present a variant of the Meissl scheme to relate surface spherical harmonic coefficients of the disturbing potential of the Earth’s gravity field on the surface of the geodetic ellipsoid to surface spherical harmonic coefficients of its first- and second-order normal derivatives on the same or any other ellipsoid. It extends the original (spherical) Meissl scheme, which only holds for harmonic coefficients computed from geodetic data on a sphere. In our scheme, a vector of solid spherical harmonic coefficients of one quantity is transformed into spherical harmonic coefficients of another quantity by pre-multiplication with a transformation matrix. This matrix is diagonal for transformations between spheres, but block-diagonal for transformations involving the ellipsoid. The computation of the transformation matrix involves an inversion if the original coefficients are defined on the ellipsoid. This inversion can be performed accurately and efficiently (i.e., without regularisation) for transformation among different gravity field quantities on the same ellipsoid, due to diagonal dominance of the matrices. However, transformations from the ellipsoid to another surface can only be performed accurately and efficiently for coefficients up to degree and order 520 due to numerical instabilities in the inversion.     

全球电离层TEC的轮胎调和分析与预报

将轮胎调和分析引入电离层TEC的模型化过程中,建立了基于轮胎调和分析的电离层TEC球谐系数模型,并对该模型进行详细的验证和分析。结果表明,本文模型计算精度高,系数截断为15阶时,恢复误差全年统计不超过4%,且除南、北极区外球谐模型具有很好的适用性。然后对该模型系数的时间序列特性进行了函数估计：引入逐级余差建模方法,使用趋势函数、功率谱分析、ARMA模型对球谐系数的时间序列进行分析,找出了模型系数时间序列变化的规律,构建了预报模型,实现了基于模型系数的预报,并对预报系数的精度变化问题和系数本身短期预报的数据积累时间进行分析,最终通过TEC的预报,验证了模型的精度。     

Basic equations for constructing geopotential models from the gravitational potential derivatives of the first and second orders in the terrestrial reference frame

This research represents a continuation of the investigation carried out in the paper of Petrovskaya and Vershkov (J Geod 84(3):165–178, 2010) where conventional spherical harmonic series are constructed for arbitrary order derivatives of the Earth gravitational potential in the terrestrial reference frame. The problem of converting the potential derivatives of the first and second orders into geopotential models is studied. Two kinds of basic equations for solving this problem are derived. The equations of the first kind represent new non-singular non-orthogonal series for the geopotential derivatives, which are constructed by means of transforming the intermediate expressions for these derivatives from the above-mentioned paper. In contrast to the spherical harmonic expansions, these alternative series directly depend on the geopotential coefficients ${\bar{{C}}_{n,m}}$ and ${\bar{{S}}_{n,m}}$ . Each term of the series for the first-order derivatives is represented by a sum of these coefficients, which are multiplied by linear combinations of at most two spherical harmonics. For the second-order derivatives, the geopotential coefficients are multiplied by linear combinations of at most three spherical harmonics. As compared to existing non-singular expressions for the geopotential derivatives, the new expressions have a more simple structure. They depend only on the conventional spherical harmonics and do not depend on the first- and second-order derivatives of the associated Legendre functions. The basic equations of the second kind are inferred from the linear equations, constructed in the cited paper, which express the coefficients of the spherical harmonic series for the first- and second-order derivatives in terms of the geopotential coefficients. These equations are converted into recurrent relations from which the coefficients ${\bar{{C}}_{n,m}}$ and ${\bar{{S}}_{n,m}}$ are determined on the basis of the spherical harmonic coefficients of each derivative. The latter coefficients can be estimated from the values of the geopotential derivatives by the quadrature formulas or the least-squares approach. The new expressions of two kinds can be applied for spherical harmonic synthesis and analysis. In particular, they might be incorporated in geopotential modeling on the basis of the orbit data from the CHAMP, GRACE and GOCE missions, and the gradiometry data from the GOCE mission.     

基于不同潮汐类型的不同时间尺度分潮结果精度分析

选取4种不同潮汐类型的验潮站实测潮位资料,对验潮零点在一定范围内变动对调和分析结果的影响进行分析,结果表明,一年的潮位数据中若后续6个月的数据发生10 cm、20 cm、30 cm零点漂移,各分潮振幅最大变化量为0.47 cm,迟角为0.16°。采用1~12个月不同中期尺度的实测数据,对调和结果的精度进行分析,结果表明,当潮位观测数据时间尺度小于6个月时,其4个主要分潮O_1、K_1、M_2、S_2振幅综合中误差在2 cm以上,迟角中误差在2°以上,具有显著的不稳定性;相同时间尺度的不同潮汐类型的潮位站潮位资料调和分析得到的主要分潮的振幅精度差异较大,迟角分布相似。当潮位观测数据时间尺度达到或多于6个月时,澳门站、汕尾站的分潮综合中误差分别为1.16 cm、0.61 cm,厦门站、北海站分别为2.90 cm和2.51 cm;各分潮迟角中误差均在2°内。当时间尺度超过9个月后,4个验潮站分潮振幅综合中误差均在2 cm以内,各分潮迟角中误差均在1°左右。     

各向异性组合滤波法反演陆地水储量变化

地球时变重力场模型反演陆地水储量变化已为全球气候变化研究作出巨大贡献,考虑到时变重力场模型球谐系数中存在相关性,其高阶次项具有较大的误差,需采用最优的滤波方法进行空间平滑。本文提出一种新的各向异性组合滤波方法,其基本思想是将改进的高斯滤波法与均方根(root mean square,RMS)滤波法组合,即对球谐系数的低阶次采用改进的高斯滤波法,而高阶次采用RMS滤波法。首先分析了最新的GRACE RL05系列时变重力场模型系数误差特性,基于全球水储量变化反演结果,分析比较了高斯滤波、改进的高斯滤波、RMS滤波和DDK滤波与本文提出的组合滤波法的有效性及精度,并利用模型结果进行了验证,计算结果表明,组合滤波法的中误差最小。研究结果表明,本文提出的组合法相比于先前的滤波方法,可有效地过滤高阶次的噪声,消除南北条带误差,同时减少信号泄漏,提高信噪比,保留更多有效的地球物理信号,进而提高反演精度。     

一种基于小波相位分析的InSAR干涉图滤波算法

为了提高InSAR干涉图的滤波质量,本文在考虑干涉信号在同一尺度内具有相关性的基础上,提出并实现了一种基于小波相位分析的干涉图滤波算法.为检验算法的有效性,选取美国Phoenix地区高信噪比和伊朗Bam地区低信噪比两个真实C波段干涉图进行了滤波实验,通过与中值滤波算法、小波-中值滤波算法比较,证实了小波相位滤波算法是有效的且具有一定的优势.研究表明,不论是高信噪比干涉图还是低信噪比干涉图,小波相位滤波算法都能有效抑制其噪声.     

椭球谐和球谐系数之间一个简单的转换关系

本文推导的椭球谐系数和球谐系数相互之间转换关系的核心思想是在ε~2量级下利用Legendre函数的正交性,从球谐系数求解的积分表示出发,将积分中的椭球坐标变量与球坐标变量相互转换,从而得出椭球谐系数与球谐系数之间的转换关系。本文导出的转换关系有以下优点:①对于第二类Legendre函数的计算采用Laurent级数表示,使计算第二类Legendre函数更为简单;②保留了ε~2量级下,导出的转换关系相比文献[2]的形式更简单,满足物理大地测量边值问题线性化的要求;③顾及了余纬和归化余纬的区别。     

应用Slepian局部谱方法解算中国大陆重力场球谐模型

根据经典的球谐函数方法,为满足正交化要求,观测数据需要覆盖整个球面,而对于地表局部测量数据,则无法应用球谐方法解算重力场模型。针对此问题,采用Slepian局部谱分析方法解算中国大陆范围内的实测重力场变化数据,并以GOCE卫星球谐函数解作为已知模型,评估由于实际陆地重力测点的非均匀分布对球谐函数解的误差影响。通过计算多个阶次中国大陆局部范围的Slepian基函数分布;采用GOCE卫星获得重力场模型的前72阶球谐系数作为已知结果,评价实际测点非均匀分布的解算有效性,并针对中国大陆地区采用Slepian基函数进行解算,通过模型对比选择最优截段项数;针对2005—2008年中国大陆地区流动重力测量获得的重力场变化信号进行解算,获得了72阶重力场变化模型。     

Rectangular rotation of spherical harmonic expansion of arbitrary high degree and order

In order to move the polar singularity of arbitrary spherical harmonic expansion to a point on the equator, we rotate the expansion around the y-axis by \(90^{\circ }\) such that the x-axis becomes a new pole. The expansion coefficients are transformed by multiplying a special value of Wigner D-matrix and a normalization factor. The transformation matrix is unchanged whether the coefficients are \(4 \pi \) fully normalized or Schmidt quasi-normalized. The matrix is recursively computed by the so-called X-number formulation (Fukushima in J Geodesy 86: 271–285, 2012a). As an example, we obtained \(2190\times 2190\) coefficients of the rectangular rotated spherical harmonic expansion of EGM2008. A proper combination of the original and the rotated expansions will be useful in (i) integrating the polar orbits of artificial satellites precisely and (ii) synthesizing/analyzing the gravitational/geomagnetic potentials and their derivatives accurately in the high latitude regions including the arctic and antarctic area.