航空重(磁)多参量梯度探测与反演技术研究进展

马国庆,王君楠,孟庆发,孟兆海,秦朋波,王泰涵,李丽丽

1.吉林大学地球探测科学与技术学院,长春 130026

2.天津航海仪器研究所,天津 300131

3.广州海洋地质调查局,广州 510000

0 引言

近些年,我国资源勘探开始趋向大区域、深部勘探。而常规地面测量在交通不便、人迹罕至的区域难以实施,还存在测量精度不能满足要求、效率低等问题。随着卫星导航技术的逐渐成熟,航空重磁勘探开始投入实际应用。航空重磁勘探是指将测量仪器(重力仪、磁力仪等)和一些辅助设备装载在飞行器上,在测定区域上按照固定的测线对重力(梯度)数据或者磁(梯度、张量)数据进行测量的地球物理方法。与地面重磁勘探相比,航空重磁勘探具有特殊优势:1)在一些对于地面重磁勘探来说探测困难的地区(森林、山脉、沙漠、沼泽等),航空重磁勘探仍然可以获得高精度数据;2)对于较大的探测区域,航空重磁勘探方法效率更高;3)航空重磁勘探的覆盖区域更加完整。近些年,为了更加准确地描述地下构造,航空重磁勘探开始走向高精度、多参量测量发展,为地球多圈层、深部资源探测和国防等领域提供技术支撑[1-4]。

20世纪90年代,我国开始关注重磁多参量梯度探测技术的发展,各单位研制重力梯度仪和磁梯度仪[3, 5-10]。航空重磁多参量梯度探测将仪器搭载在飞行平台上,在测量过程中不可避免地会产生各种误差,影响后续的反演和解释结果,为了提升测量结果的精度,更多学者对航空重磁多参量梯度探测的数据处理展开研究[3-4]。针对单独反演的多解性问题,开始进行重力和磁数据联合反演。由于梯度数据有更高的水平分辨率,重磁及其梯度的协同反演可以更好地反映地下物性结构[11]。

基于航空重(磁)多参量梯度探测技术的研究背景,本文首先系统地总结了国内外探测装备的研究进展,然后针对航空重磁梯度和张量数据的特点,简单介绍了数据处理流程和反演方法,接着结合应用实例,分析了航空重(磁)多参量梯度探测技术的应用前景,最后提出展望,以期推动我国航空重磁梯度探测技术向实用化发展。

1 航空重磁多参量梯度探测技术

1.1 航空磁力探测

航空磁力探测是最早发展起来、应用最广泛的航空物探方法。航空磁力探测方法通常是在飞机的前方或者四周伸出一个探杆,探杆上装有磁力仪,在飞机内部搭载日变仪、卫星导航系统、飞行测高仪等辅助设备进行测量。最早的航空磁测出现在1936年,苏联使用了旋转线圈感应式航空磁力仪开展探测工作,但是其结果精度不满足要求[3]。二战期间,美国因为探测潜水艇的需要,发明了磁通门式航空磁力仪,灵敏度在1 nT左右,二战结束后,航空磁测开始应用于地质探测领域,包括地质调查和矿产探测等方面[3,6]。经过几十年的发展逐步发展成熟,航空磁测技术目前已经成功应用于商业勘探[3,5]。

航空磁测的测量对象主要包括:1)地球总磁场强度(T)及其梯度(Tx、Ty、Tz);2)磁感应强度三分量(Bx、By、Bz);3)磁全张量梯度(Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz)。总磁场强度梯度代表磁场模量的空间变化率,磁全张量梯度代表磁感应强度三分量的空间变化率(图1)。与传统的航磁总场测量相比,航磁梯度测量有以下优点:1)消除地磁日变影响,提高空间分辨率;2)提供更多的地下磁场信息,提高测量结果质量[5]。

图1 磁梯度及其张量梯度测量数据示意图

1.2 航空重力测量

航空重力测量是以飞机为载体,综合应用重力仪或者加速度计和辅助设备,包括惯性导航系统、卫星导航系统和测高、测姿设备来测定近地空重力加速度的测量方法[4]。航空重力勘探主要得到的测量数据有重力异常(gz)和重力张量梯度(gxx、gxy、gxz、gyx、gyy、gyz、gzx、gzy、gzz)。常规重力勘探得到的结果是重力位垂向一阶导数,用来表示重力异常,而重力梯度测量得到的是重力位的二阶导数,即重力张量,用来表示重力异常的变化率(图2)。

图2 重力梯度及其张量梯度测量数据示意图

航空重力梯度测量技术始于19世纪末,匈牙利物理学家厄缶发明了扭秤,重力梯度测量开始进入地球物理勘探邻域[7]。20世纪90年代后,航空重力梯度测量开始逐步实现商业化。随着仪器精度的不断提升,航空重力梯度测量技术因其高分辨率、高效率的优势快速发展。而重力全张量梯度测量最开始是为了满足美国海军的需要,绘制水下潜艇和导弹的位置。重力全张量梯度测量在20世纪90年代中期开始发展,随着冷战的结束,全张量梯度系统技术被允许商业化[8-9]。

重力梯度测量相比于重力测量,具有以下优点:1)重力梯度数据有更高的分辨率,可以反映地质体的细节;2)重力测量得到的是重力场的铅垂分量,重力梯度测量可以获得重力梯度张量中相互独立的5个参数,蕴含更多的地下信息;3)重力梯度测量有较强的抗动态干扰的能力,不易受到梯度仪搭载平台加速度的影响,可以在运动情况下进行测量,更适用于航空和海洋重力测量[10-11]。

2 航空重磁多参量梯度探测装备

2.1 航空磁梯度仪

按照发展历史,航空磁法探测技术可以分为3个阶段[12-15](表1):1)航磁总场测量阶段(20世纪40—70年代),学者们主要采用感应式磁力仪等磁力仪进行航磁总场测量(图3);2)航磁总场以及三分量测量阶段(20世纪70—90年代),此阶段开始出现航磁三分量测量,学者们一般采用光泵磁力仪测量地磁场的水平和垂直梯度(图4);3)航磁全张量梯度测量阶段(20世纪90年代至今),除了航磁三分量测量,学者们开始进行全张量梯度测量,一般采用磁通门磁力仪、高温超导磁力仪、低温超导磁力仪进行。

表1 航空磁测发展历程

据文献[12]修改。

据文献[12]修改。

2.1.1 国外发展

航磁梯度测量仪包括航磁三分量测量仪和航磁全张量梯度测量仪。航磁三分量测量系统包括吊舱式和硬架式两种。吊舱式航磁三分量测量系统包含三轴梯度测量系统和水平梯度测量系统,主要用来寻找与弱磁梯度异常相关的金、银等贵金属矿床,进行地质填图和工程地质测绘等工作。目前进入商业化阶段的三轴梯度测量系统包括3-axis AMG系统和Bluebird系统:3-axis AMG系统来自于加拿大Geotech公司,灵敏度在4 pT左右,绝对测量精度小于3 nT[16];Bluebird系统由加拿大GEM公司研制,灵敏度可以达到0.7 pT。硬架式航磁三分量测量系统有直升机硬架式航磁总场水平梯度测量系统和固定翼航磁总场梯度测量系统,前者适合在山地等复杂地形区使用,后者搭配于固定翼飞机,更适合远程的高海拔作业。

航磁全张量梯度测量系统的研究是当前的热点问题,相比于总磁场梯度测量,航磁全张量梯度测量的结果包含更多信息。航磁全张量梯度测量工作一般采用磁通门磁力仪和超导磁力仪(SQUID),美国、澳大利亚、德国、俄罗斯等国正在积极地研究和开发超导磁力仪[17-20]。超导磁力仪可以分为高温超导磁力仪(HTS-SQUID)和低温超导磁力仪(LTS-SQUID)两种,均使用液氮冷却,温度分别在77和4 K左右[9]。2004年,德国IPHT(institute for physical high-technology)首次采用低温超导磁力仪(图5)在南非进行了飞行实验[21-23];同年,美国Oak Ridge国家实验室和Tristan技术公司开始将高温超导磁力仪作为航空探测仪器。世界上第一台高温超导磁力仪设备于2010年被研发出来,澳大利亚CSIRO(Commonwealth Scientific and Industrial Research Organization)和五矿公司联合研制了航空张量磁梯度仪GETMAG,而后澳大利亚DSTO(Defence Science and Technology Organisation)和CSIRO联合研制了新型Mad系统,命名为MAGSAFE[23]。该设备装载的高温超导磁探头有较高的信噪比,可以测量极其微弱的磁场,而且具有较高的灵敏度,因此可以用来探测金刚石、铁矿等稀有金属矿床。该套系统被应用于2013年南非探测中,测得的梯度数据噪声水平达到10 pT/m[23]。

据文献[17]修改。

航空磁测的搭载平台分为有人机和无人机[24-28]。有人机飞行难度高、成本高、不适合低空工作;随着无人机技术的成熟,无人机(unmanned aerial vehicle, UAV)开始应用于航空磁测(图6)。近年来,无人机航磁测量系统因其安全性高、成本低、续航能力强等特点发展迅速。相比于有人机,无人机航磁测量可以保持极低的飞行高度,获得更高分辨率的磁及其梯度数据。许多国家都开展了无人机航磁测量设备技术的研发,并取得了显著成果[29]。

据文献[28]修改。

目前用于航空磁力测量的无人机类型包括固定翼无人机、无人直升机和多旋翼无人机等(图7),其中固定翼无人机和多旋翼无人机较为常用[30]。

a. 单电机固定翼;b. 单旋翼直升机;c. 四旋翼直升机;d. 六旋翼直升机。据文献[30]修改。

固定翼无人机具有续航时间长、速度快等优点,适用于大面积的快速测量,但是,它需要跑道才能起飞和降落,不适合低速和高分辨率测量;多旋翼无人机可以自动执行任务,并且易于操作,具有地形跟随功能,适用于小规模、高分辨率的测量,与中型和大型固定翼无人机相比,它相对便宜,但有效载荷能力差,飞行时间短;无人直升机可以垂直起降,可以根据任务需要改变飞行速度,适用于复杂地形或危险区域执行任务,但因其复杂的机械结构,发生故障时会带来更高的操作风险和维护成本。

英国Magsurvey公司在2003年首次进行了无人机航磁测量,开发了Prion无人机航空磁力测量系统[29]。此后,全球许多公司开发了基于无人机平台的航磁测量系统,如荷兰Fugre公司的GeoRanger-I[31]、加拿大UWG公司的Venturer[32]、日本YamahaMotor Co. Ltd公司的RMAX-G1[33]、瑞士和德国联合开发的Scout B1-100[32]、德国的MD4-1000[34]以及加拿大女王大学集成的GJI S900[35-36]等(表2)。

表2 各国无人机航磁测量系统

2.1.2 国内发展

建国初期,我国开始重视地学仪器研制事业,1959年我国第一个地质仪器专业生产厂——北京地质仪器厂成立[37]。1960年,长春地质学院开立了地质仪器系,研制了我国第一台航空核子旋进式磁力仪和光泵磁力仪。而后,学者们在1970年开始研发地球物理超导仪器,在1987年进行了高温超导磁力仪的研究[38]。

虽然我国的航磁全张量梯度测量起步较晚,与世界先进水平差距较大,但是也在逐步开展工作。目前国内有多家机构在进行超导磁力仪的研究,低温超导磁力仪的研制较多。在国家重大科研装备研制项目“深部资源探测核心装备研发”支持下,中国科学院上海微系统与信息研究所研制出了航空低温超导磁全张量测量系统,在内蒙古某区域进行了飞行试验,获得了全张量磁梯度分布图[39-40]。而高温超导磁力仪研制难度较高,国内研究较少,只有吉林大学和中国自然资源航空物探遥感中心组成的科研团队在开展工作,研发了吊舱式高温超导全张量磁梯度测量仪,该仪器适用于直升机,在江苏省某区域开展了飞行试验,获得了全张量磁梯度分布图,精度优于±30 pT/m[2, 39-40]。

对于无人机航磁测量,为了满足无人机对小型轻型航空磁力计的需求,Lu等[41]在2020年开发了国内首个功能性固定翼和多旋翼无人机航空地球物理测量系统——微型航空磁力系统(图8),并于2021年应用于CH-3无人机,成功完成了中国南方已知铜金矿床周围的实验测量,测量结果为矿区周边的探矿工作提供了直观有效的信息。微型航空磁力系统由全球导航卫星系统(global navigation satellite system,GNSS)天线、微型航空磁力系统(the miniature aerial magnetic system of IGGE (institute of geophysical and geochemical exploration), iMAMS)、遥控天线、雷达高度计、磁通门磁力仪等组成,其中iMAMS是一个自行设计的系统,包括以下主要组件:前端模拟组件、信号转换组件、全球导航卫星系统的同步和触发组件、主卡和电源组件。

2.2 航空重力梯度仪

2.2.1 国外发展

截止到现在,国际上正处于研究阶段的航空重力梯度仪可以分为旋转加速计重力梯度仪、超导重力梯度仪、冷原子重力梯度仪等[42-44]。

1970年,为了响应美国军方号召,美国Hughes公司、Draper实验室和Bell Aerospace公司参与了动态重力梯度仪的竞标,分别研制出旋转重力梯度仪、液浮重力梯度仪和旋转加速度计重力梯度仪,最终Bell Aerospace公司在竞争中胜出。从此旋转加速度计航空重力梯度仪成为近地表移动平台重力梯度测量系统(图9)。

a. 重力梯度仪及惯性平台;b. 全套重力梯度仪;c. 重力梯度仪航空测量平台。据文献[42]修改。

旋转加速度计重力梯度仪有4个加速度计,它们到中心的距离相等,位置相对的2个加速度计敏感轴方向相反,2对加速度计连线正交,圆盘以固定的角速度旋转。2个相对的加速度计输出求和可以消除共模线加速度的影响,再对输出的和信号求差可以消除角速度的影响,对重力梯度仪输出信号进行滤波和解调,就可以得到重力梯度分量[45]。

20世纪90年代,澳大利亚BHP Billiton公司和美国Lockheed Martin公司联合研制了基于加速度计的部分张量航空重力梯度测量系统FALCON(图10),其旋转轴接近垂直,并在1999年投入实际使用[46-47]。自重力梯度测量实际应用以来的第二个突破是FALCON系统的小型化。2005年,FALCON系统升级为数字电子设备,因此重量和尺寸大大减少,使FALCON系统成功地在小型直升机、单引擎和多引擎固定翼飞机(图10)上进行了常规飞行,测量结果的灵敏度和空间分辨率显著提高,对陡峭崎岖地形的适用性更高[46-48]。

a. 塞斯纳大篷车;b. Dash-7;c.Aerospatiale 350-B3直升机。据文献[46]修改。

2002年,Lockheed Martin公司研制了全张量重力梯度测量系统Air-FTG,而后美国Bell Aerospace公司在2003年首次将其用于商业勘探[48-49]。英国ARKeX公司于2005年从Lockheed Martin公司获得了重力梯度全张量(full tensor gradient,FTG)技术,研制了全张量重力梯度测量系统FTGeX[42, 50]。

20世纪90年代,来自美国斯坦福大学的Paik等人开始研制航空超导重力梯度仪。为了突破旋转加速计梯度仪的分辨率极限,马里兰大学(图11)和加拿大的Gedex公司等多家机构在2002年竞相研制航空超导重力梯度仪,然而过程并不顺利,迄今为止,各家机构都处于攻关阶段,航空超导重力梯度仪的研制仍具有广阔的发展前景[52-53]。

据文献[51]修改。

冷原子重力梯度仪在近些年迅速发展,它以冷原子干涉技术为核心,可以测量出重力梯度绝对值,具有高精度、无偏差、低漂移、自校准等优势[51-54]。

美国斯坦福大学于1991年首次实现冷原子物质波干涉,先后研制了冷原子重力仪和冷原子重力梯度仪。

2.2.2 国内发展

航空重力梯度仪遭到了国外全面的技术封锁,在国内发展较慢,在国家“十三五”被列为重点项目[55-56]。目前,航空重力梯度仪的研制效仿国外的路线,现在有多家单位参与研制,包括华中科技大学、天津航海仪器研究所和北京航天控制仪器研究所等。天津航海仪器研究所研发了基于旋转式加速度计的重力梯度仪样机,在2017—2021年对黑龙江哈尔滨、内蒙古鄂托克前旗、海南三亚、内蒙古通辽等地开展了航空重力测量试验(图12),实验结果达到国家先进水平[57-59]。随后华中科技大学在2022年也成功研制了旋转加速度计式重力梯度仪样机,实验室静态条件下达到了10 E分辨率[60]。

据文献[57]修改。

3 航空重磁多参量梯度探测数据处理技术

3.1 航空磁梯度数据处理技术

航磁全张量梯度数据处理流程主要包括噪声滤波、数据调平、坐标转换、航磁补偿、姿态校正等步骤(图13)。

图13 航磁全张量梯度数据处理流程图

航磁补偿是指去除地磁场和载体平台的磁干扰以及传感器误差的过程。磁力仪搭载在飞机上,飞机本身就是磁性体,不可避免地会对测量结果产生磁干扰。在航磁测量中,飞机自身姿态变化也会产生磁干扰。磁干扰的频率在磁力仪的测量范围内,不能通过滤波的方式解决,因此必须补偿掉这部分磁干扰的影响值。航磁补偿包含3个部分:1)飞行器上磁性物质的干扰;2)飞行器产生的感应磁场的干扰;3)飞行器闭合电路切割地磁场产生的涡流干扰[60-67]。

航磁总场的补偿公式可以表示为

Ht=H1+H2+H3。

(1)

式中:Ht为磁干扰的总强度;H1、H2、H3分别为恒定磁场、感应磁场、涡流磁场的强度。H1、H2、H3可以分别表示为

H1=c1cosX+c2cosY+c3cosZ;

(2)

(3)

(4)

式中:T为地球总磁场强度;ci为待定系数,i=1,2,…,18;X、Y、Z分别为地磁场矢量与飞机轴横向、纵向、垂向的夹角,可以通过惯导系统求得; (cosX)′、(cosY)′、(cosZ)′分别为夹角相对于基准面的微分算子。从航磁总场的补偿公式可以推得航磁三分量测量和航磁张量梯度测量的补偿公式。当确定ci后,得到补偿值,从测量值中减去该值,即完成飞机姿态补偿校正[68]。

航磁全张量梯度测量的数据处理不需要做日变改正。指航磁总场数据会受到周日变化、短周期扰动的变化磁场的影响,为了消除这种影响,航磁总场数据需要做日变改正[67-71]。全张量磁梯度数据几乎不受日变干扰的影响,所以不需要做日变改正。

数据调平是指磁场水平调整。在测量过程中,飞行器飞行高度和方向改变、起算点磁场选择误差等因素会对结果产生干扰,调平可以将数据值范围调整到正常水平,调平结果会直接影响反演和解释的精度。现在全张量磁梯度数据调平的相关研究较少,主要采用切割线调平、微调平和伪切割线调平等方法[69-71]。

坐标转换是指地理坐标系、地心坐标系、导航坐标系之间的互相转换,目的是将不同坐标系的数据转换为统一坐标系,方便后续计算[71]。

3.2 航空重力梯度数据处理技术

航空重力梯度测量是通过飞行平台上的重力梯度仪进行测量工作,测量得到的重力张量梯度结果会受到载体侧线误差、载体姿态误差、重力梯度仪观测误差以及载体自引力误差等误差的影响。为了减小测量结果误差,航空重力梯度数据处理主要包括数据检查、坐标转换、偏心校正、地形改正、数据滤波等[72-76](图14)。

图14 航空重力数据处理示意图

常用坐标系有实用惯性坐标系、地固坐标系、当地水平坐标系、载体坐标系等,它们之间可以进行坐标转换。如果采用GPS定位,需要将GPS WGS84坐标系的坐标数据转换到用户目标坐标系[77-83]。

在航空重力梯度测量过程中,载体的姿态很难稳定,三方向姿态角极易发生变化[84-88]。由于位于稳定平台上的重力梯度仪并不随着飞机飞行姿态的改变发生方向变化,因此在飞行过程中,飞行姿态改变和燃油消耗会导致梯度仪周围产生体积分布和质量分布的变化,使得重力梯度仪的周围梯度发生改变,严重影响重力梯度测量的精度和准确性。因此,需要对飞机的姿态变化和燃油消耗产生的梯度变化值进行计算和补偿,即偏心校正,使航空重力梯度仪的结果更接近实用化标准[1,89-92]。

航空重力梯度测量以高速飞行器作为测量支持,这意味着数据采集是在飞机移动时进行的,因此数据或多或少会受到各种外部因素的干扰,例如飞机的起伏、外部气流、内部电流不稳定等[93-94]。同时,在长时间测量中,传感器在温度等环境因素的影响下可能会发生变化。数据处理也会产生额外噪声,包括地形校正、调平和网格化的误差。因此,航空重力异常数据不可避免地包含大量噪声,这些噪声大多数是高频噪声。地球表面的重力梯度为3 500～4 500 E,而地下异常体引起的重力梯度异常大约为10 E,为了保证测量数据的可用性,重力梯度仪的噪声不能超过14 E/Hz,所以需要对重力梯度仪的输出信号进行滤波。目前采用的滤波方法有FIR(finite impulse response)数字带通滤波、卡尔曼滤波、克里金分析方法等[93-96]。FIR滤波器是航空重力梯度测量的数据处理常用的低通滤波器,可以用于去除重力梯度异常信号的所有高频信息,以便仅保留低频信号。FIR滤波器的原理可以用差分方程表示:

(5)

式中:x(n-k)为输入序列;y(n)为滤波器系列;h(k)为滤波器系数;N为滤波器长度。如果有用的重力信号和噪声之间没有明显的过渡,就很难消除重叠部分的噪声。为了保留更多有用信息,Zou等[92]在2016年加入卡尔曼滤波器,通过和FIR低通滤波器级联处理重力异常数据,提高了计算精度。传统的低通滤波器需要选择参数,而克里金分析方法对于地质特征明显的区域处理效果较好。赵德军等[96]在2021年将该方法首次应用于重力梯度数据降噪,明显提高了精度。

4 重(磁)原始及其张量梯度的多参量数据协同反演技术

4.1 磁及其梯度协同反演

磁总场数据可以很好地保留深源场的异常特征,磁梯度数据能够突出浅源地质体特征,而且具有更高的水平分辨率。金属矿产资源一般具有较高的磁异常响应,因此在金属资源勘察中常常采用磁及其梯度的协同反演方法。磁及其梯度协同反演主要包含2种方法:第一种是数据约束,将原始磁总场数据和磁梯度数据放在同一个矩阵中,通过改变权重来提升反演精度;第二种是结构约束,一般是在目标函数中加入磁梯度的交叉梯度项,可以提高结果分辨率[2, 97-98]。常规的磁总场和梯度协同反演的目标函数为

(6)

式中:dΔT为磁总场数据;dΔTx、dΔTy、dΔTz分别为x、y、z方向的磁总场梯度数据;GΔT、GΔTx、GΔTy、GΔTz为连接磁总场数据、磁总场梯度数据和物性参数的核函数矩阵;δ为正则化系数;κ为磁化率;W为磁总场及其梯度协同反演的模型加权矩阵。W可以表示为

(7)

对于磁及其梯度的结构约束,相当于在目标函数中加入自结构约束项,可以写为

(8)

4.2 重力及其梯度协同反演

航空重力勘探可以得到重力和重力梯度两种数据,虽然两者都反映地质体密度,但是重力数据反映深层信息,重力梯度数据反映浅层信息,因此将两种数据进行协同反演可以提高对地质体的分辨率。对于重力及其梯度的协同反演,现在多基于数据约束来实现协同反演,将两种数据放在同一个矩阵内[99]。

Zhadanov等[100]在2002年提出一种基于正则化聚焦反演的重力张量梯度反演方法。Capriotti等[101]在2014年将自适应加权函数作为目标函数来实现协同反演。2017年,Zhdanov等[102]提出了基于共轭梯度的正则化反演方法。而后张镕哲[98]在2021年实现了三维重力及其梯度协同反演方法。2019年,高秀鹤[2]等实现了基于协克里金的重力及其梯度协同反演方法。但是数据联合的方式会使得结果分辨率偏低,所以Ma等[103]在2021年推导了重力及其梯度协同密度反演方法(图15),并将其应用于庐枞实际地区。

据文献[104]修改。

基于交叉梯度的重力及其梯度协同反演方法的目标函数可以表示为

(9)

式中:m(1)和m(2)分别为重力及其梯度的密度参数;G(1)和G(2)分别为重力及其梯度的核函数矩阵;d(1)和d(2)分别为重力及其梯度异常;Wd和Wm分别为重力和重力梯度数据的权矩阵;λ为交叉梯度项系数,在每次迭代时进行优化;Φcross(m(1),m(2))为交叉梯度项。首先经过密度反演得到2个模型。然后采用取平均值、波数域转换等方法将两个模型进行数据融合。对于波数域方法,小波变换需要人为选择阶数,可以采用傅里叶变换方法。最后对融合模型进行正则化,也可以直接将两种模型取平均值。与传统数据组合联合反演方法相比,重力及其梯度的协同密度反演方法可同时获得更高分辨率、更精确的深层目标密度分布[104]。

5 航空重磁数据联合反演技术

对不同的航空地球物理测量数据进行反演,可以得到不同的地质信息,例如重力数据可以反映密度信息,磁性数据可以反映磁性信息;但是由于实际应用中观测点数远远少于剖分块数,会导致反演的多解性问题。针对地球物理数据单独反演的多解性问题,可以通过联合反演来解决。对比单独反演,联合反演的多种数据所包含的噪声也可以互补[105-107]。

联合反演根据下降方法可以分为同时联合反演和顺序联合反演。将两种观测数据的影响项加入到同一目标函数,共同下降,此方法称为同时联合反演(图16);也可以选择把观测数据分为两种目标函数,交替下降,此方法被称为顺序联合反演方法(图17)[108-110]。

据文献[108]修改。

据文献[108]修改。

以重磁两种观测数据为例,同时联合反演的目标函数公式为

(10)

式中:m1和m2为不同的物性参数;Φd1和Φd2为数据拟合项;Φm1和Φm2为模型约束项;Φj为联合反演的结构耦合项;λm1、λm2、λj1、λj2、λd1、λd2、λj分别为各项约束的所占比重。顺序联合反演公式为:

(11)

联合反演根据不同类观测数据对应的物理模型是否相同、相关,可以分为3类。首先,地球物理数据联合反演的实现必须保证不同观测数据对应的是同源场。最简单的联合反演就是不同地球物理数据反映的物性相同,例如重力及其梯度联合反演都是密度属性,磁及其梯度联合反演都是磁化强度属性。但是很多时候需要建立对应不同物性参数观测数据的联合反演,这类联合反演根据耦合方式可以分为2类:岩石物理耦合方式和结构耦合方式。岩石物理耦合的联合反演是指不同测量数据虽然对应不同的物性参数,但是可以建立数学相关性,这种岩石物理来自于经验公式、概率统计、理论关系;但是此方法有局限性,当地下结构过于复杂时,很难建立准确的物性之间的表达式。现在的联合反演大多数是结构耦合的联合反演,测量数据对应的物理模型不同,也没有一定的数学相关性,此时,由于不同的异常产生于同一个地质体,数据有结构相关性,例如基于交叉梯度约束和格拉姆约束的联合反演方法[108-111]。

5.1 基于交叉梯度约束的航空重磁联合反演

传统的联合反演方法需要提前知道不同物性参数的关系式,但是在实际应用中,不同物性参数的关系式往往过于复杂,难以写出,从而提出基于交叉梯度约束的航空重磁联合反演。交叉梯度约束通过使用交叉梯度项最小化参数函数来加强不同物理性质之间的结构相似性[112]。交叉梯度函数形式为

t(x,y,z)=∇m1(x,y,z)×∇m2(x,y,z)。

(12)

交叉梯度方法最初由Gallardo等[113]在2003年引入,用于电阻率断层扫描和地震数据的联合反演,并已成功应用于许多其他反演。第一个用于重力和磁性数据的交叉梯度联合反演算法由Fregoso等[114]在2009年提出,对于重力和总磁强度数据的交叉梯度联合反演,他们结合了奇异值分解(SVD)和其他常规正则化约束,以确定密度和磁化率的三维分布,并增强了结构相似性。而后Fregoso等[115]在2015年将欧拉反卷积方法引入重力和磁数据的交叉梯度联合反演中,并与单独反演结果进行比较,增加了反演结果的深部分辨率,将该方法成功应用于墨西哥Sebastian Vizcaino湾。Joulidehsar等[116]在2018年开发了一种增强的基于重力和总磁强度数据的联合反演算法,提高了计算效率。2019年Gross[117]介绍了在二维平面上的重力和磁异常数据的加权交叉梯度联合反演,正则化项是二阶相关函数,而交叉梯度是四阶相关函数,为了克服交叉梯度项相对于正则化项衰减更快从而导致失去相关性的问题,通过交叉梯度的局部加权,显著提高了结构相似性。为了提高反演效率,Zhang等[118]在2020年通过求解由数据加权项、正则化约束项和由结构耦合引起的交叉梯度项组成的目标函数,提出一种重力和磁数据快速结构联合反演算法,可以有效地提高异常位置反演的准确性。这些交叉梯度联合反演的研究使用不精确的结构相似性方法。

5.2 基于格拉姆约束的航空重磁联合反演

Zhdanov等[119]在2012年提出了基于格拉姆约束的正则化重磁联合反演方法,该方法不需要先验地了解不同模型参数和/或其属性之间的具体经验或统计关系,基本思想是通过格拉姆矩阵给出不同物性模型的结构相关性表达,再对基于不同物性模型的格拉姆行列式进行最小化[109, 120-121]。

Lin等[122]在2018年提出一种基于格拉姆约束的重磁数据联合反演方法,在密度和磁化率分布中提供统一的空间边界,有效确定异常体的形状、位置和物理性质。结果表明,该方法可以有效地应用于航空重磁数据的联合反演。格拉姆行列式表示为

(13)

式中:m(1)和m(2)为不同的物性向量;L为转换因子,可以将模型参数m(1)和m(2)转到其他空间,表示为Lm(1)和Lm(2);(·,·)代表向量内积。当L=1时,将物性参数放入格拉姆矩阵中,此时为物性格拉姆矩阵;当L=∇时,将物性求取梯度后放入矩阵,此时为物性梯度格拉姆矩阵[123-125]。

6 航空重磁及其多参量技术应用现状

随着社会需求和科学发展的不断增长,地下水、矿产和能源的近地表勘探前景正在迅速变化,几十年来对自然资源的不断探索迫使我们在更深的地下、海洋和泥火山周围寻找这些资源。目前对深层地下资源的认识还很少,特别是对自然资源的勘探[27,126],航空重磁多参量梯度探测可以帮助我们更好地了解深层构造。

航空重力梯度探测在地球物理勘探方面具有广泛的应用,比如可以用来识别地下构造、寻找矿产资源[127]。Zhdanov等[128]在2004年研究了一种基于正则化聚焦反演的张量重力场分量数据解释新方法,并应用于澳大利亚昆士兰州坎宁顿Ag-Pb-Zn矿体的重力梯度数据,与钻探结果的对比表明重力梯度测量可以显著提高矿产勘查的结果分辨率。Beiki[129]在2010年基于张量欧拉反卷积算法研究了南非弗里德堡的撞击结构,估计了源的位置。王婷一[110]在2021年推导了重力及其垂直梯度交叉约束联合密度结构反演方法,对辽源采空区应用此方法,找出了6个采空区,与已知采空区位置一致,同时对长白山的航空重力及其梯度数据进行反演,确立了岩浆囊的位置,与其他地球物理数据反演位置大概一致。针对不同类型的矿床,Dubey等[130]在2023年开发了一种基于粒子群算法的航空重力和重力梯度联合反演方法,分别应用于3个地区:美国路易斯安那州海上盐丘、加拿大魁北克省诺兰达矿区和瑞典韦斯特曼省的Karrbo地区,并验证了算法的有效性。乔中坤等[131]于2023年使用加拿大CG-6石英弹簧重力仪和美国的Burris重力仪在浙江工业大学某处后山进行实验,验证了重力梯度测量对地下空间探测的有效性。重力梯度测量也可以应用于计算大地水准面[132]、检测地下二氧化碳的分布,以及军事领域,如航空器和潜水艇的导航系统等。

对于航空磁梯度探测,磁梯度数据可以用于探测小型目标体,包括查明水下或者地下铁磁管线和电缆线路的分布、走向等[133],例如崔海波等[134]在2023年将磁梯度法应用于天津某深埋平行燃气管线探测,成功地定位了管线位置。航空磁及其梯度测量也可以用来寻找矿产,例如氧化铁-铜-金(FeO-Cu-Au)矿床、矽卡岩、块状硫化物和重矿砂等;还有用于定位有利的母岩或环境,例如碳酸盐岩、金伯利岩、斑岩侵入体、断层和热液蚀变,以及远景区域的一般地质测绘[41]。航空磁及其梯度测量也可以在军事领域中应用,磁性物体定位技术在航空飞机或水下航行器等自动化监视和安全系统中有着重要的应用,可以用于定位未引爆弹药的位置等[135]。Hu等[136]在2018年提出了一种基于磁场矢量及其梯度张量(MGT)的水下多物体定位方法,通过求解一组加权正则化的非线性方程,自动检测多个不同磁矩大小的运动物体的中心位置和多余的磁矩。航空磁测的另一个重要经济用途是绘制埋藏火成岩体的地图,另外,在区域勘探中,磁梯度探测对于了解构造环境非常重要[5,137-138]。

7 展望

航空重磁梯度、张量测量能更准确地描述地质体分布,未来重磁勘探会逐渐走向多参量化测量发展,因此在设备研制、数据处理、数据反演和勘探应用方面会迎来诸多挑战。

1)在设备研制方面:发展基于量子、超导、MEMS(micro electromechanical system)等新一代航空重磁多梯度测量装备,突破关键器件精度和灵敏度的限制;面向便携化测量,开发重力和磁张量梯度小型化、轻量化装备,从而有效搭载在无人机平台开展测量,扩展其应用范围。

2)在数据处理方面,针对运动环境下重磁张量梯度数据噪声频带范围大的特点,建立智能化数据处理方法,实现数据的高质量获取。

3)在数据反演方面,挖掘重磁张量梯度数据反演分辨率的数据潜力,突破深度与分辨率的限制,有效实现反演精度的提升。

4)在勘探应用方面,应用卫星重磁梯度、张量测量,面向深空探测,攻关卫星测量环境下重磁张量梯度测量设备以及智能化数据处理与反演算法。