2. 中国自然资源航空物探遥感中心, 北京 100083;
3. 自然资源部航空地球物理与遥感地质重点实验室, 北京 100083
2. China Aero Geophysical Survey and Remote Sensing Center for Natural Resources, Beijing 100083, China;
3. Key Laboratory of Airborne Geophysics and Remote Sensing Geology, Ministry of Natural Resources, Beijing 100083, China
重、磁勘探以其经济快速、绿色环保,覆盖范围广、横向分辨率高等特点在解决地质问题中得到了广泛应用.为了更好地解决某一地质问题,往往在不同时期采集了不同尺度(比例尺、精度)、不同"维度"(卫星、航空、地面、海洋及井中等)的重、磁位场数据.然而,由于存在观测精度和观测位置及数据基准等差异,同一研究区采集的位场数据存在着精度差异、高度差异及数据基准差异等问题.如何将这些不同精度、不同比例尺、不同"维度"和不同基准的重、磁位场数据融合成为统一基准、同一"维度"(后简写为维度)的位场数据就成为重、磁位场数据融合的核心问题之一.根据位场数据融合中数据观测比例尺和观测维度的差异,本文将位场数据融合分为单尺度和多尺度、单维和多维融合问题.单尺度指相同网格距、相同比例尺和相同精度的位场数据,多尺度指不同网格距、不同比例尺或不同精度的位场数据;单维指同一个连续观测面的位场数据,多维指不同高度连续观测面或间断观测面的位场数据.
单尺度重、磁位场数据具有相同的观测位置和观测精度,不同来源的单尺度位场数仅存在数据基准的差异,其融合主要考虑数据基准的统一;而多尺度数据由于存在观测精度或观测比例尺的差异,需要在单尺度数据基准统一的基础上考虑精度校正.
单尺度数据融合方法主要有加权平均和主成分分析法等加权类方法(薛重生等,1997;张丽莉等,2003),其在没有基准差异的单尺度数据融合中计算效果良好,若数据基准存在较大差异,不考虑背景差校正的加权融合则不满足位场数据融合要求.多尺度数据融合方法主要分为两大类,第一类是插值方法,不同精度的数据同时进行插值计算,主要方法有Kriging插值、加权最小曲率插值、Shepard插值、样条插值等(Beylkin and Cramer, 2002;刘翔等,2016;刘善伟等,2018a;王虎彪等,2018;柯宝贵等,2018;支澳威和陈华根,2019),此类方法与单尺度数据中的加权计算方法思想类似,能够考虑不同精度数据的差异,但直接插值计算结果会受到多尺度数据基准差影响,不同精度数据边界处形成梯级带;第二类是由不同分辨率的图像融合方法延伸而来,包括基于小波变换(Ranchin and Wald, 1993;Deighan and Watts, 1997;郭达志和方涛,1997;Kuroishi and Keller, 2005;汪海洪,2005;李双成等,2006;Douma and De Hoop,2007;Herrmann et al., 2008;吴健生和刘苗,2008;王海江,2012;杨扬,2013;毕奔腾等,2016;谭沛森,2018;刘善伟等,2018b;陈勐韬等,2020)和Curvelet变换(Candes and Donoho, 2005;Nencinia et al., 2007;张恒磊和刘天佑,2011) 的融合方法,这类方法将数据分解为不同频段的信息,通过对多来源数据频段的选择和重新组合完成数据背景差的校正和不同精度数据的融合,其在不同分辨率的图像融合中应用效果良好,但由于不同精度的位场数据都包含全段频率的信息,进行频率分解再组合的方法会损失部分高精度数据细节信息,不符合位场数据优选高精度数据的原则.
单维数据不存在观测位置差异,其融合按照上述分类直接考虑单尺度或多尺度融合方法,而多维数据则需在单尺度融合和多尺度融合的基础上再考虑观测空间位置的差异.在实际的数据处理工作中,有学者不考虑观测位置的差异,将多维数据当作单维数据进行处理,这不符合位场数据随空间高度位置变化的规律.
多维数据融合可以通过位场延拓方法实现维度统一后,再进行上述单尺度和多尺度数据的融合计算(吴招才等,2018;余春平,2019);而目前针对多维数据的融合方法主要分为统计法和解析法两大类.统计法主要有最小二乘配置法,利用多维数据建立协方差函数后进行系统偏差估计和精度估计,并在此基础上逐渐发展了一些迭代算法、正则化方法等(Tscherning and Rapp, 1974;Forsberg,1987;杨元喜,1996;Tscherning et al., 1998;Strykowski and Forsberg, 1998;Tziavos et al., 1998;Olesen et al., 2002;Kern et al., 2003;Huang et al., 2006;翟振和和孙中苗,2010;刘晓刚等,2012;黄谟涛等,2013b;潘宗鹏等,2013;孙文等,2016;王灼华等,2017;赵启龙,2017),该方法能够将不同观测面上的多种类型数据融合起来,但需以足够分辨率的观测数据为基础,且经验协方差函数的选取复杂、大型矩阵的求解过程使得该方法适用性和稳定性受到限制.许多学者致力于解析法的研究(Lehmann,1993;Sideris,1996;Li and Sideris, 1997;宁津生等,1998;Kern et al., 2003;郝燕玲等,2007;Eicker,2008;Wittwer,2010;黄谟涛等, 2013a, 2015;王燚等,2015;吴怿昊等,2016;刘金钊等,2020),该方法在迭代中利用不同观测面上解析计算的位场数据残差修正观测数据或模型系数,引入分步平差、推估内插、迭代下延、点质量等技术,在大区域的重、磁力场计算中应用广泛,能够省去统计法中求解大型矩阵逆的问题,在迭代解析计算中完成多维多尺度位场数据的融合,但迭代中修正模型数据的方法无法全部保留高精度数据,且仍存在迭代稳定性的问题.
针对上述研究现状,本文提出了一套适用于多维多尺度重、磁位场数据融合的实用性方法;对数据维度、基准、尺度的差异,分别采用位场延拓、回归分析、精度加权平均等方法来解决高度不同、基准不同、精度不同或维度不同、尺度不同的重、磁位场数据融合问题,以期得到统一维度、基准及尺度的重、磁位场数据,为进一步处理、反演和解释奠定数据基础.
1 多维多尺度重、磁位场数据融合方法为了解决多维多尺度位场数据融合问题,需研究尺度统一方法、基准差统一方法(单尺度、多尺度位场数据融合方法)和维度统一方法(多维单尺度位场数据融合方法).
单尺度数据融合方法,采用相关分析和回归分析校正基准差解决数据基准差异;多尺度数据融合方法,采用不同精度数据加权解决观测尺度和观测精度的差异;多维多尺度位场数据融合方法,在以上两类数据融合的基础上采用位场延拓方法解决观测维度的差异(图 1).
本文采用位场延拓实现多维多尺度重、磁位场数据融合中的维度统一.本次所使用的位场延拓方法以空间域积分-迭代法(徐世浙,2006)为基础算法,将不稳定的向下延拓问题转换为迭代的稳定向上延拓问题,每次迭代中向上延拓和偏差叠加更新过程是影响计算精度和效率的关键,这两步骤也影响维度统一过程的计算精度和效率,因此成为进行多维数据融合的基础.基于此,本文采用滑动窗口加曲面切片技术(刘芬和王万银,2019)提高计算效率,采用向上延拓滤波技术提高计算稳定性.
(1) 位场向上延拓
空间域由原始点(ξ, η, ζ)延拓至计算点(x, y, z)的基本公式(曾华霖,2005)为:
(1) |
以上解析积分在程序实现中使用的是数值积分的矩形积分公式,则(1)式变为:
(2) |
其中,Δξ和Δη分别为x和y方向的点距,H(x, y, z, ξi, ηj, ζij)为核函数,其值如(3)式,
(3) |
本文采用滑动窗口和曲面切片的向上延拓计算技术(刘芬,2019)提高计算效率.根据位场延拓计算式(式(1))看出,距离越远的场源对计算点产生的异常影响越小.因此根据延拓高度在计算点一定窗口范围内计算,一般选择20倍的延拓高度以上,且小于剖面长度的一半(刘芬等,2019),可在保证精度的前提下减少计算量.
对于延拓高度相同的点(图 2中P点和Q点),在相同的窗口大小内,其核函数矩阵是相等的;平面到平面延拓(图 2由ΓA到ΓS),可以计算得到窗口内核函数矩阵,然后直接代入式(2)完成每个计算点的延拓.
若延拓面为曲面(图 2中ΓC),用一组平面切片Slice将其切割后,采用平面的滑动窗口核函数计算方法,首先计算得到每个切片高度上的固定窗口内的核函数矩阵,落在切片上的计算点(Q点和M点)直接代入,分布在切片之间的计算点(N点)的位场值UN则采用反距离加权插值得到,其加权公式(刘芬,2019)为:
(4) |
其中,UM、UQ分别为M点和Q点的位场值,hMN、hQN为N点分别与M点、Q点在垂直方向的距离.
由此,通过在空间域向上延拓计算中加入滑动窗口和曲面切片技术可在保证精度的前提下大大减少计算量,提高多维重、磁位场数据融合计算的效率.
设计理论模型检验向上延拓计算方法技术的正确性和有效性.场源由五个直立六面体组成,空间位置和物性参数(表 1)及其水平展布位置(图 3a)已知;模型基础观测曲面高度起伏为-1000~-2500 m(观测曲面一,z坐标向下为正,图 3b),观测平面为z=0 m,x方向坐标范围0~26000 m,y方向坐标范围0~26000 m,两个方向网格间距都为100 m,在观测曲面一和z=0 m平面上正演计算得到重力异常和磁力异常(图 3c、d、e、f).
本次测试z=0 m平面数据向上延拓至观测曲面一(-1000~-2500 m)的计算精度和效率.分别进行无窗口的全区域延拓计算、加滑动窗口的延拓计算以及滑动窗口结合曲面切片的延拓计算,分别统计了延拓计算至观测曲面一的重力异常(图 4a、b、c)和磁力异常(图 4d、e、f)与理论重、磁力异常(图 3c、d)的均方根误差、相对均方根误差以及计算时间(表 2,配置Intel(R) Core(i5) CPU 2.80 GHz,内存12.0 G的计算机).由表 2可以看出延拓计算的重、磁力异常与理论值的偏差都较小,滑动窗口和曲面切片的计算技术能够保证延拓计算的精度;无窗口的全区域的延拓计算耗时最长,加滑动窗口的延拓计算耗时次之,滑动窗口结合曲面切片的延拓计算耗时最短,验证了本次采用的滑动窗口结合曲面切片计算技术能够提高计算效率.这为迭代法向下延拓的多维重、磁位场数据融合中维度统一处理奠定了计算基础.
(2) 位场向下延拓
本文采用积分-迭代法(徐世浙,2006),将向下延拓转换为迭代法的向上延拓.如图 5所示,将观测面ΓB上的位场值UB(x, y, z)向下延拓至ΓA主要进行以下几个步骤:
S1:将ΓB上的位场值赋给ΓA,作为其初值UA(1)(x, y, z),即UA(1)(x, y, z)=UB(x, y, z);
S2:ΓA和ΓB之间无场源,利用空间域向上延拓公式(式(1)),由ΓA初值UA(1)(x, y, z)计算ΓB上的位场值UB(1)(x, y, z);
S3:计算UB(1)(x, y, z)与UA(1)(x, y, z)的偏差ΔUB(1)(x, y, z),直接加入UA(x, y, z)中进行校正更新,
(5) |
S4:当UB(x, y, z)和UB(n)(x, y, z)的残差不满足精度要求时,重复S2步和S3步,得到迭代公式,
(6) |
当|UB(x, y, z)-UB(n)(x, y, z)|<ε时,可以认为UA(x, y, z)≈UA(n+1)(x, y, z),其中ε是一个很小的数,即认为向下延拓求解完成.
若数据质量较高,采用该积分-迭代延拓可得到稳定延拓计算结果;若原始数据存在较大的噪声干扰,则在迭代过程中会存在随机误差逐步累积导致计算结果不稳定的问题,因此,本文提出对每次迭代中计算值与原始值的偏差采用向上延拓(式(1))进行低通滤波的技术提高计算精度.即在每一次迭代中对S3步的偏差值ΔUB(n)(x, y, z)进行低通滤波如(7)式,
(7) |
其中滤波器H就是向上延拓的核函数(式(3));相应地,迭代中的校正更新公式(式(6))变为:
(8) |
采用理论模型测试向下延拓方法技术的准确性和精度.观测曲面二与观测曲面一坐标范围一致,z坐标在观测曲面一的基础上抬升3000 m(z坐标从-3000 m到-4500 m变化,图 6),在其正演的重、磁力异常(图 7a、b)基础上加入各自数据幅值5%的随机噪声,即重力异常数据噪声均值为0 mGal、幅值为-0.049~0.055 mGal,磁力异常数据噪声均值为0 nT、幅值为-27.488~27.888 nT,利用加噪后的数据(图 7c、d)向下延拓至不同的平面(z=0 m、z=-1000 m、z=-2000 m、z=-3000 m),得到的向下延拓结果(图 8),计算的均方根误差和相对均方根误差见表 3.由图 8可以看出,向下延拓计算结果符合位场值随空间变化的规律,随着观测面高度与模型体距离的增加,重、磁力异常的幅值减小.由表 3可以看出,原始数据存在噪声时,进行低通滤波处理后向下延拓计算效果良好,各个延拓结果的计算误差小;且随着延拓距离的增加,向下延拓计算精度逐渐降低,这与位场随空间变化规律相符合.
同时,为检验上述位场延拓计算方法在多维重、磁位场数据融合中的应用效果,本文将向下延拓的计算结果(图 8)分别与原始理论重、磁力异常进行融合(图 9a、c),同时对比不进行位场延拓直接"融合"的结果(即将多维数据当作单维数据直接处理,图 9b、d),并统计融合计算误差(表 4).结果显示,在多维数据融合计算中,不考虑数据之间的维度差异直接进行"融合"是不正确,而本文采用的位场延拓方法可以有效统一不同数据的维度差异,计算误差小,是进行多维多尺度位场数据融合的维度统一的有效方法.
本文采用回归分析校正基准的方法技术解决单尺度数据和多尺度数据的基准差异问题.不同组数据间进行回归分析的前提是其存在显著性相关,故需对重叠区的多来源数据进行相关显著性检验.
在单尺度位场数据融合中,由于不存在精度或比例尺的差异,在理论上选取任一组数据作为基准数据都可,其余数据都向此基准数据校正;而对于存在精度或比例尺差异的多尺度数据,须选择精度最高或比例尺最大的一组数据作为校正的基准.
统计学领域中,相关分析是用来考查两组地位平等且无因果关系的变量相互变化的关联关系;回归分析是研究一个因变量与一个或多个自变量之间数量关系的统计方法.应用至位场数据融合中,相关分析中为两组重、磁位场数据,例如不同时期采集的同一地区的维度相同、尺度相同但存在数据基准差异的重力异常,记为gA、gB,采用Pearson相关计算公式,即(9)式,通过计算相关系数r判定两组数据之间的线性相关程度(Pearson, 1895):
(9) |
其中,gAi、gBi为两组数据在对应第i个数据点上的位场值,gA、gB为两组数据分别的平均值,n则为两组数据对应的点数.一般地,|r|越接近1,gA与gB的关系越密切,但严格意义上要作相关系数的显著性检验,如计算得到两组数据显著相关,则进行回归分析校正数据基准差.
不同来源的重、磁位场数据存在的数据基准差主要源自观测中基点不同,因此在建立回归模型过程中考虑构建较简单的回归模型,如一元线性回归或一元二次回归即可拟合不同数据基准的差异.求解方法是在建立合适的回归模型条件下,以最小二乘法求解模型统计数,获得回归方程,从而依据其中一组作为基准的数据对其他位场数据进行校正.
例如,对于N元线性回归,简单回归模型表示为(Galton,1886):
(10) |
模型中,gA是gB的函数部分加误差项,函数部分反映了数据gB与gA之间的变化对应关系,误差项ε为随机变量,反映了不同时期采集的位场数据gB和gA之间的线性关系之外的随机干扰,β0、β1、…、βN为回归模型的参数,由最小二乘估计得到.
用样本统计量
(11) |
其中,参数
(12) |
用最小二乘拟合的曲线来代表gA与gB之间的关系,与实际数据的误差比其他任何曲线都小,根据最小二乘法的要求,可求解
根据这一基本原理,进行多组数据之间的一元线性回归分析的处理是选择一个回归的数据体gA(x, y),依次用其他数据gB(x, y), gC(x, y), …, gM-1(x, y)与其进行回归校正,并建立简单线性回归模型形式为:
(13) |
即只用数据本身进行校正,得到每个数据体与选取的基准数据的回归系数后,进行相应的计算即完成了数据基准的回归分析校正.
选取单尺度理论模型数据测试该回归分析方法在消除基准差中的正确性和有效性.理论模型场源几何参数与物性参数见表 1,为模拟不同时期采集的存在数据基准差异的不同范围的单尺度重力、磁力异常数据,设计两个不同观测曲面(观测曲面一和观测曲面三),观测曲面三(图 10)是在观测曲面一(图 3b)的基础上存在空白区,两观测曲面的重复区占总范围的37%,设计观测面三原始观测异常(数据体2)和观测面一原始观测异常(数据体1)之间存在线性的基准差,其关系式如式(14)所示:
(14) |
对两组重力异常(图 11a、b)和磁力异常(图 11c、d)进行回归分析的基准校正,分别得到校正后重力异常和磁力异常(图 12a、b)和计算误差(图 12c、d),可以看到其校正后重、磁力异常形态和变化趋势符合位场数据的特点,且形态和幅值均与理论重力异常(图 3c)和磁力异常(图 3d)吻合很好,计算误差小;重力异常校正计算的均方根误差小于2.008×10-4mGal,相对均方根误差为0.0079%,磁力异常校正计算的均方根误差小于0.135 nT,相对均方根误差为0.0087%,计算精度高,说明本次采用的回归分析能够对单尺度位场数据基准差进行很好的校正.
多尺度理论模型场源几何参数与物性参数见表 1.为模拟多尺度数据特点,在观测曲面一的基础上,改变数据间隔为50 m并设置空白区形成观测曲面四(图 13a),改变数据间隔为200 m形成观测曲面五(图 13e),数据基准差关系式如式(14),同时加入不同大小的随机噪声模拟两者观测精度的差异,重力异常数据噪声均值为0 mGal、幅值为-0.028~0.027 mGal(图 13b)和-0.048~0.045 mGal(图 13c),磁力异常数据噪声均值为0 nT、幅值为-5.826~5.944 nT(图 13f)和-11.861~12.982 nT(图 13g);两组数据覆盖范围不同,存在数据重复区(数据重复区占总区域37%).
对重力异常、磁力异常分别进行回归分析校正得到原观测面上背景差校正后的重、磁力异常(图 13d、h),可以看到其与理论重力异常(图 3c)和磁力异常(图 3d)的形态和幅值均吻合很好,基准校正的效果受观测比例尺和随机误差的影响小,重力异常计算的均方根误差为0.013 mGal、相对均方根误差为0.433%,磁力异常计算的均方根误差为3.336 nT、相对均方根误差为0.175%,计算精度较高.同时,进行了不同空白范围的测试,即保持观测面五的范围不变,观测面四范围占总区域比例为9.46%、25.00%、37.86%、47.92%、71.59%,并进行误差统计(表 5),可以看出两组数据重叠区不同时,多尺度位场数据基准校正计算结果稳定性好,重叠区覆盖范围越广,校正计算精度越高.
尺度统一需考虑数据观测精度或比例尺的不同,在重叠区需根据精度或比例尺加权计算;为满足后期数据处理和解释的数据需求以及高精度数据的完全保留,需将多尺度数据统一至最高精度的位场数据尺度,本文采用最小曲率方法进行网格化(王万银和邱之云,2011).
重、磁位场数据是地下密度、磁性不均匀体的综合反映,高精度数据包含区域场信息和更丰富的细节信息,为在数据融合能保留高精度数据,在数据重叠区精度加权的原则为:高精度数据权重为1,低精度权重为0;若精度相同,则权重相同(即取平均值).需要注意的是,在进行位场延拓时,延拓的距离越大,计算的误差相应也越大,因此不同高度数据经位场延拓换算至同一个平面或曲面后可信度也有所差别,所以,多维单尺度数据重叠区也存在精度差异,需要加权处理.
根据多尺度位场数据融合的方法原理,选取已经完成维度统一和基准差统一的单尺度数据(图 12)和多尺度数据(图 13d、h),进行单尺度和多尺度融合(图 14).由图 14可以看出,位场数据融合计算结果能够保留高精度、大比例尺数据的细节信息,同时对研究区内空白区进行了有效补充,不同数据接边位置光滑,没有梯级带和突跳.同时,为测试重叠区范围对多尺度数据融合结果的影响,设计不同重叠区模型测试并统计了融合计算结果的误差(表 6),重叠区覆盖范围不同时,融合计算误差均很小,相对均方根误差均小于1%,计算稳定;且重叠区范围越广,融合精度越高.因此,本文提出的多尺度位场数据融合方法在数据重叠区范围不同时均可得到较好的融合计算结果.
最终设计多维多尺度理论模型测试其位场数据融合计算的正确性和有效性.场源参数如表 1,设计3个不同高度、不同尺度、数据基准有差异的不同覆盖范围的观测面,观测面一(图 3b),观测面六(在观测面一的基础上z坐标增大1000 m,z坐标范围为0~-1500 m,设置空白区),观测面七(在观测面一的基础上z坐标减小2000 m,z坐标范围为-3000~-4500 m;网格间距调整为200 m并设置空白区),以其理论重力异常(图 15b)和磁力异常(图 15c)进行融合测试.
采用本文的多维多尺度重、磁位场数据融合方法.首先将多维数据以维度统一方法计算至同一个观测面上变为单维数据,接着分别统一单维多尺度数据的基准差异和尺度差异,得到该观测面上的融合结果,利用位场延拓得到了z=0 m、z=-2000 m、z=-4000 m和z=-6000 m四个平面的重、磁力异常(图 16).可以看到其融合结果与各个计算面的理论重、磁力异常均吻合较好,重、磁力异常等值线光滑、异常变化趋势一致;随着延拓高度的增加,异常形态变化趋于平缓,符合位场数据随高度变化的规律;融合计算结果与理论异常的误差统计表(表 7)显示,其计算至不同高度平面结果均较稳定,计算精度较高,验证了该方法在多维多尺度重磁位场数据融合中的有效性.
为检验本文多维多尺度重、磁位场数据融合方法实用性,选取在中国某地航空和地面观测的不同比例尺重力异常数据进行多维多尺度位场数据融合试验.
2.1 单维多尺度重力数据融合试验选取中国某地1∶50000地面实测重力异常(数据体一,图 17b)和1 ∶200000地面实测重力异常数据(数据体二,图 17c)检验单维多尺度重、磁位场数据融合方法的实用效果.两个不同比例尺的地面数据的观测面均为地形面(图 17a),地形高度起伏变化范围是747~2623 m,两次重力观测得到了不同比例尺的数据体一(幅值-8~100 mGal)和数据体二(幅值-43~96 mGal),网格化后X和Y方向的网格间距都为150 m,其异常趋势和幅值变化在相同位置上大体一致,但1∶50000的数据反映更多细节信息,而1∶200000数据的变化趋势相对光滑.
利用单维多尺度重、磁位场数据融合方法将两组不同比例尺的重力数据进行融合,以高精度大比例尺数据为基准,将小比例尺数据进行数据基准差校正,相关系数为0.883,回归模型的斜率为0.996,截距(即数据基准差)为5.636 mGal;校正后的小比例尺数据与大比例尺数据进行偏差统计,其偏差平均值为0.017 mGal,均方根偏差为7.546 mGal,相对均方根偏差6.968%,其偏差较小能够保证校正后数据与高精度数据的一致性;之后进行数据融合得到单维多尺度重力数据融合结果(图 17d),融合结果重力异常的形态和幅值(33~100 mGal)均与原观测面基本一致,保留了地面高精度的大比例尺数据(图 17b)所包含的细节信息,同时用较低精度的小比例尺数据有效补充了高精度大比例尺数据的空白区,融合后的重力异常也符合位场的特点,验证了该融合方法在不同尺度实测重力异常的融合中的实用性.
2.2 多维多尺度重力数据融合试验为检验多维多尺度重、磁位场数据融合方法的实用性,在地面不同比例尺重力异常的基础上增加1∶50000航空实测重力异常数据,航空观测面GPS高度为1638~2705 m(图 18a),测线间隔为500 m,测点间隔为20 m,重力异常幅值为-43~96 mGal(图 18b),X和Y方向的网格化尺寸为150 m,航空重力异常与地面重力异常变化趋势大体一致,但总体幅值较低;两者比较看出,地面数据显示更多细节信息,航空数据观测比例尺虽与地面1∶50000重力数据一致,但其观测精度低于地面1∶50000数据,航空数据主要在覆盖范围和工作效率上更具有优势.
采用本文多维多尺度重、磁位场数据融合方法对数据体一(1∶50000地面重力异常,图 17b)、数据体二(1∶200000地面重力异常,图 17c)和数据体三(1∶50000航空重力异常,图 18b)三组航空和地面实测的不同精度重力异常数据进行融合试验.
数据体一和数据体三融合过程中,航空数据采用位场延拓计算至地面(维度统一)后进行基准差统一,与地面数据对比,其平均偏差小于10-5mGal,均方根偏差为7.379 mGal,相对均方根偏差为6.813%. 可以看出,数据体三计算至地面后的精度略高于数据体二(相对均方根偏差为6.968%),但精度相差不大;数据融合重力异常如图 19a所示,融合重力异常符合位场数据特征,形态和幅值均与原观测面基本一致,幅值范围为-44~124 mGal,保留了地面高精度数据(图 17b)包含的细节信息,同时用效率高、覆盖范围广的航空重力数据有效地补充了地面重力数据的空白区,两者数据相接的位置趋势一致且光滑无梯级带.
根据上述两组数据融合的精度比较结果,该地区三组多维多尺度数据融合时,首先将计算至地面后精度相差不大的航空1∶50000和地面1∶200000重力异常数据融合,加权平均的权重为0.5,得到初步融合结果(图 19b),再将其与精度最高的地面1∶50000重力数据进行多尺度融合,即为此三组多维多尺度数据融合重力异常(图 19c).
数据体二和三的融合结果(图 19c)已经结合了大范围的航空重力数据与小比例尺小范围地面数据的各自优势,经校正后异常形态趋势一致,最终融合结果(图 19c)相对于初步融合结果又包含更多地面高精度细节信息;统计初步融合结果校正至地面高精度数据后的偏差,平均值小于10-5mGal,均方根偏差为6.556 mGal,相对均方根偏差为6.054%,由此看出三组数据整体融合的效果优于任意两组数据,同样验证了该思路和方法在多维多尺度位场数据融合中的有效性.
另外,为了进一步验证和对比多维数据融合中维度统一的必要性和有效性,本文测试将多维数据当作单维数据不进行位场延拓直接"融合"(图 19d),其结果显示在不同数据的接边处存在明显梯级带,且异常值的变化趋势无法达到有效一致,不符合位场变化规律,由此也验证不考虑观测位置高度变化而将多维数据当作单维数据直接"融合"的做法是错误的.本文多维多尺度重、磁位场数据融合方法原理及计算技术对于不同观测比例尺的地面和航空数据进行了有效融合,验证了该方法在实测地面和航空数据融合中的实用性.
3 结论与建议本文研究了多维多尺度重、磁位场数据融合问题,采用改进的积分-迭代位场延拓方法实现维度统一,采用相关分析和回归分析实现数据基准统一,采用数据精度加权平均的方法实现尺度统一,解决了多维多尺度重、磁位场数据融合中维度、尺度和数据基准差异的问题,提出了一套适用于多维多尺度重、磁位场数据融合的方法.
理论模型测试和实测资料试验显示,该方法适用于航空、地面或海洋观测的重、磁位场数据融合,融合重、磁力异常能够作为后续位场数据处理、反演和解释的基础数据,具有良好的实用价值.本文提出的重、磁位场数据融合的思路也可推广至重、磁力异常各分量或张量数据的融合处理中.
致谢 感谢论文评审专家提出的修改意见和编辑对论文的编辑、加工.
Beylkin G, Cramer R. 2002. Toward multiresolution estimation and efficient representation of gravitational fields. Celestial Mechanics and Dynamical Astronomy, 84(1): 87-104. DOI:10.1023/A:1019941111529 |
Bi B T, Hu X Y, Li L Q, et al. 2016. Multi-scale analysis to the gravity field of the northeastern Tibetan plateau and its geodynamic implications. Chinese Journal of Geophysics (in Chinese), 59(2): 543-555. DOI:10.6038/cjg20160213 |
Candes E J, Donoho D L. 2005. Continuous curvelet transform: Ⅱ. discretization and frames. Applied and Computational Harmonic Analysis, 19(2): 198-222. |
Deighan A J, Watts D R. 1997. Ground-roll suppressing using the wavelet transform. Geophysics, 62(6): 1896-1903. DOI:10.1190/1.1444290 |
Douma H, De Hoop M V. 2007. Leading-order seismic imaging using curvelets. Geophysics, 72(6): S231-S248. DOI:10.1190/1.2785047 |
Eicker A. 2008. Gravity field refinement by radial basis functions from in-situ satellite data[Ph. D. thesis]. Bonn, Germany: Bonn University.
|
Forsberg R. 1987. A new covariance model for inertial gravimetry and gradiometry. Journal of Geophysical Research: Atmospheres, 92(B2): 1305-1310. DOI:10.1029/JB092iB02p01305 |
Galton F. 1886. Regression towards mediocrity in hereditary stature. The Journal of the Anthropological Institute of Great Britain and Ireland, 15: 246-263. DOI:10.2307/2841583 |
Guo D Z, Fang T. 1997. Space-scale adaptive fusion of multi-source image information using wavelet transform. Journal of China University of Mining & Technology (in Chinese), 26(4): 49-53. |
Hao Y L, Cheng Y, Liu F M, et al. 2007. Simulation of combination algorithm for heterogeneous marine gravity data. Journal of System Simulation (in Chinese), 19(21): 4897-4900. DOI:10.3969/j.issn.1004-731X.2007.21.013 |
Herrmann F J, Wang D L, Hennenfent G, et al. 2008. Curvelet-based seismic data processing: A multiscale and nonlinear approach. Geophysics, 73(1): A1-A5. |
Huang M T, Ouyang Y Z, Zhai G J, et al. 2013a. Analytical methods of multi-source gravity data fusion processing in the sea area. Geomatics and Information Science of Wuhan University (in Chinese), 38(11): 1261-1265. |
Huang M T, Ouyang Y Z, Zhai G J, et al. 2013b. Tikhonov Regularization collocation for multi-source gravity data fusion processing. Hydrographic Surveying and Charting (in Chinese), 33(3): 6-12. DOI:10.3969/j.issn.1671-3044.2013.03.002 |
Huang M T, Ouyang Y Z, Liu M, et al. 2015. Regularization of point-mass model for multi-source gravity data fusion processing. Geomatics and Information Science of Wuhan University (in Chinese), 40(2): 170-175. DOI:10.13203/j.whugis20130087 |
Hwang C, Guo J Y, Deng X L, et al. 2006. Coastal gravity anomalies from retracked geosat/GM altimetry: improvement, limitation and the role of airborne gravity data. Journal of Geodesy, 80(4): 204-216. DOI:10.1007/s00190-006-0052-x |
Ke B G, Zhang L M, Zhang C Y, et al. 2018. Fusion the altimetric and shipborne gravity data based on point mass fit method. Acta Geodaetica et Cartographica Sinica (in Chinese), 47(7): 924-929. DOI:10.11947/j.AGCS.2018.20170299 |
Kern M, Schwarz K P, Sneeuw N. 2003. A study on the combination of satellite, airborne and terrestrial gravity data. Journal of Geodesy, 77(3/4): 217-225. DOI:10.1007/s00190-003-0313-x |
Kuroishi Y, Keller W. 2005. Wavelet approach to improvement of gravity field-geoid modeling for Japan. Journal of Geophysical Research: Solid Earth, 110(B3): B03402. DOI:10.1029/2004JB003371 |
Lehmann R. 1993. The method of free-positioned point masses-geoid studies on the Gulf of Bothnia. Bulletin Géodésique, 67(1): 31-40. |
Li J W, Sideris M G. 1997. Marine gravity and geoid determination by optimal combination of satellite altimetry and shipborne gravimetry data. Journal of Geodesy, 71(4): 209-216. DOI:10.1007/s001900050088 |
Li S C, Gao W M, Zhou Q F, et al. 2006. Multi-scale spatial analysis on NDVI and topographical factors using wavelet transform. Acta Ecologica Sinica (in Chinese), 26(12): 4198-4203. DOI:10.3321/j.issn:1000-0933.2006.12.037 |
Liu F. 2019. Research on the indirect method of surface potential field data processing and transformation[Master's thesis] (in Chinese). Xi'an: Chang'an University.
|
Liu F, Wang W Y, Ji X L. 2019. Influence factors and stability analysis of plane potential field continuation in space and frequency domains. Geophysical and Geochemical Exploration (in Chinese), 43(2): 320-328. DOI:10.11720/wtyht.2019.1397 |
Liu J Z, Liang X H, Ye Z R, et al. 2020. Combining multi-source data to construct full tensor of regional airborne gravity gradient disturbance. Chinese Journal of Geophysics (in Chinese), 63(8): 3131-3143. DOI:10.6038/cjg2020O0044 |
Liu S W, Wan J H, Fu Y H. 2018a. Satellite altimetry gravity data and shipborne measurement gravity data fusion method: China, 107589464A. 2018a-01-16.
|
Liu S W, Kong X X, Wan J H. 2018b. Satellite height measurement gravity data and gravity satellite data fusion method: China, 107578068A. 2018b-01-12.
|
Liu X, Wang J X, Chen L, et al. 2016. A combination method of geomagnetism based on multi-surface function. Science of Surveying and Mapping (in Chinese), 41(2): 30-33. |
Liu X G, Pang Z X, Wu J. 2012. Earth's gravitational field model determination from different types of gravimetric data based on iteration method. Progress in Geophysics (in Chinese), 27(6): 2342-2347. DOI:10.6038/j.issn.1004-2903.2012.06.009 |
Nencinia F, Garzelli A, Baronti S, et al. 2007. Remote sensing image fusion using the curvelet transform. Information Fusion, 8(2): 143-156. DOI:10.1016/j.inffus.2006.02.001 |
Ning J S, Li J W, Chao D B. 1998. The spectral methods for computation of non isotropic local gravity field. Journal of Wuhan Technical University of Surveying and Mapping (in Chinese), 23(3): 227-229. |
Olesen A V, Andersen O B, Tscherning C C. 2002. Merging of airborne gravity and gravity derived from satellite altimetry: test cases along the coast of Greenland. Studia Geophysica et Geodaetica, 46(3): 387-394. DOI:10.1023/A:1019577232253 |
Pan Z P, Wang W, Zu A R. 2013. Fusion methods of multi-source gravity field information. Journal of Geomatics Science and Technology (in Chinese), 30(3): 236-240. DOI:10.3969/j.issn.1673-6338.2013.03.005 |
Pearson K. 1895. Note on regression and inheritance in the case of two parents. Proceedings of the Royal Society of London, 58(1): 240-242. |
Ranchin T, Wald L. 1993. The wavelet transform for the analysis of remotely sensed images. International Journal of Remote Sensing, 14(3): 615-619. DOI:10.1080/01431169308904362 |
Sideris M G. 1996. On the use of heterogeneous noisy data in spectral gravity field modeling methods. Journal of Geodesy, 70(8): 470-479. DOI:10.1007/BF00863619 |
Strykowski G, Forsberg R. 1998. Operational merging of satellite, airborne and surface gravity data by draping techniques. //Forsberg R, Feissel M, Dietrich R eds. International Association of Geodesy Symposia. Berlin, Heidelberg: Springer, 19: 243-248.
|
Sun W, Wu X P, Wang Q B, et al. 2016. Normalized collocation based on variance component estimate and its application in multi-source gravity data fusion. Geomatics and Information Science of Wuhan University (in Chinese), 41(8): 1087-1092. DOI:10.13203/j.whugis20140159 |
Tan P S. 2018. The multi-source gravity field fusion method and its application to the South China Sea. Geophysical and Geochemical Exploration (in Chinese), 42(1): 104-110. DOI:10.11720/wtyht.2018.1.12 |
Tscherning C C, Rapp R H. 1974. Closed covariance expressions for gravity anomalies, geoid undulations, and deflections of the vertical implied by anomaly degree variance models. Columbus: The Ohio State University.
|
Tscherning C C, Rubek F, Forsberg R. 1998. Combining airborne and ground gravity using collocation//Forsberg R, Feissel M, Dietrich R eds. Geodesy on the Move. Berlin, Heidelberg: Springer, 18-23.
|
Tziavos I N, Sideris M G, Forsberg R. 1998. Combined satellite altimetry and shipborne gravimetry data processing. Marine Geodesy, 21(4): 299-317. DOI:10.1080/01490419809388144 |
Wang H B, Xiao Y F, Wu L, et al. 2018. The fusion of gravity data and inversion of gravity vertical gradient anomaly. Hydrographic Surveying and Charting (in Chinese), 38(1): 1-4, 17. DOI:10.3969/j.issn.1671-3044.2018.01.001 |
Wang H H. 2005. Research on applications of wavelet multiscale analysis in the earth's gravity field[Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
|
Wang H J. 2012. Application of multi-resolution analysis methods to DEM multi-scale representation[Ph. D. thesis] (in Chinese). University of Chinese Academy of Sciences (Research Center of Soil and Water Conservation and Ecological Environment).
|
Wang W Y, Qiu Z Y. 2011. The research to a stable minimum curvature gridding method in potential data processing. Progress in Geophysics (in Chinese), 26(6): 2003-2010. DOI:10.3969/j.issn.1004-2903.2011.06.014 |
Wang Y, Jiang X D, Li D Y. 2015. The robust fusion of multi-source gravity data based on the spherical cap harmonic model. Acta Geodaetica et Cartographica Sinica (in Chinese), 44(9): 952-957. DOI:10.11947/j.AGCS.2015.20140345 |
Wang Z H, Jin H L, Fu F Y, et al. 2017. EGM2008 model and the measured ground gravity data fusion and its application in the western Sichuan region. Progress in Geophysics (in Chinese), 32(3): 1051-1058. DOI:10.6038/pg20170314 |
Wittwer T B. 2010. Regional gravity field modelling with radial basis functions[Ph. D. thesis]. The Netherlands: Delft University of Technology.
|
Wu J S, Liu M. 2008. Potential field data fusion based on wavelet multi-resolution decomposition. Journal of TongJi University (Natural Science) (in Chinese), 36(8): 1133-1137. DOI:10.3321/j.issn:0253-374X.2008.08.024 |
Wu Y H, Luo Z C, Zhou B Y. 2016. Regional gravity modeling based on heterogeneous data sets by using Poisson wavelets radial basis functions. Chinese Journal of Geophysics (in Chinese), 59(3): 852-864. DOI:10.6038/cjg20160308 |
Wu Z C, Gao J Y, Shen Z Y, et al. 2018. Magnetic anomaly in eastern China and adjacent sea, and its tectonic significance. Earth Science Frontiers (in Chinese), 25(1): 210-217. DOI:10.13745/j.esf.yx.2016-11-49 |
Xu S Z. 2006. The integral-iteration method for continuation of potential fields. Chinese Journal of Geophysics (in Chinese), 49(4): 1176-1182. |
Xue C S, Fu X L, Wang J M. 1997. Fusion processing of remote sensing and geophysical data and its application in geology-an example in Shangrao area. Geological Science and Technology Information (in Chinese), 16(S1): 36-43. |
Yang Y. 2013. Research on image fusion algorithms using multiscale analysis[Ph. D. thesis] (in Chinese). Beijing: University of Chinese Academy of Sciences (Changchun Institute of Optics, Fine Mechanics and Physics).
|
Yang Y X. 1996. Adaptively robust least squares estimation. Acta Geodaetica et Cartographica Sinica (in Chinese), 25(3): 206-211. |
Yu C P. 2019. Establishment of marine quasi-geoid model by multi-source gravity observation data. Geotechnical Investigation & Surveying (in Chinese), 47(9): 75-78. |
Zeng H L. 2005. Gravity Field and Gravity Exploration (in Chinese). Beijing: Geological Publishing House.
|
Zhai Z H, Sun Z M. 2010. The adaptive fusion of multi-source gravity data in Bohai Gulf. Acta Geodaetica et Cartographica Sinica (in Chinese), 39(5): 444-449. |
Zhang H L, Liu T Y. 2011. Potential field data fusion in curvelet domain. Oil Geophysical Prospecting (in Chinese), 46(4): 648-653. DOI:10.13810/j.cnki.issn.1000-7210.2011.04.006 |
Zhang L L, Wu J S, Wang J L. 2003. Application of image processing techniques to geophysics. Oil Geophysical Prospecting (in Chinese), 38(3): 317-323. DOI:10.3321/j.issn:1000-7210.2003.03.020 |
Zhao Q L. 2017. Data processing and application research of airborne gravimetry[Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.
|
Zhi A W, Chen H G. 2019. A fine fusion method for magnetic data of different platforms. Science of Surveying and Mapping (in Chinese), 44(12): 14-20, 34. DOI:10.16251/j.cnki.1009-2307.2019.12.003 |
毕奔腾, 胡祥云, 李丽清, 等. 2016. 青藏高原东北部多尺度重力场及其地球动力学意义. 地球物理学报, 59(2): 543-555. DOI:10.6038/cjg20160213 |
郭达志, 方涛. 1997. 多源影像信息的小波空间-尺度自适应融合. 中国矿业大学学报, 26(4): 49-53. |
郝燕玲, 成怡, 刘繁明, 等. 2007. 融合多类型海洋重力数据算法仿真研究. 系统仿真学报, 19(21): 4897-4900. DOI:10.3969/j.issn.1004-731X.2007.21.013 |
黄谟涛, 欧阳永忠, 翟国君, 等. 2013a. 海域多源重力数据融合处理的解析方法. 武汉大学学报·信息科学版, 38(11): 1261-1265. |
黄谟涛, 欧阳永忠, 翟国君, 等. 2013b. 融合多源重力数据的Tikhonov正则化配置法. 海洋测绘, 33(3): 6-12. DOI:10.3969/j.issn.1671-3044.2013.03.002 |
黄谟涛, 欧阳永忠, 刘敏, 等. 2015. 融合海域多源重力数据的正则化点质量方法. 武汉大学学报·信息科学版, 40(2): 170-175. DOI:10.13203/j.whugis20130087 |
柯宝贵, 张利明, 章传银, 等. 2018. 卫星测高与船载重力测量数据融合的点质量拟合法. 测绘学报, 47(7): 924-929. DOI:10.11947/j.AGCS.2018.20170299 |
李双成, 高伟明, 周巧富, 等. 2006. 基于小波变换的NDVI与地形因子多尺度空间相关分析. 生态学报, 26(12): 4198-4203. DOI:10.3321/j.issn:1000-0933.2006.12.037 |
刘芬. 2019. 间接的曲面位场数据处理与转换方法研究[硕士论文]. 西安: 长安大学.
|
刘芬, 王万银, 纪晓琳. 2019. 空间域和频率域平面位场延拓影响因素和稳定性分析. 物探与化探, 43(2): 320-328. DOI:10.11720/wtyht.2019.1397 |
刘金钊, 梁星辉, 叶周润, 等. 2020. 融合多源数据构建区域航空重力梯度扰动全张量. 地球物理学报, 63(8): 3131-3143. DOI:10.6038/cjg2020O0044 |
刘善伟, 万剑华, 符艺瀚. 2018a. 一种卫星测高重力数据与船测重力数据融合方法: 中国, 107589464A. 2018a-01-16.
|
刘善伟, 孔旭旭, 万剑华. 2018b. 一种卫星测高重力数据与重力卫星数据融合方法: 中国, 107578068A. 2018b-01-12.
|
刘翔, 王解先, 陈丽, 等. 2016. 一种多源地磁数据融合处理方法. 测绘科学, 41(2): 30-33. |
刘晓刚, 庞振兴, 吴娟. 2012. 联合不同类型重力测量数据确定地球重力场模型的迭代法. 地球物理学进展, 27(6): 2342-2347. DOI:10.6038/j.issn.1004-2903.2012.06.009 |
宁津生, 李金文, 晁定波. 1998. 各向异性局部重力场计算的谱方法. 武汉测绘科技大学学报, 23(3): 227-229. DOI:10.3321/j.issn:1671-8860.1998.03.010 |
潘宗鹏, 王伟, 祖安然. 2013. 多源重力场信息融合处理方法. 测绘科学技术学报, 30(3): 236-240. DOI:10.3969/j.issn.1673-6338.2013.03.005 |
孙文, 吴晓平, 王庆宾, 等. 2016. 基于方差分量估计的正则化配置法及其在多源重力数据融合中的应用. 武汉大学学报·信息科学版, 41(8): 1087-1092. DOI:10.13203/j.whugis20140159 |
谭沛森. 2018. 多源重力数据融合方法及其在南海的应用. 物探与化探, 42(1): 104-110. DOI:10.11720/wtyht.2018.1.12 |
王燚, 姜效典, 李德勇. 2015. 多源重力数据球冠谐模型抗差融合法. 测绘学报, 44(9): 952-957. DOI:10.11947/j.AGCS.2015.20140345 |
王虎彪, 肖耀飞, 武凛, 等. 2018. 重力数据融合与重力垂直梯度异常反演. 海洋测绘, 38(1): 1-4, 17. DOI:10.3969/j.issn.1671-3044.2018.01.001 |
汪海洪. 2005. 小波多尺度分析在地球重力场中的应用研究[博士论文]. 武汉: 武汉大学.
|
王海江. 2012. 多分辨率分析方法在DEM多尺度表达中的应用[博士论文]. 中国科学院研究生院(教育部水土保持与生态环境研究中心).
|
王万银, 邱之云. 2011. 一种稳定的位场数据最小曲率网格化方法研究. 地球物理学进展, 26(6): 2003-2010. DOI:10.3969/j.issn.1004-2903.2011.06.014 |
王灼华, 金红林, 付广裕, 等. 2017. EGM2008模型与地表观测重力数据的融合及其在川西地区的应用. 地球物理学进展, 32(3): 1051-1058. DOI:10.6038/pg20170314 |
吴健生, 刘苗. 2008. 基于小波的位场数据融合. 同济大学学报(自然科学版), 36(8): 1133-1137. DOI:10.3321/j.issn:0253-374X.2008.08.024 |
吴怿昊, 罗志才, 周波阳. 2016. 基于泊松小波径向基函数融合多源数据的局部重力场建模. 地球物理学报, 59(3): 852-864. DOI:10.6038/cjg20160308 |
吴招才, 高金耀, 沈中延, 等. 2018. 中国东部及近海磁异常特征及大地构造意义. 地学前缘, 25(1): 210-217. DOI:10.13745/j.esf.yx.2016-11-49 |
徐世浙. 2006. 位场延拓的积分-迭代法. 地球物理学报, 49(4): 1176-1182. DOI:10.3321/j.issn:0001-5733.2006.04.033 |
薛重生, 傅小林, 王京名. 1997. 遥感与地球物理数据的融合处理及其地质应用——以上饶地区为例. 地质科技情报, 16(S1): 36-43. |
杨扬. 2013. 基于多尺度分析的图像融合算法研究[博士论文]. 北京: 中国科学院研究生院(长春光学精密机械与物理研究所).
|
杨元喜. 1996. 自适应抗差最小二乘估计. 测绘学报, 25(3): 206-211. DOI:10.3321/j.issn:1001-1595.1996.03.009 |
余春平. 2019. 多源数据融合建立海域似大地水准面模型研究. 工程勘察, 47(9): 75-78. |
曾华霖. 2005. 重力场与重力勘探. 北京: 地质出版社.
|
翟振, 孙中苗. 2010. 渤海湾多源重力数据的自适应融合处理. 测绘学报, 39(5): 444-449. |
张恒磊, 刘天佑. 2011. 基于Curvelet域的位场多源数据融合方法. 石油地球物理勘探, 46(4): 648-653. DOI:10.13810/j.cnki.issn.1000-7210.2011.04.006 |
张丽莉, 吴健生, 王家林. 2003. 图像处理在地球物理学中的应用. 石油地球物理勘探, 38(3): 317-323. DOI:10.3321/j.issn:1000-7210.2003.03.020 |
赵启龙. 2017. 航空重力测量数据处理及应用研究[博士论文]. 武汉: 武汉大学.
|
支澳威, 陈华根. 2019. 一种不同平台磁力数据精细化融合方法. 测绘科学, 44(12): 14-20, 34. DOI:10.16251/j.cnki.1009-2307.2019.12.003 |