航空重力测量中差分器的设计与比较

孙中苗， 孔　超， 温　波

1.西安测绘研究所，陕西 西安，710054；2.地理信息工程国家重点实验室，陕西 西安，710054；3.信息工程大学地理空间信息学院，河南 郑州，450052



航空重力测量中差分器的设计与比较

孙中苗1,2， 孔超3， 温波3

1.西安测绘研究所，陕西 西安，710054；2.地理信息工程国家重点实验室，陕西 西安，710054；3.信息工程大学地理空间信息学院，河南 郑州，450052

航空重力测量中垂直加速度的计算需要差分运算。本文给出了5类常用差分器的系数及传递函数，以两点中心差分、2阶牛顿差分、3阶牛顿差分、2阶Lanczos、移动平均差分和49阶等波纹差分器为例，对差分器的性能做了分析比较。与理想差分器比较，结果表明，6种差分器的幅频响应在0～0.01Hz的低频段均与理想差分器相一致，而在0.01～0.10Hz 频段只有3阶牛顿差分器和等波纹差分器吻合较好；与静态GPS测量数据比较，结果表明，6种差分器获得的垂直加速度经20s低通滤波后标准差几近相同，鉴于两点中心差分器的边界效应最小，可在实践中加以采用。

航空重力测量，垂直加速度，差分器，传递函数

1　引　言

航空重力测量在获取地面重力测量难以到达地区的地球重力场信息方面已得到广泛应用[1-3]。目前实用的航空标量重力测量系统包括平台式和捷联式两类。平台式将航空重力仪固定在阻尼平台上进行测量，如LCR航空重力仪[4]；捷联式是捷联式惯性导航和GPS组合的重力测量系统[5]。无论是平台式系统还是捷联式系统，获取高质量成果的前提是具备良好的飞行作业条件和完善的数据处理算法，其中，提高因飞机运动产生的垂直加速度的确定精度一直是重中之重。

垂直加速度目前一般利用GPS和数字滤波技术相结合进行确定。惯用的做法是，先利用GPS观测数据计算出运动载体的位置或速度，然后差分(二次或一次)获取加速度，最后作低通滤波以消除高频噪声的影响。大量文献研究表明[6-10]，在良好的环境条件下，对于90s～300s的滤波周期，利用上述方法可以1～2mGal的精度导出运动平台的加速度。但由垂直加速度的频谱特性分析表明[10]，垂直加速度不仅占据相当宽的频带，而且频谱特性随测量环境(如不同的运载平台和大气湍流条件)的变化也不尽相同。因此，若欲同时提高航空重力测量的精度和分辨率，在与航空重力相关的频带上，垂直加速度的精密确定仍是主要障碍。

大体而言，垂直加速度的确定精度取决于三个方面：一是GPS定位或测速的精度，主要体现在GPS各种误差源的影响；二是差分方法；三是低通滤波器的性能，它反映了垂直加速度精度所对应的分辨率。本文主要对其中差分方法的影响进行分析，重点比较分析了六种常用差分器的性能，并以此给出了建议的差分方法。

2　垂直加速度计算方法

假设利用GPS观测数据得到的等距离散位置序列为xn(n=1,…,NP)，则位置序列的一次差分即速度序列v(n)可表示为如下卷积形式：

(1)

式中，NP是位置序列的长度；h(k)表示长度为(2M+1)的有限冲击响应(FIR)差分器的系数。

同样，由式(1)的速度序列可得加速度序列a(n)为：

(2)

显然，由式(2)确定的加速度的精度除与位置或速度精度有关外，还取决于所用差分器的性能。

3　常用差分器的系数及传递函数

3.1理想差分器

一阶理想差分器的传递函数为：

Hd(ejω)=jω

(3)

其幅频响应为过原点的一条斜线。

航空重力测量数据中含有许多高频噪声，故希望使用低通差分滤波器，其传递函数为一截断函数：

(4)

式中απ是低通差分器的截止频率。

3.2单纯M次差分

这是最简单的差分算法，记：

Δk=x(n+k)-x(n-k)

(5)

则这种差分器的输入输出关系为：

(6)

其传递函数为：

(7)

若M=1，即为两点中心差分：

(8)

(9)

3.3牛顿-柯斯特差分

这种差分器源自牛顿-柯斯特积分公式，其传递函数是：

(10)

当M=2时，

(11)

(12)

当M=3时，

(13)

(14)

3.4多项式拟合差分器(Lanczos差分器)

(15)

当M=2时，

(16)

(17)

3.5平滑化差分

这是一种数据平滑和差分相结合的算法，目的在于在低频段H(ejω)能更好地近似于jω，而在高频段能有更大的衰减。

(18)

其频率特性为：

(19)

可知H(ejω)是两项的乘积，第一项表现了加权平滑特性，第二项体现了差分特性。

(20)

(21)

当N=1，L=1时，

(22)

(23)

3.6等波纹差分器

在FIR滤波器的等波纹设计法中[11]，利用契比雪夫“交错点组定理”和数值分析中的Remez算法，可惟一地确定滤波器系数组{a(n)}，由此可得到滤波器的传递函数H(ejω)。设计低通滤波器时逼近的是幅频响应在通带和阻带的函数值，差分器的设计过程与低通滤波器完全一样，只是其逼近的是频率响应在各带的斜率。

以下讨论中本文利用等波纹设计法设计了49阶的等波纹差分器。

4　差分器的性能比较与分析

4.1与理想差分器的直观比较

以常用的两点中心差分、2阶牛顿差分、3阶牛顿差分、2阶Lanczos、移动平均差分和49阶等波纹差分器为例，为直观地评价上述各类差分器的性能，图1(a)示出了各类差分器的幅频响应，其与理想差分器幅频响应的差值示于图1(b)。显而易见，在0～0.01Hz的低频段，各个差分器的幅频响应均与理想差分器趋于一致，但在0.01～0.10Hz仅有3次牛顿-柯斯特差分器、等波纹差分器与其相吻合。因此，若目标是差分获取低于0.01Hz频段内的信号，可采用上述差分器中的任何一种；若需恢复0.01Hz以上频谱，应当采用3次以上牛顿-柯斯特差分器或高阶等波纹差分器。

图1　各类差分器的幅频响应及与理想差分器的比较

4.2实测数据比较分析

按照航空重力测量作业要求，飞机起飞前需要进行10～15m的静态测量。静态测量时，可以认为飞机在作速度/加速度为零的运动，因此将所求速度/垂直加速度与真值'0'进行比较，可以估算速度/垂直加速度的确定精度。实验数据源于2002年3月进行的大同航空重力测量的静态GPS观测数据[10]。

为了对直观比较结果作进一步说明，图2(a)、(b)分别示出了利用飞机起飞前静态测量数据计算的高程和速度的幅度谱。表1统计了由高程数据导出的速度的误差，显然，差分器E、F的精度最高，差分器B、C精度次之，而差分器D、G精度最差。这是由于高程的能量主要集中在0.01Hz以内，差分器E、F、B、C代之以在>0.01Hz频段内产生差分误差，其所具有的低通特性部分地将高频噪声进行了阻尼；而差分器D、G尽管在更宽的频带上接近理想差分器，却在高频段的差分运算中放大了高频噪声。

表2统计了由垂直速度导出的垂直加速度的标准差，未滤波时可得出上述类似结论，20s低通滤波后，6种差分器的精度相当。

航空重力测量中，经差分运算获得加速度后，通常需要采用90～300s滤波尺度的低通滤波器进一步消除高频噪声的影响。这一过程实质上是差分器和低通滤波器的卷积运算。从式(8)可以看出，两点中心差分器的边界效应只有1，且其使用简便，建议在实践中加以采用。

图2　高程和速度的幅频响应(静态测量)

表1静态高程数据计算的速度的误差(mm/s)

比较项目╲差分器BCDEFG差值均值-0.06-0.06-0.06-0.06-0.06-0.07标准偏差2.222.823.141.301.363.99

表2静态垂直速度导出的加速度的误差(mGal)

滤波尺度╲差分器BCDEFG0s175.6233.7261.868.677.7319.75s71.8783.1385.7446.7650.6486.8210s20.421.621.717.818.221.720s4.84.94.94.74.74.9

5　结束语

航空重力测量中垂直加速度的精确计算至关重要。垂直加速度通常由高程或速度通过差分运算获得。本文比较分析了六种常用差分器的性能，初步结论如下。

(1)与理想差分器幅频响应进行直观比较，可以看出，两点中心差分、2阶牛顿差分等6种差分器的幅频响应在0～0.01Hz的低频段均与理想差分器趋于一致，但在0.01～0.10Hz仅有3阶牛顿差分器和等波纹差分器吻合较好；如欲获取0.01Hz以上高频重力信号，应当采用3次以上牛顿差分器或高阶等波纹差分器。

(2)从静态GPS测量数据比较分析看，6种差分器获得的垂直加速度经20s低通滤波后，其标准差相差很小。

(3)航空重力测量中垂直加速度的计算实质是差分器和低通滤波器的卷积运算，鉴于两点中心差分器的边界效应最小，可在实践中加以采用。

[1]Brozena J, Childers V A. The NRL Airborne Geophysics Program [A]. In International Association of Geodesy Symposia, Geodesy Beyond 2000-The Challenges of the First Decade. Heidelburg: Springer Verlag,2000.

[2]Forsberg R, Olesen A V. Airborne Gravity and Geoid Surveys in the Arctic and Baltic Seas [A]. In Proceedings of IAG International Symposium on Gravity, Geoid and Geodynamics 2001. Banff, Canada,2001.

[3]Klingele E E, Cocard M, Halliday M E, Kahle H G. The Airborne Gravimetric Survey of Switzerland [M]. Contribution to the Geology of Switzerland Series, Swiss Geophysical Commission, Report No.31,1995.

[4]Valliant H D. The LaCoste and Romberg Air/Sea Gravity Meter: an Overview [M]. CRC Handbook of Geophysical Exploration at Sea, 2nd edition,1991.[5]Bruton A M. Improving the Accuracy and Resolution of SINS/GPS Airborne Gravimetry [D]. Calgary: Department of Geomatics Engineering at University of Calgary,2001.[6]Czompo J. Airborne Scalar Gravimetry Systems in the Spectral Domain [D]. Calgary: Department of Geomatics Engineering at the University of Calgary,2004.[7]Bruton A M. Improving the Accuracy and Resolution of SINS/GPS Airborne Gravimetry [D]. Calgary: Department of Geomatics Engineering at University of Calgary,2001.[8]李晓斌. 航空重力测量中GPS定位、速度及加速度解算技术与应用[D]. 北京：中国地质大学,2010.

[9]欧阳永忠. 海空重力测量数据处理关健技术研究[D]. 武汉：武汉大学,2013.

[10]孙中苗. 航空重力测量理论、方法及应用研究[D]. 郑州：信息工程大学,2004.

[11]胡广书. 数字信号处理：理论,算法与实现(第二版)[M]. 北京:清华大学出版社,2003.

Design and Analysis of the Differentiator in Airborne Gravimetry

Sun Zhongmiao1,2，Kong Chao3, Wen Bo3

1 Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China 2 State Key Laboratory of Geo-information Engineering, Xi’an 710054, China 3.Institute of Geospatial Information, Information Engineering University, Zhengzhou 450052, China

Difference is needed to calculate the vertical acceleration calculation in airborne gravimetry. The coefficients and transfer functions of five types of common differentiators are provided in this paper. Six differentiators are set as examples to compare their performance, including differentiator of two centers, Newton differentiator of second order, Newton differentiator of third order, Lanczos differentiator of second order, moving average differentiator and equal-ripple differentiator of forty-ninth order. Compared with the ideal differentiator, it is indicated that the amplitude frequency responses of all the six differentiators are consistent with those of ideal differentiator in the frequency bands of 0～0.01Hz, while in the frequency bands of 0.01～0.10Hz only the Newton differentiator of third order and the equal-ripple differentiator will coincide well. The static GPS data indicates that the standard deviations of the vertical accelerations achieved by above six differentiators respectively are almost the same after using a lowpass filter with cutoff frequency of 20 seconds. Since the edge effect of the differentiator of two centers is the least, it can be applied to practice.

airborne gravimetry; vertical acceleration; differentiator; transfer function

2015-12-24。

孙中苗(1968—)，男，研究员，主要从事航空重力测量理论和方法研究。

P223

A