CN101600975A - 重力勘测数据处理 - Google Patents
重力勘测数据处理 Download PDFInfo
- Publication number
- CN101600975A CN101600975A CNA2008800021674A CN200880002167A CN101600975A CN 101600975 A CN101600975 A CN 101600975A CN A2008800021674 A CNA2008800021674 A CN A2008800021674A CN 200880002167 A CN200880002167 A CN 200880002167A CN 101600975 A CN101600975 A CN 101600975A
- Authority
- CN
- China
- Prior art keywords
- data
- potential field
- model
- time
- measured potential
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 67
- 238000012545 processing Methods 0.000 title description 7
- 238000000034 method Methods 0.000 claims abstract description 80
- 238000005259 measurement Methods 0.000 claims abstract description 59
- 238000013507 mapping Methods 0.000 claims abstract description 38
- 230000008676 import Effects 0.000 claims abstract description 5
- 239000013598 vector Substances 0.000 claims description 28
- 238000003860 storage Methods 0.000 claims description 5
- 238000005516 engineering process Methods 0.000 abstract description 22
- 230000006872 improvement Effects 0.000 abstract description 6
- 238000004590 computer program Methods 0.000 abstract description 2
- 239000011159 matrix material Substances 0.000 description 23
- 230000006870 function Effects 0.000 description 22
- 230000008569 process Effects 0.000 description 18
- 238000005553 drilling Methods 0.000 description 12
- 230000008859 change Effects 0.000 description 10
- 238000009826 distribution Methods 0.000 description 9
- 238000004458 analytical method Methods 0.000 description 7
- 238000005457 optimization Methods 0.000 description 7
- 230000001133 acceleration Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 5
- 238000004613 tight binding model Methods 0.000 description 5
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 230000006399 behavior Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000036039 immunity Effects 0.000 description 2
- 238000012797 qualification Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000011282 treatment Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 241000271460 Crotalus cerastes Species 0.000 description 1
- 241001597008 Nomeidae Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000004836 empirical method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012821 model calculation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005381 potential energy Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000013517 stratification Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
- G01V7/02—Details
- G01V7/06—Analysis or interpretation of gravimetric records
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Environmental & Geological Engineering (AREA)
- Engineering & Computer Science (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
- Indication And Recording Devices For Special Purposes And Tariff Metering Devices (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及用于对来自机载勘测(如重力勘测)的势场勘测数据进行处理的改进技术,以及涉及针对这种技术的方法、设备和计算机程序代码。描述了一种对来自机载势场勘测或舰载势场勘测的、所测量的势场数据进行处理以确定用于映射场的场映射参数集合的方法,所述方法包括:输入所述所测量的势场数据,所述所测量的势场数据包括限定多个势场测量的数据,每个势场测量具有关联的测量位置和测量时间;以及使用模型来确定所述场映射参数集合,其中,所述模型包括空间部分与时间部分的组合,所述空间部分表示所述势场的空间变化,所述时间部分表示在所述所测量的势场数据中的时域噪声,所述确定包括使所述所测量的势场数据与所述模型的所述空间部分和所述时间部分都拟合。
Description
技术领域
本发明涉及用于处理来自诸如重力勘测等机载勘测的势场测量数据的改进技术,以及涉及针对这样的技术的方法、设备以及计算机程序代码。
背景技术
在本说明书中将涉及机载勘测,更具体地,涉及重力勘测。然而,所描述的技术不限于这些类型的勘测,而可用于其它势场勘测,所述其他势场勘测包括但不限于诸如大地电磁勘测、电磁勘测等磁场勘测。
通过测量势场数据来执行势场勘测,对于重力勘测所述势场数据可以包括一个或更多重力计数据(测量重力场)或重力梯度计数据(测量重力场梯度)、矢量磁力计数据、真实磁梯度计数据以及本领域技术人员熟知的其它类型的数据。地球物理势场勘测的共同目的是搜索潜在地指示有价值矿藏的特征。
传统地,诸如重力勘测等机载势场勘测是通过在网格图形上飞行而实现的。所述网格是由覆盖在下面地形上的二维表面上的正交平行线(飞行路径)集合来限定的。然而,所覆盖的表面受到允许飞行器飞行的离地最低点以及该飞行器爬升/下降的最大速率的限制。在本申请人共同悬而未决的PCT专利申请“Gravity Survey Data Processing”,PCT/GB2006/050211中描述了一些针对机载势场勘测的、便于接近地面采集数据的改进技术,其全部内容一并在此作为参考。
在采集了势场数据之后而在解释该数据之前,通常应用地形校正,以补偿表面高度变化。可以以数字地形海拔数据的形式置得表面轮廓,或根据(D)GPS((差分)全球定位系统)和/或诸如LIDAR(激光成像检测及测距)和SAR(合成孔径雷达)等的机载技术来确定表面轮廓。还可以使用飞行器的加速度、空间方位角、角速度和角加速度来校正势场仪器的输出数据。在于2006年1月25日提交的编号为0601482.3的共同悬而未决的UK专利申请“Terrain CorrectionSystems”中,描述了一些针对在地球物理勘测中的地形校正的改进技术,同样地其全部内容一并在此作为参考。
然而,不管先前描述过的改进技术如何,机载勘测测量中的低频噪声或漂移仍是个问题。(虽然术语漂移有时指的是随机单调变化,然而在本说明书中术语漂移用于表示特性频率比感兴趣的主信号小得多的、任何形式的相关噪声)。
在本领域中将术语“水准测量”用作通称以覆盖传统的数据修整技术。这些技术包括:去除低频漂移、对相邻线的低频内容进行匹配以及参考数据到固定高度的平面。例如,标准网格式勘测的交点可以用于交叉水准测量,其中对沿着勘测线的数据进行调整以使这些点的差异最小化。
然而,仍需要改进的数据处理,以及具体地需要改进的噪声处理。
发明内容
因此,根据本发明的第一方面,提供了一种对来自机载势场勘测或舰载势场勘测的、所测量的势场数据进行处理以确定用于映射场的场映射参数集合的方法,所述方法包括:输入所述所测量的势场数据,所述所测量的势场数据包括限定多个势场测量的数据,每个势场测量具有关联的测量位置和测量时间;以及使用模型来确定所述场映射参数集合,其中,所述模型包括空间部分与时间部分的组合,所述空间部分表示所述势场的空间变化,所述时间部分表示在所述所测量的势场数据中的时域噪声,所述确定包括使所述所测量的势场数据与所述模型的所述空间部分和所述时间部分都拟合。
在一些优选实施例中,模型的空间部分和时间部分包括模型方程的空间项和时间项,其中所述空间项和所述时间项被联合地确定,以确定场映射参数。优选地,所述模型具有线性形式;这简化了本发明。具体地,模型可以具有如下形式:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中,ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,f包括从模型估计出的测量的正演计算矢量。
优选地,模型针对测量之间的间隔具有充分高的时间分辨率,该时间分辨率小于仪器或测量噪声的预期相关时间。因而优选地响应于所测量的或预期的、所述测量数据中噪声的特性时间或相关时间,来选择节点之间的时间;这可从测量仪器上生产商的数据表单中获得。
在一些优选实施例中,模型的空间部分满足Laplace方程。优选地,模型的空间部分包括等效源模型。优选地,模型的时间部分包括由时间上的离散点(节点)所限定的线。
在实施例中,该线包括将这些节点相连接的逐片线性部分集合;在其他实施例中,可以使用高阶内插函数。优选地,节点在至少60、120、180、240、300、360或420秒的时间间隔处。优选地,噪声模型表示测量漂移(不需要是单调的)。
优选地,确定场映射参数包括:对联立方程系求解。由于等式可能是病态的,所以优选地使用正则化,然而优选地分别对模型的空间部分和时间部分进行正则化,以便于区分漂移与真正空间变化。因而,在一些优选实施例中,针对模型的空间部分使用梯度正则化,而针对模型的时间部分使用另一种类型的正则化,例如Tikhonov正则化。例如,这有助于区分真正空间变化与仪器漂移-比如可以在相对的方向上飞行的相邻勘测线中的不同空间方向上。应该意识到,在所描述的本技术的实施例中,通过同时使用来自多个不同勘测线的势场测量的数据与模型拟合,来确定场映射参数。
因而,在实施例中,拟合包括:对来自(机载)势场勘测的多个不同勘测线的、沿着这些线的长度而延伸的势场测量的数据进行联合拟合。这些勘测线可以实质上平行,或实质上不平行,例如实质上正交。在实施例中,勘测线包括多于两条勘测飞行线,优选地实质上包括所有的勘测飞行线-换而言之,在优选实施例中,实质上将从勘测中采集到的所有数据都提供给模型,其中所述模型针对完全勘测而对空间变化和时间变化进行建模。(在后一种情况中,勘测优选地包括多条飞行线,这些飞行线被布置成覆盖所勘测的区域,这些线可以进行任何合理的物理延伸。)
优选地,所测量的势场数据包括重力数据,如重力计数据(测量重力场)和/或重力梯度计数据(测量重力场梯度),并且可以包括多个分量。例如,所测量的势场数据可以包括重力梯度分量Gzz。因而,在实施例中,所映射的场包括重力场,尽管这可以由重力梯度数据来表示。然而,在其他实施例中,所测量的势场数据可以包括矢量磁力计数据、真实磁梯度计数据或其它势场数据(例如,标量重力测定数据)。
优选地,场映射参数包括可以描述势场特征以及具体地在实施例中可以用于根据正演计算而产生类似势场数据的参数或系数。
因而,在实施例中,参数值限定势场(例如,重力场或重力梯度场),其表示与所测量的势场数据的最佳拟合。优选地,此最佳拟合使残差最小化,所述残差依赖于在来自模型的势场的正演计算与所测量的势场数据之间的差异(例如均方差、模数差值或某种其他残差)。
在实施例中,场映射参数得自于模型的空间部分,并且可以包括针对等效源质量元的值,或更一般地针对具有关联的源强度的等效源元的值,在此情况下,场映射参数可包括这些源强度。
在实施例中,本方法还包括使用场映射参数集合来确定地图,例如在附加处理之后,这可以描述势场数据或从势场数据中获得的数据。这样的地图可以显示二维或三维数据,在实施例中,通过执行正演计算以根据场映射参数集合来确定例如表面上的势场,从而构造出所述地图。例如,可将场映射参数用在正演计算中,以使得可以确定和/或映射势场,如重力场或重力梯度场或关联的标量势。
技术人员将理解,通常,当测量势场数据时,实际测量的是重力和/或重力梯度(尽管附加地或备选地,从势场的空间倒数得到的其它量是可以被测量的)。
在使用等效源表示的本方法优选实施例中,在正演计算中使用优选地平面上(然而这不是必须的,因为可能选择在所勘测的区域中与地球表面近似的表面)多个表面质量元的表面质量(密度)值,以在平坦的映射平面或一些其它便利的表面上映射场的一个或更多分量。
在本发明的另一相关方面,提供了一种对所测量的势场数据进行处理以确定场映射数据的等效源方法,其中所述所测量的势场数据被建模成空间域信号与时域噪声分量的组合,该方法包括:确定所述所测量的势场数据与建模后的组合的联合最佳拟合。
本发明还提供了用于实现上述方法的处理器控制代码,具体地是在诸如盘、CD-ROM或DVD-ROM、已编程存储器(如只读存储器(防火墙))之类的数据载体上提供的,或者是在诸如光信号载体或电信号载体之类的数据载体上提供的。用于实现本发明的实施例的码(和/或数据)可以包括:采用传统编程语言(例如C)形式(解释或编译的)的源码、目标代码或可执行码,或汇编码,用于设置或控制ASIC(特定用途集成电路)或FPGA(现场可编程门阵列)的码、或用于硬件描述语言的码,如Verilog(商标)或VHDL(高速集成电路硬件描述语言)。本领域技术人员应意识到,这样的码和/或数据可分布在彼此通信的多个耦合的组件之间,例如沿着网络分布。
本发明还提供了一种用于实现上述方法的实施例的数据处理系统。
因此,在本发明的另一方面,提供了一种数据处理系统,用于对来自势场勘测的、所测量的势场数据进行处理,以确定用于映射场的场映射参数集合,所述系统包括:数据存储器,用于所述所测量的势场数据,所述所测量的势场数据包括限定多个势场测量的数据,每个势场测量具有关联的测量位置和测量时间;以及程序存储器,存储处理器控制码;以及处理器,耦合至所述数据存储器和所述程序存储器,用于载入并执行所述控制码,所述码包括用于控制处理器进行以下操作的码:输入所述所测量的势场数据;以及使用模型来确定所述场映射参数集合,其中,所述模型包括空间部分和时间部分的组合,所述空间部分表示所述势场的空间变化,所述时间部分表示在所述所测量的势场数据中的时域噪声,确定所述场映射参数集合的所述码被配置为使所述所测量的势场数据与所述模型的所述空间部分和所述时间部分都相拟合。
本发明还提供了一种承载数据结构的载体,所述数据结构包括对至少A和ρ的值加以限定的数据,其中,A和ρ被限定为使得针对势场测量集合m(x,y,z,t)将L(f,m)的值最小化,其中L(f,m)表示f与m之间的差异的量,其中:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,以及其中,所述数据结构所承载的数据与所述测量集合m(x,y,z,t)相关联或包括所述测量集合m(x,y,z,t)。
本发明还另外提供了一种承载数据结构的载体,所述数据结构限定地图,所述地图是根据对至少A和ρ的值加以限定的数据而得到的,其中,A和ρ被限定为使得针对势场测量集合m(x,y,z,t)将L(f,m)的值最小化,其中L(f,m)表示f和m之间的差异的量,其中:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,以及其中,所述数据结构所承载的数据与所述测量集合m(x,y,z,t)相关联或包括所述测量集合m(x,y,z,t)。
本发明还另外提供了一种如上以所述的承载利用前述方法而确定的映射数据集合的数据载体。
附图说明
现在将参考附图仅以示例的方式对本发明的这些和其他方面作进一步的描述,在附图中:
图1示出了根据本发明实施例的、针对联合等效源和时域漂移模型的示例时域内插函数;
图2示出了合成模型和仿真飞行图案;
图3a到3c示出了使用合成测量通过等效源反演而预测的Gzz的比较:a)没有漂移建模(仅空间等效源),b)根据本发明的实施例,采用同时漂移建模(扩充模型)以及c)没有向模型添加合成测量噪声(真实信号-无噪声);
图4示出了具有飞行勘测数据的飞行器,以及用于实现根据本发明的方法的实施例的数据处理系统的示例;以及
图5示出了对所测量的势场数据进行处理以实现根据本发明的方法实施例的过程的流程图。
具体实施方式
在实施例中,我们将描述使用联合等效源和针对所测量的势场数据的时域漂移建模的技术。
当提及场特别是重力场时,这并不限于矢量场,而是包括标量和张量场、势场以及得自于势场的任何衍生物。
势场数据包括但不限于重力计数据、重力梯度计数据、矢量磁力计数据和真实磁梯度计数据。可以从标量得到势场的元素和表示。
对于重力,相对势是重力标量势Φ(r),定义为:
其中r、ρ(r′)和G分别是重力场的测量位置、在位置r′处的质量密度、以及重力常数。经历重力场时的重力加速度是标量势的空间导数。由于重力具有方向性,所以重力是矢量。相对于任何所选的笛卡尔坐标系以三个分量将重力表示为:
这三个分量中的每一个分量在三个方向中的每一个方向上变化,因此而产生的九个量形成重力梯度张量:
通过对标量势函数、其导数、其傅里叶变换以及其它数学量的特性进行分析,得到基本方程以及势场的关系。
根据格林定理之一,如果在封闭表面上已知标量势的任何空间导数(包括标量势自身),则在该表面所封闭的体积内所有点已知该空间导数的值。对此的推论是,一旦通过微分和积分在所有点处知道该量,则可以获得标量势的所有其它导数,包括标量势自身在内。因此,当在封闭体积的表面上仅已知标量势的导数之一时,在该体积内所有点处可以有效地获知该标量势及其所有导数。这表明,对标量势的任何导数的任何分量的完全测量使得可以计算出该标量势的任何导数的任何其他分量。遵循这一点,至少在理论上,无论测量哪个量,都将仪器的选择简单归结为哪个仪器以最大的信噪比测量期望的量。
最终,上述重力标量势的微分最终产生:
其在没问题的区域中将重力的重要基本关系简化为Laplace方程:
调和函数满足Laplace方程,并且具有可以在对从势场勘测中采集到的数据进行分析中使用的很多特性。
可使用一系列技术对数据进行分析和处理,这些技术作为起始点利用从勘测中采集的数据来进行工作,随后改变数据和/或其格式,以使得与所勘测的量相关联的值都出现在水平固定海拔分析平面上规则的2-D网格上(水平校准和网格化)。
一般来说,在网格化过程中,将所勘测的区域划分矩形单元,所述矩形单元的边优选地与用于勘测的原则飞行方向对齐,然后将实际测量数据替换成以下这样的数据(网格化的数据):该数据“等效于”所测量的数据,但是现在分派有在每个单元中间点处的值。可以根据在两个正交的方向上飞行的线的平均间距来选择每个单元的大小。一旦数据是这种网格化格式的,该数据在数学上就更易于处理。可将数据视作数字组并且可以例如通过统计方法或其它方法来对该数据进行处理,以给出水平分析平面上势场的最佳估计。
先前描述了改进的等效源型方法(参见PCT/GB2006/050211,并入在此作为参考)。
等效源方法
作为背景对等效源建模方法的一些方面进行描述有助于理解本发明的实施例。
观察到,可以对物外面的重力场进行建模,使其如同来自整个位于物体表面并精确沿着物体表面的、薄得难以察觉的层中的物质。这样的层限定了二维等效源,即,产生实质上(确切地说是在理论上)与物体自身所产生的重力特征相同的的重力特征的重力源。存在很多可以限定等效源的方式,这些等效源可以是如上所述的或可以是严格水平的,可以全部或部分在地表之上或之下;可以是二维或三维的。然而,这些方式的一般特性是目的在于实质地产生与所勘测的地球部分所产生的重力场相同的重力场。关于其他信息,可以参考RJ.Blakely,″Potential Theory in Gravity and Magnetic Applications″,CambridgeUniversity Press,1995。
在勘测中,可以使用重力梯度计来测量作为位置rmeasure的函数的Gzz,以及利用Gzz来工作而无需产生重力梯度张量的其它元素。这可以用于产生对下面的质量分布的表示。以上提到的格林定理表明,原则上可以从Gzz得到下面的质量分布,尽管在后面描述的技术的实施例中不需要明确使用格林定理。
在示例等效源方法中,典型地以一侧50m的等级将勘测区域的表面分成小片(可以称其为小板或质量元)。可以对模型的分辨率进行选择以大致与勘测的分辨率相对应。很容易对来自每一个小板的重力进行正演计算(例如,参见R.J.Blakely,″Potential Theory in Gravity andMagnetic Applications″,Cambridge University Press,1995),调整该小板的质量直到获得与所勘测数据的最佳全面拟合为止。这种质量确定可使用标准的最小二乘法拟合过程。通过将真实测量位置处的数据与在相同的真实测量位置处由所提出的等效源产生的重力场相匹配,可以获得这种拟合。该过程在数学上是严格的并且不需要对数据进行任何人工调节以使该数据与水平矩形勘测相符。
一旦获得拟合,就认为该模型是主要数据集合。然后,用于确定地质结构的所有后续分析可以(但不是必需)对任何给出的地质结构将产生的重力场与由等效源产生的重力场之间的差进行比较并使该差异最小化。这种技术的一个显著优点是:最佳拟合来自质量分布(即使是合成的质量分布),从而该最佳拟合方案自动满足Laplace方程。这是对产生数字最佳拟合但并未强加以下附加限制的方法的改进:数据必须满足Laplace方程,即,数据可以来自实际质量分布。
等效源方法不是必须使用符合地形的表面,可以使用覆盖任何表面的源,所述表面可以在恒定海拔处也可以不在恒定海拔处,可以高于或低于地球的真实表面,可以与地球的真实表面相交也可以不相交等等。此外,等效源方法不是必须使用表面,而是可以使用三维模型。同样地,小板可以是二维或三维的任何大小或形状(或甚至可以是类似点的形状),小板实际上甚至不需要在大小和形状上相对应,如果允许小板的大小和几何结构根据每个区域中地形和/或地质变化的快慢而变化,则这对分析的效率是有帮助的。
因为目前可以合理地推断出遵循地形甚至遵循下面地质的等效源,所以该等效源可能产生更小的独立小板质量变化,然而原理上任何可行的选择都不对总体结果造成很大影响。因此,例如,如果怀疑出现具体地质层或者说地质异常(如角砾云橄岩管道),则可以限定等效源以将考虑该异常。同样地,可以在地质结构变化更详细或快速的区域加入更多的源。因而,在某些方面,可以认为等效源形状的选择类似于将工程产品离散化以进行有限元分析。
部分地,该过程的数学复杂度是由以下因素来确定的:所使用的源的数目,以及这些源当中有多少源用在了勘测过程中每个位置处的分析中。该技术的一个优点是:针对重力或重力梯度的一些分量,可以仅使用数据点区域内的源,这样可以实质上降低分析的复杂度。阈值距离典型地为几公里,例如在1到10公里的范围内,尽管阈值距离在某种程度上取决于地理(例如,可能需要扩展该距离以包围附近的大型山脉)。例如,勘测区域一侧的源元不可能对勘测区域另一侧的所测量的场做出显著贡献,并且在源元的影响比预期的或实际的噪声小得多的地方可以忽略这种贡献(设置为0)。这是有帮助的,因为例如这样的矩阵可能包括100K×500K的矩阵元素,如果这些元素大多是零,则显著减小了处理负担。
随后定义了以下方程:在该方程中源元对所测量的信号的贡献取决于由源影响值构成的矩阵并且(如在实践中所常见的)可以忽略很多源的影响,该矩阵变成稀疏矩阵以允许更高效的数字处理。
一旦产生了等效源,就可以通过直接正演计算来预测表面上重力标量势的任何导数。从分析和从可视化的观点来看,该过程都是有帮助的。因而,一旦发现质量元集合,正演计算(即,对质量元的影响求和)使得可以得到其它标量势的分量以及通过差分得到G的其他分量。可以将所推导出的G的值与地质模型相比较(称之为“解释”)以确定下面的地质结构。
更详细地,在给出每一个源元的质量(或密度)的情况下,使用直接正演计算来预测在每个勘测点针对所测量的量(比如重力矢量的分量或重力梯度的张量)可以获得什么样的值。通常,这将是以下所示形式的求和。此处使用gg作为所测量的量的标志,如上提到的,在某些优选实施例中该符号是Gzz。
在上述的方程中,F称为格林函数(同前,Blakely,第185页,并入在此作为参考),rmass-element限定质量元的位置(例如,重心或某一其他限定的点)。
函数F是标准函数,为本领域技术人员所公知(并且其应用不限于重力勘测)。本质上,函数F是单位质量或密度并且被限定了形状的源(质量元)在相关(测量)点处可以产生的影响。该源可以是点源、球或椭球,然而在技术的实际应用中该源更普遍地是可能不规则的棱柱。众多教科书针对简单的形状列出了格林函数;在该文献中可以发现针对更复杂的源几何结构的函数。此外,源影响进进行重叠,使得如果可以将复杂形状离散成多个简单形状,则可以将针对这些离散形状的格林函数相加在一起。原则上,这使得可以确定针对任何任意形状的格林函数的数值,即使实际上相对简单的形状一般是优选的。通过示例的方式,同前在Blakely第187页定义了针对矩形棱柱的格林函数F,其一并在此作为参考;该函数有8项,每一项与棱柱的顶点相对应。
等效源反演方法可以将势场分布的多个测量(势场的可观测量实际上是空间导数)组合成单个模型。该方法包括:寻找源模型的参数(例如,密度或几何结构),使得根据该模型计算出的数据与给定的测量集合相拟合。转换后的等效源模型然后可以再生测量中的信号,并且在极限之内可以用于将数据重新投射至不同的位置。通常,可调节的模型参数的数目小于独立测量的数目,因此,通过平均法则,从模型重新计算出的信号倾向于具有优于原始测量的信噪比。
等效源构造不仅仅是测量平均技术,因为该技术提出了以公式来表示空间基函数集合的一种方便且切实的方法,所述空间函数支持作为测量的基础的信号的函数形式。例如,基函数可以允许将不同类型的势场测量相组合;重力和重力梯度张量可以全部被转换成单个模型。基函数还可以倾向于使得转换后的模型对测量中某些错误变体的抗扰度提高,因为如果这些测量没有表现出真势场,则这些测量将不会被该模型的函数形式所支持。此外,通过排除已知在物理上不可实现的方案,等效源模型还可聚焦其函数空间。例如,置于地水准平面或其下的等效源模型将自然地过滤测量中非自然的高频变体(尽管原则上该模型可以被置于任何位置)。否则,在受约束较小的模型中,这样的变体可以产生大的错误。
测量的数目与模型的自由度之比的增大提高了反演的精确度,使得得到的等效源模型对测量噪声和干扰更不敏感。然而在很多情况下,由于与获取势能变体的完整采样相关联的实际困难,该比率通常不够有利,并且等效源反演可能变差。
机载勘测测量中的低频噪声(漂移)使得反演过程变差,因为这些现象可以类似于真信号(尽管术语漂移通常指的是随机单调变化,然而在本说明书中术语漂移用于表示特性频率比感兴趣的主信号小得多的、任何形式的相关噪声)。因此,等效源模型可以容易地容纳这样的变体,尤其在未确定该变体的地方。然后,通过后续等效源正演计算而得到的结果呈现了与测量图案相关的假特征。除此之外,在反转过程中噪声的相干分量没有最终到达到平衡,从而导致了劣等的模型参数估计。
因此,以下描述了等效源技术的扩充,其引入了容纳相关时域噪声的单独模型。通过使用单独模型(与等效源模型同时操作),可以使得测量中源参数估计更不易受到这种类型的噪声的影响。在反演之后,所扩充的模型可以单独预测所需的势场空间分布以及伴随测量的长波时域噪声。
联合等效源和时域漂移建模
可以将通过地球物理勘测在位置(x、y、z)时间t得到的测量分成一系列的项,
m(x,y,z,t)=S(x,y,z)+I(t)+C(t)+n (1)
其中,S(x,y,z)表示信号,I(t)表示干扰源,C(t)相关噪声,n纯粹随机噪声。在勘测之后进行处理的主要目的是以最佳精确度确定真实信号S(x,y,z)。通过对下面干扰和对应误差耦合传递函数的适当测量,可以校正干扰源I(t)。理论上,在这些校正之后,仅剩余的不需要的项是相关噪声C(t)和白噪声n。白噪声是基本的,仅可以通过多次勘测利用平均定律来减小白噪声的影响。然而,可以利用内插函数来对相关噪声中的大量功率进行建模,所述内差函数的时间周期与噪声的特性时间周期相似或比噪声的特性时间周期短。
图1示出了简单的时域内插函数,该时域内插函数是通过以规则间隔将线性部分逐片串联在一起而构造的。更具体地,图1示出了具有500秒的特性时间、通过利用逐片线性内插器对每400秒的节点(方形)进行连接而建模的指数相关噪声。
内插器在该范围上是连续的,然而在结点处具有不连续的导数。模型完全由节点处的值来指定,从而使自由度的数目等于节点的数目。在这种情况下,在最优化中(参见下文),节点值变成模型拟合参数。如果使用更复杂的内插函数,例如高阶多项式和样条,则自由度的数目会更多。如图1所示,使用具有比噪声的时间周期稍小的时间周期的基本内插器,可以很好地对相关噪声的长波分量进行建模。所描述的同时等效源和时域漂移模型使低阶漂移模型可以在测量集合内吸收大量时域相关噪声,使得等效源模型仅与空间相关信号相拟合。测量的其余未建模部分(残差拟合)主要是白噪声,这表明反演统计量与良好的模型参数估计相一致。
在线性反演模型的情况下,扩充的模型可被写成:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t) (2)
其中f是测量的模型正演计算矢量,ρ是等效源模型参数(例如,离散质量源单元的密度)的矢量,A是使源元的响应与测量位置相关的矩阵,λ是漂移模型参数(例如,逐片线性模型中的节点值)的矢量,并且B是使漂移模型内插与测量的次数相关的矩阵。以这种形式,在给出测量集合m(x,y,z,t)的情况下,可以利用使拟合中的残差最小化的任何标准最优化技术来确定模型参数ρ和λ,
minimise[L(f-m)] (3)
其中L表示对残差的量;例如L2范数使最优化成为最小二乘拟合。使均方误差最小化的L2范数具有易于计算的优点,并且有许多商业上可用的、用于实现最小二乘拟合数学包/库组件。一种变体是使用模(L1范数)来代替。这将降低野值(outlier)据点的重要性,从而使最优化更具鲁棒性;同样,合适的数学包/库组件在商业上是可用的,然而计算开销会更大(由于基础数学运算更难)。
在完全反演之后,独立使用等效源项来估计作为测量的基础的原始信号,
Aρ=S(x,y,z) (4)
或预测可能不同的位置处的不同量。类似地,可以独立使用(2)中的第二项来估计测量中的漂移,
Bλ=C(t) (5)
可以对(2)的模型进行扩展,以对来自所测量的势场数据的多个源的信号进行联合建模。在这种情况下,可以针对所测量的势场数据的不同源使用共同的等效源模型,然而,优选地单独的噪声(漂移)项包括在针对每个不同所测量数据源的模型中。例如,所测量的势场数据的不同源可以是对G的不同分量(Gxx、Gyy、Gzz、Gxy、Gxz、Gyz——因此可以有六个这种类型的不同设备)进行测量的不同设备和/或对g(gx、gy、gz——因此可以有三个这种类型的设备)进行测量的不同设备。例如,可以利用多个测量设备(尽管可以将这些设备组合在单个仪器中)以这种方式对重力和重力梯度势场数据的一个或更多个分量进行联合测量。可以采用与以上所述(并且下面做进一步讨论)的相类似的方式来执行最佳拟合最小化,尽管以更多的漂移参数来确定。
在这一点上,有助于详细说明矩阵A和B以及矢量ρ和λ。继续参见模型的空间部分,通常初始步骤是构建遵循勘测区域内地形的模型。将该模型离散成有限元集合,其中每个有限元限定质量元,所述质量元的密度是通过反演方案来确定的。当单位密度取决于质量元时,矩阵项Ai,j表示由测量位置i处的质量元j贡献的模型正演计算信号。例如,所使用的正演计算类型遵循测量信号类型;如果在相同位置得到测量n和m,然而测量n是重力梯度Gzz,测量m是重力梯度Gxx,那么来自质量元j的相应正演计算An,j和Am,j将相应地与Gzz和Gxx相对应。
矢量ρ对限定了模型的质量元的实际密度集合加以限定,并且在最优化之后包含方案的等效源部分。因此,矢量ρ有效地提供了映射场的参数集合。如上提到的,乘积Aρ的目的在于预测测量内的实际信号。
为了便于(2)的最优化,可以在构建矩阵A期间引入近似。这些近似可以采用舍位(将可略的矩阵项替换成零)或取平均的形式,通过取平均将若干独立元素的影响组合成表示平均影响的单个元素。当与源的大小相比测量点与源元之间的距离变得大的时候,上述两种近似都成为有效的。当使用这样的近似时,空间矩阵A变得稀疏,从而允许更有效的数字矩阵计算。
矩阵B的元素限定了将测量中的预测漂移与漂移模型参数λ相关的方程。对于逐片线性内插函数的情况(参见图1),如果在时刻ti得到落入在(节点)时间Tj和Tj+1处限定的节点漂移参数λj和λj+1之间的第i个测量,则针对测量中的漂移的线性内插预测如下:
依照矩阵等式(2)进行表达,上述等式如下限定矩阵项Bi,j和Bi,j+1:
漂移节点参数的时间Tj限定漂移模型的分辨率。当测量数据集合中存在中断时,例如在勘测线的结尾,或当出于某种原因暂停数据记录时,优选的是通过在恢复测量的时刻强加漂移参数来重启漂移模型。
当在相同勘测(例如,其中有若干仪器输出通道)中处理多个测量类型时,方便的是划分矩阵方程(2)使得每个测量类型占用矩阵A和B中的一组行。这样,λ中的第一组漂移参数表示针对第一测量类型的已建模漂移,第二组漂移参数表示针对第二测量类型的已建模漂移,以此类推。重要的是,应注意,虽然每个测量类型具有它自己的漂移模型,然而仅存在一个等效源模型。
如果在时域中残差的谱(f-m)不足够白,那么这指示漂移模型分辨率可能不合适,利用更短的时间周期或更高阶的内插可以得到好处。类似地,如果在空间域残差是有色的,那么这可以指示等效源模型被过度约束。如果在这些考虑之后,残差仍然示出可察觉的相关性,那么可能存在更基本的问题,例如测量中的过度干扰或设计得很差的等效源模型。
由于信号和漂移的性质有很大差别,一个在空间域相关,另一个在时域相关,因此在最优化期间自然要将这些量划分成(2)右侧的两项。为了说明这一点,考虑来自机载勘测的两条相邻线。虽然漂移沿着这两条独立的线高度相关,然而当在空间域中来看时这两条线之间没有相关性-这两条相邻线可能是相隔数周而得到的,或甚至是由不同仪器得到的。在(3)的最优化期间,由于这种类型的变化没有与等效源模型很好地拟合,所以要被正确分派给漂移模型项。相反,真信号通常出现在漂移模型的带宽外部,因此优选地被分派给等效源模型项。
然而,有多种情况,特别是在欠定的情况下,很难区分低频仪器噪声与真长波信号。如果沿着相邻线的漂移碰巧具有相同的一般走向,则噪声可以表现为空间相干的,并且在反演期间可以进入并破坏等效源项。相反,在时域中长波地质信号可以容易地作为漂移而出现,从而发现自己归因于(2)中错误的项。结果导致了在与不期望的噪声走向耦合的等效源预测(4)中的信号功率丢失。
为了加强对该问题的控制,可以针对扩充模型(2)的两部分引入单独的正则化。对于等效源密度分布ρ,梯度正则化是优选的,因为倘若模型的几何结构类似于物理上似是而非的某物,则现实中密度分布倾向于平滑变化。因此,优选地,选择更加平滑的密度方案将会提高本发明的完整性。对于漂移模型参数λ,标准Tikhonov正则化是优选的选择。在正则化下,(3)采取如下形式
其中,a是总正则化因子,b控制两个模型之间的相对正则化(实际中,由模型矩阵将正则化项归一化)。优选地,通过将测量中的多义性变化归因于密度变化而不是漂移,增大b将使得模型更加受等效源支配。针对正则化因子,没有用于推论出正则化因子正确值的决定性方法,然而“L-曲线”分析和其它经验方法可以给出指导,以确定在同测量的拟合与解决方案的行为之间的最佳正则化权衡。最终,要求进行一系列参数测试以确定a和b的最佳值。
因此,广泛地说,在一些优选的实施例中,在模型的空间部分中引入正则化,目的在于发现拟合测量的最平滑分布(优选地,在可能的情况下,通过设计初始源元模型来考虑地下密度的可疑的、大的变化)。广泛地,a的值决定了对模型强加平滑方案的程度,并且实际上该参数可以根据预期的下面地况而变化,以便得到看似良好的结果。然而,此外,为b选择大的值使得降低了漂移模型项的重要性(因为随后方程(8)的总体最小化需要较小的漂移模型项)。
然而,解决漂移与信号之间的多义性的最有效方法是确保在勘测中有足够勘测,这些勘测在空间上封闭而在时间上充分分离。例如,考虑勘测数据集合中两条线(通常是垂直的)交叉的点。在这些位置,信号S(x,y,z)对于任一线上的两个测量而言都是共同的。(如果一对测量之间的间隔与信号的相关宽度相比足够小,与勘测的高度相当,则可以将这对测量看作是有效的交叉点)。因此,在时间t1和t2处得到的测量之间的差异应该仅仅是由于相关噪声I(t2)-I(t1)加上一些随机噪声而得到的。在(2)的矩阵构造中,这些点对在等效源项中产生复制的行,而在漂移模型项中产生独立的行。因此,它们一起提供了以下方程:所述方程仅实质影响漂移方案,并且从而能够解决长波信号与漂移之间的多义性。通过说明信号的空间变化,同时等效源和漂移模型方案比交叉水准测量更加精确,这是因为漂移模型实质上是根据不同勘测线上多个位置处的勘测中的测量(潜在地,是所有的测量而不仅仅是交点处的测量)而推导出的。
为了说明组合的等效源与漂移模型的性能,提供了合成的示例。合成模型包括对5条不同大小的角砾云橄岩管道的选择,由遵循规则勘测飞行图案的机载重力梯度计来勘测这5条角砾云橄岩管道,如图2所示,图2示出了合成模型和仿真飞行图案。
在仿真中,将重力梯度分量Gzz正演计算成勘测图案的位置以表示由重力梯度计测量的信号。然后将指数相关时域噪声加到该合成信号以表示漂移以及可能在实际测量系统中出现的其他低频噪声。所仿真的测量用于首先在没有添加漂移模型的情况下对等效源进行反演,随后在添加漂移模型的情况下对等效源进行反演。在两种情况下,都使用反演后的模型来将Gzz正演计算成网格化的点集合,这些点与原始勘测具有相同的高度。图3示出了对于使用合成测量通过等效源反演而得到的、所预测的Gzz的比较:a)没有同时漂移建模,b)采用同时漂移建模。c)在没有向模型添加合成测量噪声的情况下得到的理想结果。
在上述示例中,为了说明当前的方法,没有对数据执行额外的降噪过程,尽管实际上还可以(一般是优选地)使用传统的降噪技术,如滤波/修整。从扩充的等效源模型中得到的结果(b)清楚地示出了鲁棒性级别以及对测量中相关噪声的抗扰度。在没有该扩充的情况下,结果(a)示出与勘测图案相关的严重伪迹以及对真实信号的几乎完全的屏蔽。
因此,以时域漂移模型来扩充等效源反演方法并随后同时针对两个模型的进行求解,这提供了一种对包含非常低频的噪声在内的势场数据集合进行可靠处理的方法。虽然通过传统的水准测量过程可以去除一定程度的低频噪声,然而例如在勘测不具足够的交叉点时这些方法是行不通的。将漂移模型与等效源模型相耦合的一个优点是:可以根据勘测中大量的或实质上所有的测量来推导出漂移行为。
现在参见图4,图4示出了飞行器10的示例,所述飞行器10用于执行势场勘测以获得用于依照上述方法进行处理的数据。飞行器10包括惯性平台12,在所述惯性平台12上安装有向数据采集系统16提供势场勘测数据的重力梯度计14(和/或矢量磁力计)。惯性平台12与惯性测量单元(IMU)18相拟合,所述惯性测量单元(IMU)18也向数据采集系统16提供数据,该数据典型地包括:姿态数据(例如,俯仰、侧滚以及偏航数据)、角速度和角加速度数据、以及飞行器加速度数据。飞行器还配备有差分GPS系统20和LIDAR系统22或类似物,以提供与下面地形上方飞行器的高度有关的数据。优选地从(D)GPS获得位置和时间数据,可选地,为了精确,与IMU相结合从(D)GPS获得位置和时间数据。
飞行器10还可以配备有同样向数据采集系统进行馈送的其它仪器24,如磁力计、TDEM(时域电磁系统)系统和/或超光谱成像系统。数据采集系统16还具有来自一般飞行器仪器26的输入,所述一般飞行器仪器26包括例如高度计、空气和/或地面速度数据等。例如,数据采集系统16可提供某种初始数据预处理,以校正针对飞行器运动的LIDAR数据以及将来自IMU 18和DGPS 20的数据相组合。可以为数据采集系统16提供通信链路16a和/或非易失性存储器16b,以使得可以存储采集到的势场和位置数据以供后续的处理。还可以提供网络接口(未示出)。
用于针对势场勘测产生地图数据的数据处理通常(但不是必须)是离线执行的,有时是在与采集勘测数据的所在地不同的地方执行的。如所示出的,数据处理系统50包括:耦合至码和数据存储器54的处理器52、输入/输出系统56(例如包括针对网络和/或存储介质和/或其它通信的接口)以及用户接口58(例如包括键盘和/或鼠标)。可以在可拆卸存储介质60上提供存储在存储器54中的码和/或数据。在依照以下所述的、图5所示的流程的实施例中,在操作中,所述数据包括从势场勘测采集到的数据,所述码包括用于对数据进行处理以产生地图数据的码,下面将进行描述。
参见图5,示出了用于在数据处理器上实现的过程的示例,在实施例中,所述数据处理器包括:通用计算机系统,用于依照前述技术对通过飞行勘测而得到的数据进行处理。因此,在步骤S200中,该过程输入所测量的势场数据,例如重力梯度计数据和关联的3D位置数据。可选地,在步骤S200a中,例如,可应用某种预处理来去除异常和/或减少(或增加)或选择要处理的数据。
在步骤S202中,该过程如上所述组成矩阵A和B,然后构造与所测量的势场数据(例如Gzz、3D测量位置以及时间数据)、空间场映射参数(ρ)和时域漂移参数(λ)有关的稀疏矩阵方程。可以特制所勘测的地形的等效源模型(S202a),或可以使用标准的规则网格。然后,该过程对矩阵方程求解(S206)以确定场映射参数集合,更具体地,等效源元强度集合,这些映射参数可以提供来自该过程的输出作为地图数据。备选地,可以通过使用这些参数执行正演计算,来作为该过程的一部分明确地计算出地图(S208)。例如,可以构造如图3所示类型的二维地图或如图2所示类型的三维地图(更准确地,地质模型)。
毫无疑问本领域普通技术人员可想到很多其他有效备选方案。例如,沿着例如在同上共同悬而未决PCT申请(PCT/GB2006/050211,并入在此作为参考)中所描述的线,可以在傅里叶域实现本文所描述的技术。
虽然使用机载势场勘测的优选示例描述了本技术,然而对于从船只得到的舰载势场勘测,更一般地,对于从其它移动平台或交通工具得到的势场勘测,也可以应用本实施例。
本技术不限于处理重力数据,例如,还可以用在处理磁场数据中。因此,例如,可使用由磁力梯度计所做的测量通过测量磁场和/或通量密度矢量和/或其幅度来获得所测量的势场数据。那么,等效源元可以具有例如表面电流密度或磁极强度。
应该理解,本发明不限于所述的实施例,本发明包括了对于本领域普通技术人员来说显而易见的、落入所附权利要求的精神和范围之内的修改。
Claims (29)
1.一种对来自机载势场勘测或舰载势场勘测的、所测量的势场数据进行处理以确定用于映射场的场映射参数集合的方法,所述方法包括:
输入所述所测量的势场数据,所述所测量的势场数据包括限定多个势场测量的数据,每个势场测量具有关联的测量位置和测量时间;以及
使用模型来确定所述场映射参数集合,其中,所述模型包括空间部分与时间部分的组合,所述空间部分表示所述势场的空间变化,所述时间部分表示在所述所测量的势场数据中的时域噪声,所述确定包括使所述所测量的势场数据与所述模型的所述空间部分和所述时间部分都拟合。
2.根据权利要求1所述的方法,其中,所述模型的所述空间部分和所述时间部分包括模型方程的空间项和时间项,所述确定包括联合地确定对所述空间项和所述时间项的估计,以确定所述场映射参数。
3.根据权利要求1或2所述的模型,其中,所述模型具有如下形式:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中,ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,f包括从模型估计出的测量的正演计算矢量。
4.根据权利要求1、2或3所述的方法,其中,所述模型的所述空间部分包括等效源模型。
5.根据权利要求1至4中任意一项所述的方法,其中,所述模型的所述时间部分包括内插函数,用于以大于60秒、120秒或300秒的时间间隔在节点之间进行内插。
6.根据权利要求1至5中任意一项所述的方法,其中,所述模型的所述时间部分包括逐片线性模型。
7.根据前述任意一项权利要求所述的方法,其中,所述噪声包括测量漂移。
8.根据前述任意一项权利要求所述的方法,其中,所述确定包括:将所述所测量的势场数据与利用所述场映射参数集合而预测的数据之间的差异的量最小化,以确定所述场映射参数集合。
9.根据前述任意一项权利要求所述的方法,其中,所述确定包括对所述模型的所述空间部分和所述时间部分进行正则化。
10.根据前述任意一项权利要求所述的方法,其中,所述确定包括:以第一正则化对所述模型的所述空间部分进行正则化,以及以不同的第二正则化对所述模型的所述时间部分进行正则化。
11.根据权利要求9或10所述的方法,包括针对所述模型的所述空间部分使用梯度正则化。
12.根据权利要求9、10或11所述的方法,包括针对所述模型的所述时间部分使用Tikhonov正则化。
13.根据前述任意一项权利要求所述的方法,其中,所述拟合包括对来自所述势场勘测中多条不同勘测线的势场测量的数据进行联合拟合。
14.根据前述任意一项权利要求所述的方法,其中,所述所测量的势场数据包括重力计数据和重力梯度计数据当中的一个或更多。
15.根据权利要求1至14中任意一项所述的方法,其中,所述场映射参数包括针对磁等效源元的值。
16.根据权利要求1至14中任意一项所述的方法,其中,所述场映射参数包括针对等效源质量元的值。
17.根据前述任意一项权利要求所述的方法,其中,所述所测量的势场数据包括来自多个势场测量源的数据,其中所述模型的所述空间部分对于所述多个测量源是共同的,以及其中针对每个所述势场测量源提供单独的、所述模型的所述时间部分。
18.根据前述任意一项权利要求所述的方法,还包括使用所述场映射参数来确定地图。
19.一种承载处理器控制码的载体,所述处理器控制码用于在运行时执行前述任意一项权利要求所述的方法。
20.一种数据处理系统,用于对来自势场勘测的、所测量的势场数据进行处理,以确定用于映射场的场映射参数集合,所述系统包括:
数据存储器,用于所述所测量的势场数据,所述所测量的势场数据包括限定多个势场测量的数据,每个势场测量具有关联的测量位置和测量时间;以及
程序存储器,存储处理器控制码;以及
处理器,耦合至所述数据存储器和所述程序存储器,用于载入并执行所述控制码,所述码包括用于控制处理器进行以下操作的码:
输入所述所测量的势场数据;以及
使用模型来确定所述场映射参数集合,其中,所述模型包括空间部分和时间部分的组合,所述空间部分表示所述势场的空间变化,所述时间部分表示在所述所测量的势场数据中的时域噪声,其中,确定所述场映射参数集合的所述码被配置为使所述所测量的势场数据与所述模型的所述空间部分和所述时间部分都相拟合。
21.一种对所测量的势场数据进行处理以确定场映射数据的等效源方法,其中所述所测量的势场数据被建模成空间域信号与时域噪声分量的组合,该方法包括:确定所述所测量的势场数据与建模后的组合的联合最佳拟合。
22.根据权利要求21所述的方法,还包括:将空间正则化和时间正则化之一或两者包括进建模后的、所述空间域信号和所述时域噪声的组合中。
23.根据权利要求21或22所述的方法,其中,所述所测量的势场数据包括重力数据。
24.一种承载处理器控制码的载体,所述处理器控制码用于在运行时执行权利要求21、22或23所述的方法。
25.一种数据处理系统,包括根据权利要求24所述的载体。
26.一种承载数据结构的载体,所述数据结构包括对至少A和ρ的值加以限定的数据,其中,A和ρ被限定为使得针对势场测量集合m(x,y,z,t)将L(f,m)的值最小化,其中L(f,m)表示f与m之间的差异的量,其中:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,以及其中,所述数据结构所承载的数据与所述测量集合m(x,y,z,t)相关联或包括所述测量集合m(x,y,z,t)。
27.一种承载数据结构的载体,所述数据结构限定地图,所述地图是根据对至少A和ρ的值加以限定的数据而得到的,其中,A和ρ被限定为使得针对势场测量集合m(x,y,z,t)将L(f,m)的值最小化,其中L(f,m)表示f和m之间的差异的量,其中:
f(x,y,z,t)=Aρ(x,y,z)+Bλ(t)
其中ρ和λ分别是空间模型参数的矢量和时间模型参数的矢量,A和B是矩阵,以及其中,所述数据结构所承载的数据与所述测量集合m(x,y,z,t)相关联或包括所述测量集合m(x,y,z,t)。
28.一种载体,承载利用权利要求1至17中任意一项所述的方法而确定的场映射参数集合。
29.一种载体,承载利用权利要求21至23中任意一项所述的方法而确定的映射数据。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
GB0701725.4 | 2007-01-30 | ||
GB0701725A GB2446174B (en) | 2007-01-30 | 2007-01-30 | Gravity survey data processing |
PCT/GB2008/050059 WO2008093139A2 (en) | 2007-01-30 | 2008-01-30 | Gravity survey data processing |
Publications (2)
Publication Number | Publication Date |
---|---|
CN101600975A true CN101600975A (zh) | 2009-12-09 |
CN101600975B CN101600975B (zh) | 2013-01-23 |
Family
ID=37873010
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN2008800021674A Expired - Fee Related CN101600975B (zh) | 2007-01-30 | 2008-01-30 | 重力勘测数据处理 |
Country Status (7)
Country | Link |
---|---|
US (1) | US8332184B2 (zh) |
CN (1) | CN101600975B (zh) |
AU (1) | AU2008211708B2 (zh) |
CA (1) | CA2674535A1 (zh) |
GB (1) | GB2446174B (zh) |
WO (1) | WO2008093139A2 (zh) |
ZA (1) | ZA200902158B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116794735A (zh) * | 2023-06-02 | 2023-09-22 | 中国自然资源航空物探遥感中心 | 航空磁矢量梯度数据等效源多分量联合去噪方法及装置 |
Families Citing this family (21)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB2447699B (en) | 2007-03-23 | 2011-07-13 | Arkex Ltd | Terrain correction systems |
GB2451807B (en) | 2007-08-02 | 2012-01-18 | Arkex Ltd | Geophysical data processing systems |
CA2761863A1 (en) * | 2009-05-15 | 2010-11-18 | Arkex Limited | Geophysical data processing systems |
GB2471682B (en) | 2009-07-07 | 2014-01-01 | Arkex Ltd | Potential field data survey |
US20110029291A1 (en) * | 2009-07-31 | 2011-02-03 | Xiaowei Weng | Method for fracture surface extraction from microseismic events cloud |
WO2011098821A2 (en) | 2010-02-12 | 2011-08-18 | Arkex Limited | Geophysical data processing systems |
GB201008993D0 (en) | 2010-05-28 | 2010-07-14 | Arkex Ltd | Processing geophysical data |
GB2481643A (en) | 2010-07-02 | 2012-01-04 | Arkex Ltd | Gravity survey data processing |
GB2482505A (en) | 2010-08-04 | 2012-02-08 | Arkex Ltd | Using potential field data and seismic data to generate a representation of the underlying geology of an underground formation |
GB2489230B (en) | 2011-03-21 | 2014-06-18 | Arkex Ltd | Gravity gradiometer survey techniques |
GB2503371B (en) * | 2011-03-25 | 2016-11-02 | Baker Hughes Inc | Use of frequency standards for gravitational surveys |
US9964653B2 (en) | 2011-12-21 | 2018-05-08 | Technoimaging, Llc | Method of terrain correction for potential field geophysical survey data |
CA2882128C (en) * | 2012-08-15 | 2018-02-13 | Bell Geospace Inc. | Directional filter for processing full tensor gradiometer data |
CN103399350B (zh) * | 2013-07-29 | 2016-02-24 | 中国人民解放军国防科学技术大学 | 一种基于积分迭代算法的航空重力向下延拓方法 |
KR102126507B1 (ko) * | 2013-12-09 | 2020-06-24 | 삼성전자주식회사 | 센서 데이터 스트림을 처리하는 단말기, 시스템 및 방법 |
CN106772621B (zh) * | 2017-01-24 | 2019-11-19 | 山东大学 | 一种近全方位电阻率隧道超前地质预报方法 |
WO2021021097A1 (en) * | 2019-07-29 | 2021-02-04 | Landmark Graphics Corporation | Gas saturation distribution monitoring in hydrocarbon reservoir |
CN111458761B (zh) * | 2020-04-16 | 2022-09-09 | 自然资源部第二海洋研究所 | 海上重力比对场建设方法 |
CN113606752B (zh) * | 2021-06-29 | 2023-03-03 | 宁波德业日用电器科技有限公司 | 避免跳变的除湿机湿度显示方法 |
US11906694B2 (en) * | 2022-03-24 | 2024-02-20 | Saudi Arabian Oil Company | Method to estimate the depth of the weathering layer using gravity response |
CN115421212B (zh) * | 2022-06-21 | 2024-04-02 | 中国科学技术大学 | 重磁图件补全与去噪同步处理方法及系统 |
Family Cites Families (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
RU2091820C1 (ru) | 1994-02-15 | 1997-09-27 | Научная станция Института высоких температур РАН | Геофизическая система сбора и обработки информации |
US5673191A (en) | 1995-04-10 | 1997-09-30 | Atlantic Richfield Company | Method and apparatus for identifying geological structures using wavelet analysis of potential fields |
US5886255A (en) * | 1997-10-14 | 1999-03-23 | Western Atlas International, Inc. | Method and apparatus for monitoring mineral production |
US6675097B2 (en) * | 1999-04-02 | 2004-01-06 | Conoco Inc. | Nonlinear constrained inversion method to determine base of salt interface from gravity and gravity tensor data |
US6993433B2 (en) * | 1999-04-02 | 2006-01-31 | Conocophillips Company | Modeling gravity and tensor gravity data using poisson's equation for airborne, surface and borehole applications |
GB0121719D0 (en) * | 2001-09-07 | 2001-10-31 | Univ Edinburgh | Method for detection fo subsurface resistivity contrasts |
DE102004005676A1 (de) | 2004-02-05 | 2005-08-25 | Giesecke & Devrient Gmbh | Datenträger mit plattformunabhängigem Anwendungs-Programmcode |
US7065449B2 (en) * | 2004-03-05 | 2006-06-20 | Bell Geospace, Inc. | Method and system for evaluating geophysical survey data |
US7542850B2 (en) * | 2004-06-24 | 2009-06-02 | Bell Geospace, Inc. | Method and system for synchronizing geophysical survey data |
CN1797035A (zh) * | 2004-12-29 | 2006-07-05 | 中国石油天然气集团公司 | 三维复杂地质体重构的技术 |
AU2006273791B2 (en) * | 2005-07-27 | 2011-10-06 | Arkex Limited | Gravity survey data processing |
GB2435523B (en) | 2006-01-25 | 2010-06-23 | Arkex Ltd | Terrain correction systems |
-
2007
- 2007-01-30 GB GB0701725A patent/GB2446174B/en not_active Expired - Fee Related
-
2008
- 2008-01-30 CN CN2008800021674A patent/CN101600975B/zh not_active Expired - Fee Related
- 2008-01-30 CA CA002674535A patent/CA2674535A1/en not_active Abandoned
- 2008-01-30 WO PCT/GB2008/050059 patent/WO2008093139A2/en active Application Filing
- 2008-01-30 AU AU2008211708A patent/AU2008211708B2/en not_active Ceased
- 2008-01-30 ZA ZA200902158A patent/ZA200902158B/xx unknown
- 2008-01-30 US US12/311,750 patent/US8332184B2/en not_active Expired - Fee Related
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116794735A (zh) * | 2023-06-02 | 2023-09-22 | 中国自然资源航空物探遥感中心 | 航空磁矢量梯度数据等效源多分量联合去噪方法及装置 |
CN116794735B (zh) * | 2023-06-02 | 2024-03-19 | 中国自然资源航空物探遥感中心 | 航空磁矢量梯度数据等效源多分量联合去噪方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
US8332184B2 (en) | 2012-12-11 |
US20090287464A1 (en) | 2009-11-19 |
AU2008211708B2 (en) | 2012-12-13 |
RU2009132490A (ru) | 2011-03-10 |
AU2008211708A1 (en) | 2008-08-07 |
ZA200902158B (en) | 2010-07-28 |
WO2008093139A2 (en) | 2008-08-07 |
GB2446174A (en) | 2008-08-06 |
GB2446174B (en) | 2011-07-13 |
WO2008093139A3 (en) | 2009-04-23 |
CN101600975B (zh) | 2013-01-23 |
GB0701725D0 (en) | 2007-03-07 |
CA2674535A1 (en) | 2008-08-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101600975B (zh) | 重力勘测数据处理 | |
CA2920499C (en) | Stratigraphic function | |
CN101999086B (zh) | 用于确定地震数据质量的方法 | |
CA2576586C (en) | Method and system for processing geophysical survey data | |
EP3413092B1 (en) | Method for evaluating a geophysical survey acquisition geometry over a region of interest, related process, system and computer program product | |
Shamsipour et al. | 3D stochastic inversion of gravity data using cokriging and cosimulation | |
US7345951B2 (en) | Method of determining specular information after prestack seismic imaging | |
US20180113235A1 (en) | Geologic Stratigraphy Via Implicit and Jump Functions | |
EP3149517A2 (en) | Properties link for simultaneous joint inversion | |
CN101772716A (zh) | 地球物理数据处理系统 | |
Rossi et al. | Integrating geological prior information into the inverse gravimetric problem: the Bayesian approach | |
GB2465715A (en) | Gravity survey data processing | |
Berdyshev et al. | Mapping Problems of Geophysical Fields in Ocean and Extremum Problems of Underwater Objects Navigation | |
Medagoda et al. | Water column current profile aided localisation for autonomous underwater vehicles | |
CN103140777A (zh) | 用于处理地球物理数据的系统和方法 | |
Rossi | Bayesian gravity inversion by Monte Carlo methods | |
Ashena et al. | A Novel 2.5 D Deep Network Inversion of Gravity Anomalies to Estimate Basement Topography | |
US20190146108A1 (en) | System and method for assessing the presence of hydrocarbons in a subterranean reservoir based on seismic data | |
Bedada | Absolute geopotential height system for Ethiopia | |
Garcia et al. | Interfacing an ensemble Data Assimilation system with a 3D nonhydrostatic Coastal Ocean Model, an OSSE experiment | |
Stonebridge | Diagonal approximations to the observation error covariance matrix in sea ice thickness data assimilation | |
RU2486549C2 (ru) | Обработка данных гравиметрической съемки | |
Sambolian et al. | Revisiting the hypocenter-velocity problem through a slope tomography inspiration | |
ROSSI | Bayesian interpretation of gravity data with geological prior information | |
Barzaghi et al. | Exploitation of GOCE data for a local estimate of gravity field and geoid in the Piemonte area (northern Italy) |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20130123 Termination date: 20150130 |
|
EXPY | Termination of patent right or utility model |