专利名称：用于地球物理勘测的地形校正方法
技术领域：
本发明涉及用于针对地球物理勘测的地形校正的方法、设备及计 算机程序代码，具体涉及用于针对势场勘测的地形校正的方法、设备 及计算机程序代码。
背景技术：
为了方便起见，在说明书中，我们将主要涉及航空勘测，更具体 地涉及重力勘测。然而，我们所描述的技术并不局限于这些勘测类型， 并且可应用于基于空中、陆地或是海上的重力勘测，或者， 一般地说， 可应用于包括但不局限于诸如大地电磁勘测、电磁勘测等的磁场勘测 的其它势场勘测。通过测量势场数据来执行势场勘测，对于重力勘测， 该势场数据可以包括重力计数据(测量重力场)或者重力梯度计数据 (测量重力场梯度)。在其它勘测类型中，势场数据包括矢量磁力计数 据、真磁梯度计数据、或者本领域技术人员公知的其它类型的数据。 地球物理势场勘测的共同目的是寻找潜在地指示有价值矿床的特征。
在重力勘测中，附近的质量提供高和低(有关空间的)频率数据， 而更深质量的作用主要仅集中于较低频率。当查看下层异常时，中间 质量具有主导作用，并且为了提供对深层特征的精确表征，需要对表 面特征的良好表征，以便能够特别减去较高频率(在功率谱中占优势)。
例如，波长义的信号按照邵Hb)随着高度z下降，其中h2w;i，由此 可以估计来自深度为100米的质量的、波长为200米的信号分量在地 球表面下降至其初始值的大约1/20 (随着高度的增加而逐渐衰减)， 然而，可以看出波长越长则衰减得越少。通常，在给定目标大小和深 度的情况下，根据与所期望的特征相对应的波长刻度广泛地选择勘测 的大小和位置。
按照惯例，在网格图形上进行诸如重力勘测之类的航空势场勘测
(背景技术可以在US2005/0017721和EP0589554A中找到)。通过在 下层地形上覆盖的二维表面上的平行线(飞行路径)的正交集合来定 义网格。覆盖表面满足最小高度限制(由允许航空器距离地面飞行的 最小距离来定义)、以及针对航空器攀升/下降的最大速率(典型约为 百分之三)的限制。然而，在丘陵或山脉地形，航空器在其上飞行的 表面从下层山谷的底部至山脉/测区的顶部会改变二或三千米，因而需 要另 一 种方式。在申请人的共同未决的PCT专利申请"Gravity Survey Processing", PCT/GB2006/050211 (通过引用全部合并于此)中描述 了一些用于航空势场勘测的改进技术。这些有利于势场数据的收集， 尤其是接近地面的重力数据的收集。
在已经收集到势场数据之后，通常在分析数据之前实施地形校正。 这包括补偿表面高度变化，或者更具体为在勘测区域表面之上的航空 器高度的变化。可以采用许多技术来获得表面剖面图，包括但不局限
于下列操作中的一个或更多个:购买本地数字地形高程数据、获取GPS
(全球定位系统)数据、以及诸如Lidar (激光雷达)和SAR (合成 孔径雷达)之类的航空技术。为了给出对校正量级的指示，地形高度 变化的典型特征是数以百计的Eotvos，而下层地质构造的典型特征是 1-10个Eotvos。
尽管如我们之前已经描述的飞行和数据处理技术有所改进，仍然 存在对所勘测区域以下的地形作用更好地进行校正的需要。
发明内容
根据本发明，因此这里提出了一种针对地球物理勘测的地形校正
的方法，该方法包括捕获要勘测区域的多光谱或超光谱图像；使用
所述捕获的图像(即多光谱或超光谱图像)确定所述勘测区域的表面
地质组成；使用所述确定的表面地质组成确定所述地球物理勘测的地 形校正数据；以及使用所述地形校正数据执行针对所述地球物理勘测 的所述地形校正。
优选地，该地形校正包括定义勘测区域上的表面密度或多孔性的 变化的密度数据。在实施例中，通过采用超光谱成像来计算地形密度 (以及高程)，可以充分地衰减地形的不期望的"虚"图像。这有助于 有效地减少平面层上的地形。对于相对多孔的岩石，密度可以取决于 最近或平均降雨量，因而多孔性是有用的。
绘制表面地质组成可以包括将组成分类成多个类型之一(或更 多)；这可以基于表面岩石(这里通常包括粘土)的直射影像，例如在 植被很少或几乎没有植被的地方，或者分类可以基于表面植被的变化。 将理解，不需要精确的密度来提供某些益处，并且表面密度的近似估 计甚至可以是有用的。绘制的密度或多孔性可以包括相对值，而不是 绝对值。
在实施例中，可以釆用一种或多种不同技术来根据勘测区域的多 光谱或超光谱图像确定表面地质组成。这可以包括(但不局限于)对 所捕获的图像进行主分量分析；将成像的光谱与所存储的光谱数据的 库进行匹配；对所捕获的光谱数据进行匹配滤波；以及对所捕获的图 像数据进行基于人工神经网络的分析。可选地，可以使用线性或非线 性插法，像素或区域中的两个或更多个确定的地质组成之间进行内插。 例如，单个像素可以与两个相似的、但不同的岩石类型紧密匹配，和/ 或附近像素可以匹配不同的相关岩石类型。在这样的情况下，可以采 用岩石类型和/或光谱的线性组合来确定与所识别的岩石类型相关联 的那些之间的密度值。
优选地，地球物理勘测还包括重力或重力梯度勘测，尽管通常可 以利用包括根据周围场(例如大地电磁)或采用人工电磁场源的任意 势场勘测来实施上述技术。当利用重力勘测来实施地球物理勘测时， 优选地将地形校正密度数据与勘测区域的表面高程数据结合起来，该 表面高程数据可以通过任意传统技术获得或购买。
可以确定地形校正数据，并将其提供用于地球物理勘测的目的， 或者可以结合势场测量数据，直接采用地形校正数据来解决"反演问 题"，以确定勘测区域的估计的下层地质组成。
通常，这种反演问题的解决方案是受约束的，因此例如来自物理 表面勘测和/或和钻孔的任意可用先验数据都是有用的。
可以应用诸如反演问题之类的的解决解决方案的传统技术，除了
将固定密度用于勘测区域的地形校正，可以采用测量的(估计的)表 面密度。由于采用数值技术，因而可以直接修改现有的计算机软件， 以达到上述目的。同样地，表面地形的每个元素都有不同高度，也可 以分配不同的密度；对于表面密度未知或无法充分精确地估计表面密 度的部分勘测区域(如果有的话)，可以采用缺省密度值，例如2.67g/cc 的典型的上层大陆地壳密度。
在优选实施例中，多光谱或超光谱捕获图像具有至少3个波段， 更有优选地为至少5、 10、 50或100个波段。例如在优选实施例中， 采用了 393个特定波段。优选地，至少一个波段扩展到紫外波段，特 别是具有小于或等于200nm的波长。在优选实施例中，针对表面地质 组成的确定所采用的多光谱或超光谱数据主要来自于紫外(小于 400nm)的光谱区域。
本发明还提供了一种计算机程序代码，特别是在载体上提供一种 计算机程序代码，用于在运行时实现上述方法的实施例。该载体可以 包括诸如DVD或CD-rom之类的盘、程序存储器、或光或电信号载体。 该计算机程序代码可以使用任何传统的计算机语言来编写，包括诸如 C或Fortran之类的高级语言以及诸如汇编码之类的的低级语言。该代 码可以包括源、对象、或可执行代码，并且可以例如通过网络分布。
本发明还提供一种使用上述代码编写的传统计算机系统，其包括 用于确定地形校正数据的工作存储器以及处理器。优选地，该系统连 接到用于使超光谱图像与岩石类型相匹配的系统，并且包括用于使每 一岩石类型与密度和/或多孔性相关联的数据存储器。优选地，还提供 了一种终端，用于交互确定下层地质结构的三维模型，该终端包括来 自其它源的、诸如钻孔数据等的数据的输入(在可用的情况下)。
现在将参考附图，仅作为示例进一步描述本发明的这些和其它方
面，在附图中
图1示出了体现本发明的方面的系统的框图； 图2示出了用于实现根据本发明方法的实施例的过程的框图。
具体实施例方式
参照图1，航空器ioo装备了多个地球物理勘测仪器，其中地球
物理勘测仪器包括一个或更多个势场测量仪器102、至少一个超光谱 成像器(HSI) 104以及其它仪器106。例如，势场测量仪器可以包括 矢量重力计、重力梯度计、磁力计、磁梯度计等。其它仪器可以包括 GPS (全球定位系统)仪器(优选地为DGPS (差分GPS)设备)、高 度计、高度测量设备、压力测量设备等。可选地，可以包括电磁测量 (EM)仪器用于势场测量，例如时域电磁系统(TDEM)。
在一个实施例中，利用LIDAR (激光雷达)和IMU (惯性测量单 元)的组合，将航空器装配为产生精确的DEM (数字高程模型)。结 合DGPS，以便针对航空器运动校正LIDAR数据。DEM和DGPS数 据还可用于校正针对地形的测量的势场数据。同样地，还可以使用航 空器的加速度、飞行姿势、角速率和角加速度数据来校正势场仪器的 输出数据。本领域的技术人员将理解，实质上可以采用任意机上或远 程传感器来给航空器和/或势场仪器提供位置和运动信息。
接下来转向超光谱成像器104，概括地说，本领域的技术人员将 理解，该设备包括用于在宽波长范围上的多个窄波段上进行测量的一 个或更多个成像光谱仪。这里所描述的技术不局限于任何特定类型的 超光谱成像器，但是部分由于红外信号在白天被淹没，该设备优选地 具有紫外波段的多个波段。该成像器可以是无源的或有源的；在后者 的情况中可以采用类似于LIDAR的技术。
在一个示例性实施例中，在393个波段中釆用超光谱成像器来捕 获图像，在实施例中，393个波段覆盖从小于150nm或100nm到500nm 或更长的波长区域(尽管在特定条件下IR会被淹没，但在其它条件下， 光谱范围的IR部分是可用的)。，
技术人员将理解，尽管通常很方便在航空器IOO上安装超光谱成 像器104，在所描述的技术的备选实施例中，可以采用通过卫星捕获 的超光谱图像数据(尽管分辨率通常比较低，成本较高，并且光谱覆 盖也不是最佳的)。
将超光谱成像数据和(如果同时实施)势场勘测数据提供给数据
处理系统108，该数据处理系统108包括处理器108a、用于地形校正 的程序代码108b、以及用于根据势场数据确定勘测区域的地质构造的 三维模型的优选的反演问题解决程序代码108c。数据处理系统108可 选地通过网络耦合至包括光谱匹配/分类码的HSI图像库110。该库接 收来自数据处理系统108的HSI数据，并返回由地形校正程序代码所 使用的岩石密度数据。提供用户终端112，用于输入关于可从己知数 据获得的反演问题的任意约束，例如该己经数据源自表面而非机载的 勘测。可选地，其它数据源还可以提供类似的输入，通常可以提供的 关于下层地质的已知信息越多，则产生的模型越好。在航空器100没 有收集地形高程数据的情况下，可以分别提供DEM (数字高程模型) 数据114。数据处理系统108提供包括地形校正数据、以及可选地包 括三维地质模型数据的输出116。可以将该数据存储在数据库118中 和/或通过网络输出或在可移动存储介质120上提供该数据。
现在参照图2，图2示出了图1中的系统的过程的流程图。最初 从飞行勘测200 (在其它实施例中，可以采用陆地或海上勘测)中获 得数据，从而提供了用于确定岩石类型并用于绘制岩石密度和/或多孔 性202所使用的超光谱成像数据。如上所述，这包括对与所捕获的图 像中的每个像素相关联的光谱与所存储的HIS数据库进行匹配，以便 对下层表面岩石进行分类并因此获得密度/多孔性数据。上述提到的这 一点可以解决匹配与获取图像的每一像素关联的光谱图和存储在HSI 库中的数据。可以通过采用来自不同类型的下层岩石的典型反射系数 对HSI数据进行分类，或者可以在物理上本地检查勘测区域内的一个 或更多个位置的地质，然后可以相对于这些位置来绘制相对变化；备 选地，可以采用这两种技术的综合。例如，云母包括三种不同密度的 碳酸盐单元，并且可以采用本地采样来最好地匹配在超光谱图像数据 中看到的CO特征的吸收带。可选地，用于捕获超光谱图像数据的波 段可以根据所期望的岩石类型而变化，例如含有铁的砂岩最好在红外 波段勘测。可以采用主分量分析来将像素分类成岩石类型，例如石灰 石、粉砂石灰、砂石等。在这些类型的岩石具有相似的密度范围的情 况下，可以将这些类型可选地结合起来；另外或备选地，在不同、但 邻近的像素定义了不同类型、但处于相似密度范围内的岩石的情况下， 不需要精确的分类来确定勘测区域图像的相关部分的近似密度。然而， 概括地说，该技术要基于超光谱图像数据将勘测区域图像划分为岩石 类型的类别。可以采用类似的技术，从表面土壤中生长的植被的超光 谱图像中确定下层岩石类型或至少岩石类型的类别。
商业上可用的超光谱图像处理软件包括从RSI公司获得的ENVI 软件；例如，还可以通过分析成像和Colorado, USA的Geophysicals LLC提供HIS数据处理服务。
超光谱图像分析技术的更多详细资料可以在美国科罗拉多州RSI 公司的ENVI(可视图像环境)软件指南 (www.rsinc.com/envi/tutorial.pdf) 中找至il， 尤其是12至Ul7章。
一旦在步骤202已经绘制了出岩石密度/多孔性，该过程确定地形 校正数据，其实质上包括结合数字高程模型的绘制的岩石密度，该数 字高程模型源自勘测或其它源206。例如，可以使用LIDAR (激光雷 达)和IMU (惯性测量单元)的组合结合DGPS (差分全球定位系统) 来产生精确的DEM(数字高程模型)，以针对航空器运动来校正LIDAR 数据。DEM和DGPS数据也可以用于针对地形来校正测量的势场数 据。同样地，航空器加速度、飞行姿势、角速率和角加速度数据也可 以用于校正势场仪器的输出数据。任何机上或远程传感器可以用于提 供针对航空器和/或势场仪器的位置和运动信息。
因此，该地形校正数据包括在勘测区域上的X、 Y位置处定义高 度(高程)和密度(或多孔性)的一组数据。这用于校正势场数据， 目的是有效地将地形减少到平坦、均匀的密度层。在执行该校正之前， 通常对势场数据进行预处理208，例如在本技术的某些优选实施例中 通过使用我们早期PCT专利申请(同上)中所描述的等效源方法的网 格化。因此，这对于在本阶段概述这些方法是有帮助的。
势场数据包括但不局限于重力计数据、重力梯度计数据、矢量磁 力计数据和真磁梯度计数据。可以通过标量来推导出势场的元素和表 征。
对于重力，相关势能是重力标量势O)(r)，定义为
其中，r、 P(o、 G分别是重力场测量的位置、位置,处的质量密度、 以及重力常数。重力(所体验到的重力场)是标量势的空间导数。重
力是方向已知(重力向下作用)的矢: 标系统，通过三个分量来表示重力<formula>formula see original document page 11</formula> 岭)3岭)滩(。、
对于任何所选的Cartesian坐
对于这三个分量，每个分量沿三个方向中的每个方向而改变，因 而产生了九个量，形成了重力梯度张量
<formula>formula see original document page 11</formula>
势场的基本方程和关系遵循对标量势函数、其导数、其傅立叶变 换和其它数学量的属性的分析。
根据格林定理，在己知闭合表面上的标量势的空间导数(包括标 量势本身)的情况下，在该表面所包围的体积内的所有点上已知该空 间导数值。其推论是， 一旦该量对于所有点都是己知的，则通过微分 和积分，可以导出包括标量势自身的标量势的所有其他导数。因而当 在包围了该体积的表面上仅已知导数之一时，则有效地知道了该体积 内的所有点上的标量势和所有导数。这指示了，标量势的任何导数的 任何分量的全部测量允许计算标量势的任何导数的任何其他分量。据 此得出(至少在理论上)测量哪个量并不重要，装置的选择简单地归 结于哪个装置以最大信噪比测量了期望量。
以上重力标量势的微分最终产生
该方程在不重要的区域内变为拉普拉斯方程，得到重力的重要基本关
系
<formula>formula see original document page 12</formula>
谐波函数满足拉普拉斯方程，其具有可以在对从势场勘测中采集 到的数据进行分析的过程中使用的许多属性。
可以使用多种技术来分析和处理数据，这些技术与从勘测中采集 到的数据一起起作用，作为起点，但是之后改变数据和/或数据格式， 因而与所测量的量关联的值全部出现在位于水平的、固定海拔的分析 平面上的规则2-D网格上，将所使用的过程称为"分级(levelling)" 和"网格化(gridding)"。
在网格化中，概括地说，将所勘测区域划分为其侧面优选地与进 行勘测的主要方向对齐的矩形小区，然后利用"等效于"测量数据的 数据(网格化数据)来代替实际测量数据，该测量数据现在分配有每 个小区中部的点处的值。每个小区的尺寸可以基于沿两个正交方向飞 行的线路的平均间隔相关。 一旦数据采取这种网格化格式，则在数学 上容易处理得多。例如，可以通过统计或其它方法将数据示为一组数 值，并对其进行处理，以给出对水平分析平面上的势场的最佳估计。
可以将数据变为2-D傅立叶序列，其中每行数据具有2n个数据点，
以有利于使用傅立叶变换；还沿正交方向采集数据。对于重力梯度， 可以将2-D傅立叶序列表示为以下形式的2-D空间正弦波之和
<formula>formula see original document page 12</formula>
其中波数、，、与勘测的尺寸相关，x、 y方向上的4和、分别由以下
表示 ^ = — , = 一 Ac 丄，
然而这种扩展仅在恒定高度时有效，除非L是Z的函数。因此， 我们先前己经描述了 一种改进的、等效源类型方法 (PCT/GB2006/050211，通过引用合并)在示例性等效源类型方法中， 将勘测地区的表面分成小块，典型地为一侧50m的量级，这些被称为 小板或质量元素。很容易前向计算出每个小板的重力(例如，见 R.J.Blakely的"Potential Theory in Gravity and Magnetic Applications", Cambridge University Press, 1995)，调整其质量直至获得对测量数据的 最佳整体拟合。通过标准最小平方拟合过程来做出该质量确定。通过 将真实测量位置处的数据与在相同的真实测量位置处的所提出的等效 源产生的重力场匹配来获得该拟合。该过程在数学上是精确的，不会 将任何人为调整引入数据中，以便其以与水平矩形勘测一致。
一旦获得了拟合，便认为该拟合是最初的数据集合。所有用于确 定地质结构的后续分析比较任何给定地质结构将会产生的重力场与由 等效源所产生的重力场之差，并将该差最小化。这种技术的一个显著 优点在于，最佳拟合来自于质量分布(尽管是合成的)，因而最佳拟合 情况将自动满足拉普拉斯方程。这不同于上述的传统方法，上述传统 方法产生数值最佳拟合，但是并不施加该数据必须满足拉普拉斯方程 的附加限制，即数据来自真实的质量分布的限制。
等效源方法不必使用与地形共形的表面，可以使用可以覆盖位于 恒定高度、在地球真实表面之上或之下的任何表面、可以穿过地球真 实表面等的源。在各个小板的质量中，符合该地形的表面的选择可能 产生较小的变化，但是原理上整个结果并不会受到表面的任何合理选 择的显著影响。
通过网格或上述等效源方法所模拟的势场数据提供了对计算机 程序代码210的输入，用于确定造成所测量的(和模拟的)势场的三 维下层地质结构，并考虑步骤204中确定的地形校正以及任意可用先 验数据212。在文献(例如，参见Blakely (同上))描述了适当的技 术，可以另外地或备选地采用可通过商业手段获得的程序代码，例如，可以从基于美国Utah大学的电磁模拟和反演研究组(CEMI)获得的 GRMAG3D代码。
毫无疑问，本领域的技术人员将想到许多其它有效的备选。将被 理解的是，本发明不局限于所描述的实施例，并且包括对于本领域技 术人员而言显而易见的、落入所附权利要求的精神和范围内的修改。
权利要求
1.一种针对地球物理勘测的地形校正的方法，所述方法包括捕获要勘测区域的多光谱或超光谱图像；使用所述捕获的图像，多光谱或超光谱图像，确定所述勘测区域的表面地质组成；使用所述确定的表面地质组成，确定针对所述地球物理勘测的地形校正数据；以及使用所述地形校正数据来执行针对所述地球物理勘测的所述地形校正。
2. 根据权利要求l所述的方法，其中，所述地形校正包括定义所 述勘测区域上的表面密度或多孔性的变化的密度数据。
3. 根据权利要求2所述的方法，其中，所述密度数据确定包括.-针对所述捕获图像中的一个或更多个像素，在两个或更多个表面地质 组成之间进行内插。
4. 根据权利要求2或3所述的方法，其中，所述确定所述地形校正数据包括将所述密度数据与针对所述勘测区域的表面高程数据结合起 来。
5. 根据前述任一权利要求所述的方法，其中，所述表面地质组成 确定响应于所述勘测区域中的表面植被的变化。
6. 根据前述任一权利要求所述的方法，其中，所述表面地质组成 确定包括将所述勘测区域的各部分分类成多个表面地质组成类别之
7. 根据前述任一权利要求所述的方法，其中，所述表面地质组成 确定包括下列操作中的一个或更多个对所述捕获的图像进行主分量 分析、使用存储的光谱数据对所述捕获的图像进行匹配操作、对所述 捕获的图像进行匹配滤波、以及对所述捕获的图像进行基于人工神经 网络的分析。
8. 根据前述任一权利要求所述的方法，其中，所述地球物理勘测 包括重力或重力梯度勘测。
9. 根据前述任一权利要求所述的方法，还包括执行所述地形校正 的地球物理勘测，以确定所述勘测区域的估计的下层地质组成。
10. 根据前述任一权利要求所述的方法，其中，所述多光谱或超 光谱捕获图像包括至少3个波段，优选地为至少5、 10、 50或100个波段。
11. 根据前述任一权利要求所述的方法，其中，所述多光谱或超 光谱捕获图像包括多个波长等于或小于400nm的波段。
12. —种承载计算机程序代码的载体，所述计算机程序代码在运 行时实现前述任一权利要求所述的方法。
13. —种计算机系统，被编程为实现前述任一权利要求所述的方法。
全文摘要
本发明总体涉及用于针对地球物理勘测的地形校正的方法、设备和计算程序代码，尤其是用于针对势场勘测的地形校正的方法、设备和计算程序代码。针对地球物理勘测的地形校正方法包括捕获要勘测区域的多光谱或超光谱图像；使用所述捕获的图像，即多光谱或超光谱图像，确定所述勘测区域的表面地质组成；使用所述确定的表面地质组成确定针对所述地球物理勘测的地形校正数据；以及使用所述地形校正数据执行针对所述地球物理勘测的所述地形校正。
文档编号G01V3/38GK101371165SQ200780002995
公开日2009年2月18日 申请日期2007年1月24日 优先权日2006年1月25日
发明者马克·戴维斯 申请人:阿克斯有限责任公司