航空重力梯度融合多源数据构建区域的扰动全张量v分析.docx

航空重力梯度融合多源数据构建区域的扰动全张量v分析相对于传统的重力场数据,重力梯度数据含有更多的高频信息量,在油气矿产资源勘探、大地测量学、地球物理学以及国防建设中具有较强的优越性和应用价值(Bell,1998;曾华霖,1999;Hatch,2004;张昌达,2005;边少锋和纪兵,2006;舒晴等,2009;熊盛青,2009;叶周润,2010;梁星辉,2012;袁园,2015;刘金钊等,2015;叶周润等,2015).1886年,匈牙利科学家LorándEötvös设计了最早的重力梯度测量系统——扭称,随后几十年,北美、加拿大、澳大利亚、德国、意大利和英国等多个国家对重力梯度测量仪器及其关键技术进行了研制与实验,并于20世纪70到80年代研制成功了基于旋转加速度计技术的航空重力梯度测量系统以及于2007年研制成功的具有更高测量精度的超导重力梯度测量系统,同时基于冷原子技术的重力梯度测量系统也在不断测试发展中(表1)(Murphy,2004;Lumleyetal.,2004;罗嗣成,2007;Anstieetal.,2010;Dransfieldetal,2010;Kieranetal,2000;梁星辉,2012);2009年3月欧洲航天局成功发射了基于重力梯度测量技术的GOCE卫星,使得利用卫星重力梯度数据优化全球重力场模型成为大地测量研究领域的前沿和热点问题之一(Kleesetal.,2000;罗志才等,2009;郑伟等,2010;刘晓刚,2011;Pailetal.,2011;游为等,2012).因此及时开展重力梯度方面的深入研究能够更好地促进我国重力及重力梯度测量关键技术的向前发展,而高分辨率区域重力梯度场建模研究不仅能够为未来我国独立自主设计的重力梯度测量系统的外部检核提供条件,同时对于相关学科(如恢复中高频段重力场信号、区域重力场模型优化及大地水准面精化等;浅层地壳结构精细反演;水下潜艇重力梯度匹配导航、导弹精密轨道控制等)的探索研究及国民经济发展等也具有积极的科学意义(刘金钊,2017).目前,解算区域重力梯度场主要是基于地形高程数据(边少锋和张赤军,1999;张赤军,1999;JekeliandZhu,2006;RózsaandTóth,2007;MakhloofandIlk,2008;武凛等,2009;Eshaghetal,2009;刘繁明等,2013)、重力数据(MickusandHinojosa,2001;ZhuandJekeli,2006;蒋甫玉等,2011;Wangetal.,2012;Jiangetal.,2012;刘金钊等,2015)或者二者的联合(RózsaandTóth,2005;ZhuandJekeli,2007,2009;钱东等,2011;刘金钊等,2013;Buchaetal.,2016),利用频谱域或者空间域方法来实现.本文通过联合全球重力位模型(EGM2008)、航空重力扰动数据和剩余地形模型(RTM)数据,采用不同的解算方案,分别基于频谱域(二维FFT变换)和空间域(Stokes数值积分)算法,对毛乌素测区GT-2A航空重力测量系统采集的空中测线后处理重力扰动数据进行解算,构建了该地区的航空重力梯度扰动全张量.并对两种方案的解算结果进行了对比分析.1、数据和方法2015年4月17日至6月11日,自然资源部第一大地测量队采用引进的国内首台GT-2A航空重力仪在毛乌素地区实施了航空重力测量(宛家宽等,2017).毛乌素测区位于陕西省与内蒙古交界处,地貌以沙漠丘陵为主,地形起伏如图1,高程统计见表2.航空重力测量载体飞机选用塞斯纳208机型,机上搭载GT-2A航空重力测量系统和1台Ashtech双频GPS接收机,地面架设了2个GPS基准站.选取鄂尔多斯机场作为飞机起降保障机场,机场距毛乌素测区的飞行距离约为80km(蒋涛等,2018).图1中的黑色方框为本文计算结果的输出区域,该区域大小约139km×60km;灰色线条为测区飞行航线格网,航线格网间隔约0.55km,红色方框为考虑边缘效应内缩计算范围后的结果输出区域,大小约137km×58km.解算区域外地形质量影响结果分析见附录.本文航空重力梯度场构建技术路线见图2:首先在实测航空重力扰动中去除全球重力位模型EGM2008的重力扰动及剩余地形模型RTM的航空高度重力效应的影响,得到残余航空重力扰动数据,作为后续航空重力梯度场解算的输入数据之一;一方面,根据频谱域方法(二维FFT变换),沿方案一进行梯度解算;另一方面,将残余航空重力扰动向下延拓至大地水准面得到大地水准面上的残余重力扰动,由Stokes积分得到方案二梯度解算结果.具体解算过程及结果对比分析见后续章节.图1研究区域SRTM地形起伏图(黑色方框为结果输出区域,灰色线条为测区飞行航线格网,红色方框为考虑边缘效应内缩计算范围后的结果输出区域)图2本文航空重力梯度场构建技术路线图1.1航空重力扰动数据预处理实测航空重力扰动数据经插值后得到规则格网数据,这里选择的格网化插值方法是反距离加权插值,由于航空重力测线数据本身分布密集合理,简单的反距离加权插值即可达到较好的格网化效果,另外该方法对数据的处理较为简单,对数据的干预少,更能保持数据本身的特性(蒋涛等,2018).对85019个航空重力扰动数据进行格网化后,得到空间分辨率为0.3'×0.3'的重力扰动格网,相关统计量见表3.1.1.1全球重力位模型解算重力扰动基于EGM2008模型来解算研究区域的重力扰动,可以根据公式(1)来进行计算(Barthelmes,2013):是完全正则化的EGM2008扰动位球谐系数,Nmax是最大展开阶数,Pnm是m阶n次完全正则化勒让德函数.解算得到的EGM2008航空重力扰动相关统计量见表4.1.1.2剩余地形模型解算重力效应剩余地形模型RTM(ResidualTerrainModel,Forsberg,1984)通过SRTM高程数据(Jarvisetal.,2008)和DTM2006.0全球地形球谐展开模型(Pavlisetal.,2007)的差值来计算,得到的RTM数据根据公式(2)来求解相应的航空高度处的地形重力效应(Wang,1990):计算得到的RTM重力效应见表5.1.1.3残余航空重力扰动从格网化后的航空重力扰动数据中减去EGM2008得到航空高度处重力扰动及RTM的地形重力效应,就得到了残余航空重力扰动数据,见表6.作为本文利用不同解算方案构建区域梯度场的输入数据之一.1.2航空重力梯度扰动解算方案一根据1.1节得到的残余航空重力扰动数据,我们基于两种不同的解算方案(图2),分别利用频谱域算法和空间域算法求解最后的梯度张量.解算方案一选择频谱域方法,基于重力数据与重力梯度的频谱关系,通过二维傅里叶变换FFT来获得研究区域的残余航空重力梯度扰动,然后恢复EGM2008和RTM的梯度信号,最终获得完整的航空重力梯度场.1.2.1基于FFT方法计算航空高度残余梯度扰动基于残余航空重力数据,利用二维FFT方法解算得到残余航空梯度扰动见图3.1.2.2全球重力位模型解算重力梯度扰动.式中:GM是地心引力常数,r是航空高度到地心的径向距离,θ,λ为球坐标下的地心余纬和经度,CTnm和STnm是完全正则化的EGM2008扰动位球谐系数,Pnm是m阶n次完全正则化勒让德函数.解算得到的EGM2008航空重力梯度扰动相关统计量见表7.图3研究区域二维FFT方法解算的残余航空梯度扰动每个RTM点棱柱的尺寸定义为Δx,Δy及zRTM,解算得到的剩余地形模型的重力梯度效应见表8.向下延拓得到的研究区域大地水准面上残余航空重力扰动相关统计量见表9.利用该数据作为我们航空重力梯度扰动解算方案二的输入数据.同时,我们又将大地水准面上的残余航空重力扰动向上延拓回航空高度,与原残余航空重力扰动数据进行了对比,分析该算法的内符合精度,相关结果分析见第2节.基于向下延拓至大地水准面上的残余重力数据,利用Stokes数值积分方法解算得到残余航空梯度扰动见图4.与图3比较,我们发现,向下延拓过程放大了输入数据中的短波信号,突出了场源分布的细节特征.在恢复全球重力位模型的梯度扰动及剩余地形模型的重力梯度效应后,我们会对方案一、方案二得到的全张量梯度信号进行对比分析,相关结果见第2节.2、区域重力梯度标定场结果讨论2.1残余航空重力扰动延拓误差统计作为输入数据之一的残余航空重力扰动数据(原数据及其延拓值),其中解算方案一是用原残余航空重力扰动数据求解残余梯度扰动,解算方案二是用向下延拓至大地水准面上的残余重力扰动求解残余梯度扰动(见图2);向下延拓会增强原数据中的短波信号,突出细节特征,把向下延拓值再向上延拓回原航空高度,与原残余航空重力扰动进行比较,可以为我们分析解算得到的不同方案梯度扰动的定性分析提供参考(本文未就延拓过程对不同方案解算梯度扰动结果的影响进行定量分析,这部分工作将作为后续研究内容来开展).原残余航空重力扰动及其延拓值(先向下延拓至大地水准面,再向上延拓回航空高度),见图5,二者差值的统计量见表10.图4研究区域Stokes积分方法解算的残余航空梯度扰动从图5和表10中可以看出:残余航空重力扰动经向下延拓至大地水准面,再向上延拓至航空高度后与原数据差值的标准差为1.0078mGal.考虑边缘效应的影响,我们对原计算范围进行内缩,重新计算的差值统计量见表11,内缩计算范围后的差值标准差减小至0.1269mGal.2.2方案一与方案二航空重力梯度张量结果差值统计根据图2的航空重力梯度场构建技术路线,经二维FFT变换得到残余梯度扰动,恢复EGM2008解算的航空高度梯度扰动和RTM解算得到的航空高度梯度效应,最终得到方案一的梯度全张量结果,见图6.同理,由向下延拓至大地水准面上的重力扰动,经Stokes数值积分得到航空高度处的残余梯度扰动,再恢复EGM2008解算的航空高度梯度扰动和RTM解算得到的航空高度梯度效应,最终得到方案二的梯度全张量结果,见图7.图5残余航空重力扰动与其先向下再向上延拓后的值.(a)原残余航空重力扰动;(b)延拓后的残余航空重力扰动图6方案一计算航空高度梯度扰动张量对比图6和图7,我们发现两种方案解算得到的梯度各分量具有相同的变化趋势特征,所不同的是,由延拓值经Stokes数值积分得到的方案二结果表现出了更为明显的细节特征,这与我们在分析残余重力扰动延拓得到的结论是一致的.两种方案解算得到的研究区域航空高度梯度扰动张量的差值见表12,从表中可以看出两种方案得到的研究区域重力梯度扰动各分量之差的最大标准差为6.4798E(Γyz分量),最小标准差为2.6968E(Γxy分量);这里的计算结果不符值主要来源于两个方面:(1)由原残余重力扰动经延拓过程引入的误差;(2)由不同方案梯度解算算法引入的误差.这两类误差的定量统计超出了本文的研究范围,将作为后续工作进行深入分析.与残余重力扰动延拓对比分析一样,我们内缩了计算范围后发现(见表13):内缩计算范围得到的差值标准差最大值为1.8307E(Γzz分量),最小值为0.7223E(Γyz分量).目前国内没有实测重力梯度数据进行对比分析,然而由本文解算结果得到的内符合精度已经达到了2E以内,小于部分主流航空重力梯度测量系统相关技术指标(见表1),我们认为本文的解算结果是可靠的,这将为未来我国自主构建航空重力梯度标定场提供一定参考.图7方案二计算航空高度梯度扰动张量3、结论(1)本文通过联合全球重力位模型(EGM2008)、航空重力扰动数据和剩余地形模型(RTM)数据,基于频谱域(二维FFT变换)和空间域(Stokes数值积分)算法对毛乌素测区GT-2A航空重力测量系统采集的空中测线后处理重力扰动数据进行解算,构建了该地区的。