2. 山东大学空间科学研究院, 山东威海 264209;
3. 中国测绘科学研究院, 北京 100830
2. Institute of Space Sciences, Shandong University, Shandong Weihai 264209, China;
3. Chinese Academy of Surveying and Mapping, Beijing 100830, China
我国重力空白区域在我国国土面积中占有很大比重,并且多集中在难以到达的沙漠、水域、高山等区域(陈俊勇等,2001;周坚鑫等,2003;张昌达,2005).传统的静态重力测量手段难以在这些困难区域展开,卫星重力测量由于分辨率较低,也难以满足需求.而航空重力测量以其高精度、高分辨率、高机动性等特点,成为该类地区获取重力信息的有效手段.航空重力测量的研究始于20世纪50—60年代,在80年代末90年代初逐渐成熟.研究表明航空重力测量的精度为2~10 mGal(分辨率6~10 km)(Brozena et al.,1989;Klingelé et al.,1995;Brozena and Childers,2000;Forsberg et al.,2001; Verdun et al.,2002;Rene et al.,2007;孙中苗等,2013a).商用航空重力仪(例如LaCoste & Romberg、BGM、KSS等产品)具有较高的精度和较好的稳定性,但我国对该类设备的引进存在诸多限制.我国在航空重力测量的研究起步较晚,西安测绘研究所和航天科技集团16所率先研制的航空重力测量系统CHAGS,测量精度达到了2~5 mGal(分辨率10 km)(夏哲仁,2004;孙中苗,2004;孙中苗等,2004a,2004b;孙中苗等,2013a).
基于GNSS/INS的航空重力测量是目前研究最为广泛的航空重力测量技术,根据测量原理的不同,可以分为直接求差法和状态模型法(Jekeli,2001).直接求差法由加拿大Calgary大学的Bruton在2000年提出(Bruton,2000),随后Serpas、Senobar等对该方法进行了深入的研究(Serpas and Jekeli,2005;Senobari,2010).在国内,孙中苗等利用直接求差法对CHAGS系统开展了多项研究(孙中苗,2004;孙中苗等,2013b).尽管直接求差法在航空重力测量领域中有较为深入的研究,但该方法对惯性系统的校准精度和系统姿态的测量精度都有严格的要求.与直接求差法不同,状态模型法把重力扰动当作Kalman滤波系统的状态量进行最优估计,在一定程度上降低了对惯性系统的要求.Tomé、Deurloo等的研究成果表明,利用中低精度的惯性系统,根据状态模型法同样可以得到较高精度的重力扰动结果(Tomé,2002;Deurloo et al.,2011;Bastos et al.,2013).目前国内的研究多集中在直接求差法,对状态模型法的研究较少.
自适应Kalman滤波以其调节状态信息的权阵来控制状态误差对参数估计的影响的特性,在定位与导航有着广泛的研究与应用(Yang et al.,2001,2006a,2006b;Hajiyev and Soken,2013).研究表明引入自适应因子,并根据一些理论或者经验的准则对自适应因子进行合理的估计以后,组合导航系统的性能明显优于经典Kalman滤波(高广为等,2006;Yang et al.,2006a).
本文首先展开对状态模型法航空重力理论研究,给出该方法详细的实现过程,并在状态模型法的基础上引入最优自适应因子算法,调节动力学模型信息和观测信息对滤波结果的影响,改善重力扰动的计算结果,并以新疆地区的实测数据验证该方法的有效性.
2 基本原理根据牛顿第二定律,在惯性坐标系中质点的动力学方程为
其中 为质点的惯性加速度,gi为重力加速度,fi为比力加速度.
在SINS/GNSS组合导航过程中,在导航坐标系n下,速度测量值对时间求偏导可以得到(Titterton and Weston,2004):
其中ven为载体的速度,fb为加速度计的比力测量值,Cnb为b系到n系的方向余弦矩阵,ωien为地球自转角速度在n系下的投影,ωenn为n系相对于地球坐标系e系的角速度在n系下的投影,式中的第二项为哥式加速度和向心加速度在n下的投影,gn为实际重力矢量.gn通常采用正常重力γn和重力扰动δgn的矢量和表示,公式为
在物理大地测量的定义中,重力异常主要包括重力扰动(纯重力异常)和混合重力异常.其中重力扰动是指地球某一点上的重力观测值和正常重力之差.重力扰动δgn为一个三维矢量,由两个水平分量ζg,ηg和一个垂直向下分量δg组成(Heiskanen and Moritz,1967),δgn=[-ζg -ηg δg]T.本文使用的惯性单元为中等精度水平的光纤陀螺惯导,受惯性单元性能的影响,δgn水平分量很难达到期望的精度,因此这里主要讨论垂直分量上的重力扰动δg.目前估计重力扰动δg的方法有两种:直接求差法和基于Kalman滤波的状态估计法,下面对两种方法的基本原理进行介绍.
2.1 直接求差法航空重力测量直接求差法主要是通过公式将重力扰动矢量表示为
重力扰动矢量垂直向下的分量可以详细的表示为
其中(·)D为垂直向下的分量,右边第三项为厄特弗斯改正,L为纬度,h为椭球高,vE为东向速度,为垂直向下加速度,RM为子午圈半径,RN为卯酉圈半径.
利用公式(2)计算的空中重力扰动含有大量的噪声,而重力扰动信号能量通常集中在低频段很窄的部分.因此需要设计合适的低通滤波器对空中重力扰动进行噪声消除.在航空重力测量领域中,比较典型的低通滤波器有RC低通滤波器、Butterworth低通滤波器、FIR低通滤波器.FIR低通滤波器的差分方程形式可以表示为(孙中苗,2004;孙中苗等,2004b)
其中x(n)为输入信号,y(n)为输出信号,h(k)为滤波器系数(单位冲激响应),N为滤波器长度.h(k)的理想频率响应可以表示为
确定单位冲激响应函数主要有切比学夫等纹逼近法和窗函数法,主要做法请参考文献(孙中苗,2004)和(孙中苗等,2004b),本文不再赘述.
2.2 引入最优自适应因子的状态模型法航空重力测量在组合导航领域,位置、速度、姿态、惯性器件的系统误差通常被当作Kalman滤波的系统状态量,通过预计方程和更新方程估计出来.基于状态法的航空重力测量就是通过对重力扰动进行建模,将其作为系统的状态量,利用Kalman滤波进行最优化估计.
以21+1维系统状态的Kalman滤波为例,系统的状态量可以表示为
其中ren、ven分别为位置向量和速度向量,Cbn为从载体坐标系b旋转到导航坐标系n的方向余弦矩阵,ba和bω分别为惯性单元加速度计和陀螺三轴方向的随机常值误差,sa和sω分别为惯性单元加速度计和陀螺三轴方向的尺度因子.δg为重力扰动在垂直向下方向上的分量,假设在短时间内为一常量偏差,也可以采用三阶Gauss-Markov模型进行建模(Jekeli,2001),研究表明两种模型都能很好的对重力扰动进行估计(Deurloo,2011).对应的误差状态向量可以表示为
Kalman滤波的动力学预计方程为: =f(x)+Gw,其中w为过程噪声,主要由加速度计和陀螺观测量的随机噪声组成.利用一阶泰勒公式进一步展开,可以得到误差状态量的预计方程的近似表达方式: ,其中F是系统矩阵,其非线性的表达方式为 ,G和G′为噪声影响矩阵,F和G都是时变矩阵,在本文的应用中表示为
其中F11—F33由位置、速度、姿态的导航方程求偏导得出,读者可以参考相关文献或自行推导;Cbn为从惯性系到导航系的方向余弦矩阵.
一个典型的Kalman滤波的预测方程可以表示为
其中Φk-1近似的公式为
噪声矩阵Qk-1为:Qk-1≈Gk-1QΔtGk-1T.其中功率谱密度矩阵Q可以由过程噪声表示为(Gelb,1974): Q(t)δ(t-τ)=E[w(t)w(t)T],在本文的22维的Kalman滤波中,Q被定义为
其中Qa,Qω分别是加速度计和陀螺随机游走的平方,qδg=σδg2,这里σδg取 (Bruton,2000).
引入最优自适应因子的更新方程可以表示为
其中α为最优自适应因子,zk是第k时刻的更新信息,通常是指GNSS的位置和速度的解算结果,Hk为观测矩阵,Kalman增益矩阵Kk为
其中Rk为更新观测值的协方差矩阵.
tk 时刻预测状态向量的协方差可以表示为
预测向量的不符值为: ,利用Sage开窗估计法,可以求得状态预测向量协方差矩阵的新的估值为(杨元喜,2006)
采用多个历元的平均值并不能有效的反应当前历元的模型误差的大小,因此这里简单的取:
自适应因子的表达式为(Yang and He,2001;高为广等,2006;Yang et al.,2006a,2006b;杨元喜,2006)
其中 ,c在本文中取1.
3 算例及结果分析选取西部某地两天航空摄影测量任务的POS(测姿定位系统,Position and Orientation System)数据进行分析,飞行轨迹分别为图 1和图 2.搭载的POS系统为IGI AeroControl & AeroOffice,该系统包括一款中等精度的光纤陀螺惯导IMU-d和一款NovAtel OEM4 GNSS双频接收机,表 1为该系统利用本文的算法计算的位置和姿态结果与航空摄影测量的空中三角测量结果的较差统计表,从表可以看出利用该系统能提供比较可靠的位置和姿态信息.飞行载体的为运八型飞机.在测区内架设有GPS基站,利用GPS差分技术获取高精度的位置、姿态、速度、加速度信息.测试过程中IMU设备的频率为256 Hz,GNSS差分结果的频率为1 Hz.分别利用直接求差法和自适应状态模型法进行航空重力解算.解算结果与全球重力场模型EGM2008在该区域重力扰动资料进行比较,验证结果的精度.EGM2008全球重力场模型是由NGA(National Geospatial-intelligence Agency)发放全球超高阶地球重力场模型,分辨率最高可达2.5′×2.5′(约为5 km格网),研究表明该模型在我国地势平缓的区域的精度能达到10 mGal左右(章传银等,2009;杨金玉等,2012).由于目前掌握该测区的地面重力资料非常有限,本文暂时用EGM2008全球重力模型计算出的重力扰动值作为参考,对初步解算结果进行对比分析.两次测量都在飞机起飞前进行了系统初始对准,系统的初始对准包括POS系统的初始姿态对准和重力扰动的对准,对准时间约为15 min.由于停机坪附近没有重力基准点,系统采用EGM2008模型推算的重力扰动作为初始重力扰动.EGM2008提供的全球重力扰动结果为参考椭球表面的重力扰动,本文利用NGA网站提供的Fortran脚本“EGMhsynth_WGS84.f”和2.5′×2.5′的格网数据进行计算了测试所在区域椭球表面的重力扰动,然后再通过向上延拓公式延拓到航线所在的高度上,将延拓后的结果作为本文计算结果的参考值.测试中采用的向上延拓算法为直接代表法,计算公式为(石磐和王兴涛,1997;孙中苗,2004)
在k时刻的重力输入值的迭代公式为:gkair=gkh+δgk-1air,然后利用公式(8)—(15)计算εδgk,k时刻的重力扰动就可以表示为:δgkair=δgk-1air+εδgk.其中gh为在航线高度上的正常重力,公式为(NIMA,2000):
其中
这里的参数Φ为纬度,单位为弧度,h为椭球高,单位为m;a,b,e,f是重力模型的椭球参数;参数γp和γe分别为重力模型的极参数和赤道参数;ωe为地球的自转角速度,GM为地球的引力常数.这些地球参数的取值请参照文献(NIMA,2000).
3.1 算例1算例1为对新疆某沙漠地带一平均航高6000 m、航迹340 km、航速155 m·s-1的航线进行航空重力解算,飞行轨迹如图 3.飞机在飞行过程中受到大气扰动较小,飞行高度稳定在6000 m左右,起伏不超过20 m,航线经过区域为沙漠地带,地面起伏不大.由EGM2008计算该区域的重力扰动分布如图 4.
采用直接求差法、状态模型法计算的重力扰动结果与EGM2008模型计算的该航线上的重力扰动如图 5,从图中我们可以看出,利用状态模型法计算出的重力扰动和EGM2008模型的趋势是一致的,自适应状态模型法又较状态模型法的趋势更为平滑.直接求差法的结果与整体趋势偏离程度较大,可能是受POS性能的影响(与商用重力仪有一定差距),计算出的姿态角的精度不是足够高,仍对最后的重力扰动结果带来了影响.表 2为解算结果与EGM2008模型的较差统计表,从表中我们可以看出,利用状态模型法求出的重力扰动与EGM2008模型的差值中误差为8.27 mGal,利用自适应状态模型法的差值中误差为5.43 mGal,接近GB/T17944-2000《加密重力测量规范》困难地区10 mGal的限差要求.
算例2为对新疆某区平均航高5000 m、航迹140 km、航速分别为145、135、137 m·s-1的三条航线进行航空重力解算,飞行轨迹如图 6.飞机在飞行过程中受到的大气扰动较小,三条航线的行高起伏不大.航线经过的区域北部为绿洲,南部为平坦的沙漠,自北向南地面高程逐渐增加,整体无明显起伏.
三条航线的解算结果和EGM2008重力扰动进行作差比较如表 3,从表中可以看出,三条航线差值的中误差都在10 mGal左右,与算例1的结论是一致的.
(1)利用状态模型法可以有效地获取高精度、高分辨率的空中重力扰动成果,引入自适应因子以后,结果的精度也得到了一定程度的改善;与EGM2008全球重力场模型对比的差值中误差在10 mGal左右,接近国家在困难地区重力测量的精度要求.
(2)由于采用的惯性系统精度不是很高,未能完全满足直接求差法的要求.从这一点也可以看出状态模型法对设备的性能的要求要低一些.
(3)受缺少地面重力资料的限制,处理结果仅与开放性的重力场模型进行对比验证,在以后的研究中还需要在含有地面重力资料或高程异常的区域进行航空重力测量,来进一步验证本文提出的方法的有效性.
(4)利用航空摄影测量任务的POS数据进行航空重力测量分析,提高了数据的利用率,扩展了航空重力测量的方法,在一定程度上节约了航空重力测量的成本.
(5)基于状态模型法的航空重力测量方法同时还是一种矢量重力测量方法,目前国内还没有成熟的航空矢量重力测量系统,在航空矢量重力数据处理方面也尚处于起步阶段(宁津生,2014),充分吸收国内外现有研究成果,对我国的重力测量事业有非常重要的意义.
[1] | Bastos L, Yan W, Magalhāes A, et al. 2013. Assessment the performance of low cost IMUs for strapdown airborne gravimetry using UAVs.// 4th International Colloquium Scientific and Fundamental Aspects of the Galileo Programme. Prague, Czech Republic. |
[2] | Brozena J M, Childers V A. 2000. The NRL Airborne Geophysics Program. Berlin: Springer. |
[3] | Brozena J M, Mader G L, Peters M F. 1989. Interferometric global positioning system: Three-dimensional positioning source for airborne gravimetry. Journal of Geophysical Research: Solid Earth (1978—2012), 94(B9): 12153-12162. |
[4] | Bruton A M. 2000. Improving the accuracy and resolution of SINS/DGPS airborne gravimetry[Ph. D. Thesis]. Canada: University of Calgary. |
[5] | Chen J Y, Wen H J, Cheng P F. 2001. Some issues of the development of Geodesy in China. Advance in Earth Sciences (in Chinese), 16(5): 681-688. |
[6] | Deurloo R A. 2011. Development of a Kalman Filter integrating system and measurement models for a Low-cost strapdown airborne gravimetry system[Ph. D. Thesis]. Porto: University of Porto. |
[7] | Deurloo R A, Bastos M L, Geng Y, et al. 2011. Evaluation of low- and medium-cost IMUs for Airborne Gravimetry with UAVs.// American Geophysical Union, Fall Meeting 2011. San Francisco, California, USA. |
[8] | Forsberg R, Olesen A, Keller K, et al. 2001. Airborne gravity and geoid surveys in the arctic and Baltic Seas.// Proceedings of International Symposium on Kinematic Systems in Geodesy, Geomatics and Navigation (KIS-2001). Banff, 586-593. |
[9] | Gao W G, Yang Y X, Cui X Q, et al. 2006. Application of adaptive Kalman Filtering algorithm in IMU/GPS integrated navigation system. Geomatics and Information Science of Wuhan University (in Chinese), 31(5): 466-469. |
[10] | Gelb A. 1974. Applied Optimal Estimation. Massachusetts: MIT Press. |
[11] | Hajiyev C, Soken H E. 2013. Robust adaptive Kalman Filter for estimation of UAV dynamics in the presence of sensor/actuator faults. Aerospace Science and Technology, 28(1): 376-383. |
[12] | Heiskanen W A, Moritz H. 1967. Physical Geodesy. San Francisco: W. H. Freeman & Co. |
[13] | Jekeli C. 2001. Inertial Navigation Systems with Geodetic Applications. Berlin and New York: Walter de Gruyter & Co. |
[14] | Klingelé E, Halliday M, Cocard M. 1995. Airborne gravimetric survey of Switzerland.// 4th International Congress of the Brazilian Geophysical Society. |
[15] | NIMA. 2000. Department of defense world geodetic system 1984. 3rd ed. National Imagery and Mapping Agency. |
[16] | Ning J S. 2014. Research on the data processing methods of airborne vector gravimetry using SINS/GNSS. Engineering Sciences (in Chinese), 16(3): 4-13. |
[17] | Rene F, Olesen A, Munkhtsetseg D, et al. 2007. Downward continuation and geoid determination in mongolia from airborne and surface gravimetry and SRTM topography.// International Forum on Strategic Technology, 2007. IFOST 2007. Ulaanbaatar: IEEE, 470-475. |
[18] | Senobari M S. 2010. New results in airborne vector gravimetry using strapdown INS/DGPS. Journal of Geodesy, 84(5): 277-291. |
[19] | Serpas J G, Jekeli C. 2005. Local geoid determination from airborne vector gravimetry. Journal of Geodesy, 78(10): 577-587. |
[20] | Shi P, Wang X T. 1997. Determination of the terrain surface gravity field using airborne gravimetry and DEM. Acta Geodaetica Et Cartographica Sinica (in Chinese), 26(2): 117-121. |
[21] | Sun Z M. 2004. Theory, method and applications of airborne gravimetry (in Chinese). Zhengzhou: Information Engineering University. |
[22] | Sun Z M, Xia Z R, Shi P. 2004a. Progress and advance in airborne gravimetry. Progress in Geophysics (in Chinese), 19(3): 492-496. |
[23] | Sun Z M, Xia Z R, Shi P, et al. 2004b. Filtering and processing for the airborne gravimetry data. Progress in Geophysics (in Chinese), 19(1): 119-124. |
[24] | Sun Z M, Zhai Z H, Li Y C. 2013b. Status and development of airborne gravimeter. Progress in Geophysics [Ph. D. Thesis](in Chinese), 28(1): 1-8, doi: 10.6038/pg20130101. |
[25] | Sun Z M, Zhai Z H, Xiao Y, et al. 2013a. Systematic error compensation for airborne gravimetry. Chinese Journal of Geophysics (in Chinese), 56(1): 47-52, doi: 10.6038/cjg20130105. |
[26] | Titterton D H, Weston J L. 2004. Strapdown Inertial Navigation Technology. Stevenage: The Institution of Electrical Engineers. |
[27] | Tomé P. 2002. Integration of inertial and satellite navigation systems for aircraft attitude determination . Porto: University of Oporto. |
[28] | Verdun J, Bayer R, Klingelé E E, et al. 2002. Airborne gravity measurements over mountainous areas by using a LaCoste & Romberg air-sea gravity meter. Geophysics, 67(3): 807-816. |
[29] | Xia Z R, Shi P, Sun Z M, et al. 2004. Chinese airborne gravimetry system CHAGS. Acta Geodaetica et Cartographica Sinica (in Chinese), 33(3): 216-220. |
[30] | Yang J Y, Zhang X H, Zhang F F. 2012. On the accuracy of EGM2008 earth gravitational model in Chinese mainland. Progress in Geophysics (in Chinese), 27(4): 1298-1306, doi: 10.6038/j.issn.1004-2903.2012.04.003. |
[31] | Yang Y, He H, Xu G. 2001. Adaptively robust filtering for kinematic geodetic positioning. Journal of Geodesy, 75(2-3): 109-116. |
[32] | Yang Y X. 2006. Adaptive Navigation and Kinematic Positioning (in Chinese). Beijing: Surveying and Mapping Press. |
[33] | Yang Y X, Gao W G. 2006a. An Optimal Adaptive Kalman Filter. Journal of Geodesy, 80(4): 177-183. |
[34] | Yang Y X, Gao W G. 2006b. A new learning statistic for adaptive filter based on predicted residuals. Progress in Natural Science, 16(8): 833-837. |
[35] | Zhang C D. 2005. On the subject of airborne gravimetry and airborne gravity gradiometry. Chinese Journal of Engineering Geophysics (in Chinese), 2(4): 282-291. |
[36] | Zhang C Y, Guo C X, Chen J Y, et al. 2009. EGM 2008 and its application analysis in Chinese Mainland. Acta Geodaetica et Cartographica Sinica (in Chinese), 38(4): 283-289. |
[37] | Zhou J X, Liu H J, Wang S T. 2003. The prediction of China airborne gravimetry in geophysics.// China Geophysical Society 19th Annual Conference (in Chinese). Beijing: Chinese Geophysical Society. |
[38] | 陈俊勇, 文汉江, 程鹏飞. 2001. 中国大地测量科学发展的若干问 题. 地球科学进展, 16(5): 681-688. |
[39] | 高为广, 杨元喜, 崔先强等. 2006. IMU/GPS组合导航系统自适应Kalman滤波算法. 武汉大学学报(信息科学版), 31(5): 466-469. |
[40] | 宁津生. 2014. 基于SINS/GNSS的航空矢量重力测量数据处理方法研究. 中国工程科学, 16(3): 4-13. |
[41] | 石磐, 王兴涛. 1997. 利用航空重力测量和DEM确定地面重力场. 测绘学报, 26(2): 117-121. |
[42] | 孙中苗. 2004. 航空重力测量理论、方法及应用研究[博士论文]. 郑州: 解放军信息工程大学. |
[43] | 孙中苗, 夏哲仁, 石磐. 2004a. 航空重力测量研究进展. 地球物理学进展, 19(3): 492-496. |
[44] | 孙中苗, 夏哲仁, 石磐等. 2004b. 航空重力测量数据的滤波与处理. 地球物理学进展, 19(1): 119-124. |
[45] | 孙中苗, 翟振和, 李迎春. 2013b. 航空重力仪发展现状和趋势. 地球物理学进展, 28(1): 1-8, doi: 10.6038/pg20130101. |
[46] | 孙中苗, 翟振和, 肖云等. 2013a. 航空重力测量的系统误差补偿. 地球物理学报, 56(1): 47-52, doi: 10.6038/cjg20130105. |
[47] | 夏哲仁, 石磐, 孙中苗等. 2004. 航空重力测量系统CHAGS. 测绘学报, 33(3): 216-220. |
[48] | 杨金玉, 张训华, 张菲菲等. 2012. EGM2008地球重力模型数据在中国大陆地区的精度分析. 地球物理学进展, 27(4): 1298-1306, doi: 10.6038/j.issn.1004-2903.2012.04.003. |
[49] | 杨元喜. 2006. 自适应动态导航定位. 北京: 测绘出版社. |
[50] | 张昌达. 2005. 航空重力测量和航空重力梯度测量问题. 工程地球物理学报, 2(4): 282-291. |
[51] | 章传银, 郭春喜, 陈俊勇等. 2009. EGM 2008地球重力场模型在中国大陆适用性分析. 测绘学报, 38(4): 283-289. |
[52] | 周坚鑫, 刘浩军, 王守坦. 2003. 航空重力测量在我国地球物理勘探中的应用展望.// 中国地球物理学会第十九届年会论文集. 北京: 中国地球物理学会. |