CN102636819A - 重力测量数据处理 - Google Patents

Abstract

本发明涉及处理来自诸如重力测量之类的航空测量的势场测量数据的改进技术，并涉及通过这种改进的数据处理技术实现数据获取的改进技术。描述了一种处理来自地表的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的方法，所述方法包括：输入所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；定义所述势场测量和所述位置之间的多个关系，每个所述关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；以及针对所述多个关系确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据。

Description

重力测量数据处理

本申请是2006年7月17日提交的、国际申请号“PCT/GB2006/050211”、国家申请号“200680035762.9”、发明名称为“重力测量数据处理”的发明的分案申请。

技术领域

本发明涉及一种处理来自诸如重力测量(gravity survey)之类的航空测量的势场测量数据的改进技术，并涉及通过这种改进的数据处理技术来实现数据获取的改进技术。

背景技术

传统地，在网格图形上进行诸如重力测量之类的航空势场测量。通过在下层地形上覆盖的二维表面上的平行线(飞行路径)的正交集合来定义网格。覆盖表面满足最小高度限制(由允许航空器距离地面飞行的最小距离来定义)、以及针对航空器攀升/下降的最大速率(典型约为百分之三)的限制。该方式满足平坦地形的需要，但是针对丘陵或山脉地形，航空器在其上飞行的表面从下层山谷的底部至山脉/测区的顶部会改变二或三千米，因而需要另一种方式。

能够从靠近地面(即，低高度)的地方采集势场数据(尤其是重力数据)是非常有用的。在重力测量中，附近的质量(mass)提供了高和低(空间)频率数据，而较深质量的影响看似仅主要在较低频率处。当查看下层的异常布局时，居间质量具有主要作用，为了提供深层特征的准确表示，期望有表面特征的良好表示，从而能够减去尤其是(主导功率谱的)较高频率。例如，波长为λ的信号按照exp(-kz)随高度z衰减，其中k 2π/λ，从中可以估计来自深度100米的质量的波长为200米的信号分量在地球表面处降至其初始值的大约1/20(并且随不断增高的高度而进一步衰减)，而可以看出波长越长，衰减越少。通常根据与在给定目标的尺寸和深度下预期的标记相对应的波长尺寸来宽泛地选择测量的尺寸和位置。

从上述讨论中可以理解，通常期望能够在较低的高度执行飞行测量，但是实际上航空器的局限性和网格飞行规划会导致较大的限制。由于处理重力测量数据的传统技术取决于恒定高度假设，因而网格飞行规划是必要的。广义上，该假设是对于深的源，航空器的高度近似恒定，由于该假设所导致的误差按照推测仅仅是小修正。此外，传统的重力测量数据处理技术取决于均匀间隔的数据点，通常是2的幂，从而可以使用(快速)傅立叶变换，这定义了对平行飞行路径的正交集合的要求。由于现有的技术假设(例如，在两条路径交叉的情况下，它们交叉于相同高度)，需要该飞行路径在共同表面上。在所测量的区域并不是精确的矩形(例如，由于局部地形)时，现有飞行路径出现另一问题。在这种情况下，为了能够应用传统技术，例如，通过内插或外插产生矩形区域上的数据点的正则集合，来“填充”数据点。然后，通过现在填充后的矩形区域的最大x和y(长和宽)维度来确定用于傅立叶分析的波长(或更具体地，波数)。

针对传统技术的上述缺点，需要改进的势场测量数据处理技术和测量飞行图形。

发明内容

根据本发明的第一方面，提供了一种处理来自地表的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的方法，所述方法包括：输入所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；确定所述势场测量和所述位置之间的多个关系，每个所述关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；以及针对所述多个关系确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据。

势场数据可以包括重力计数据(测量重力场)、重力梯度计数据(测量磁场梯度)、矢量磁力计数据、真磁梯度计数据、或其他类似数据。场地图绘制参数(可以表征势场特性的系数或参数)可以用于产生类似的势场数据，在实施例中，广义上，通过选择定义了诸如重力场或重力梯度场之类的最接近于所测量势场数据的势场数据的参数值来确定参数。例如，这可以包括将均方差最小化或另一类似的优化。

在实施例中，所述关系包括例如可以通过矩阵表示的联立方程组，以及所述确定包括求解这些方程。优选地，该方程是过约束的，并且该求解包括通过诸如最小平方优化之类的降噪过程进行求解。优选地，所述关系具有实(重力)势场的形式，从而通过该优化过程有效地丢弃通常不具有该形式的噪声。因而，优选地，(三维)位置函数满足拉普拉斯方程。在所述方法的一些代表性示例中，可以有10K数量级个参数和100K数量级个位置。

地图数据可以只包括场地图绘制参数(势场系数或参数)集，因为该数据可以用于依据所期望的高度、覆盖区域和/或要绘制地图的表面来产生大量势场数据(例如重力场或重力梯度场数据)。然而，所述方法还可以包括根据产生地图的场地图绘制参数集来确定相关场(典型在一个表面上)，其根据场地图绘制参数集合来执行势场数据的前向计算。如之后将再次提及的，所述方法的实施例(在通过致力于将所计算和所测量的值之差最小化或进行优化的过程来确定场地图绘制参数时至少暗含地)产生与所测量的势场数据相同类型的势场数据。通过格林定理的应用，可以看出无论势场数据包括重力计数据还是重力梯度计数据，这都是有效的。此外，知道一种类型的势场数据(如重力场或重力梯度场)允许通过积分/微分来确定另外一种类型的势场数据、以及标量势场。

在方法实施例中，采用势场数据的等效源表示。在该实施例中，发现针对多个表面质量元素的表面密度或(质量)值，同时产生最接近于所测量的重力(或其他)场的重力(或其他)场。通常可以选择该表面(但是不需要是平面的)接近于所测量区域中的地球表面。一旦发现了这些质量元素的值，则可以采用简单的前向计算来预测标量势场或重力场或重力梯度场，以例如产生所测量区域的任何期望地图。典型地，做出前向计算以确定平坦地图绘制平面上的场的一个或多个分量。

概括地，在该方法的实施例中，所测量的场与根据等效源质量所确定的场(无论测量/计算的是重力还是重力梯度)之差的平方是质量误差的函数。因而通过该函数单独地相对于每个表面元素的质量的部分微分，形成了联立方程组，其中唯一的变量是每个表面质量元素的质量(假设已知相对于表面质量元素的测量位置r-r’)。方便地，可以通过矩阵表示该联立方程组，并针对表面质量元素值进行求解。由于测量数量通常远大于联立方程的数量(即质量元素的数量)，例如，至少大五或十倍，因而该方程组是过约束的。这一点是有用的，因为如果不是这种情况，则表面质量元素的值将拟合噪声；作为替代，使用优化过程(例如最小平方拟合)来确定表面质量元素的值，从而减小噪声。或者此外，由于使用了物理模型重力场，因而方程遵守拉普拉斯方程，从而趋向于拟合的是重力场而非噪声。

在所述方法的第二实施例中，采用了修改后的傅立叶基本函数技术。在该方法中，使用二维傅立叶展开来表示所测量的数据，但是利用高度相关系数

乘以每个二维分量该表示仍完全通用并遵守拉布拉斯方程(当进行微分时，该高度相关系数和2D分量的乘数抵消)。可以取决于该表示的期望精度(例如，由相关数确定)在一点上截短该展开。可以采用其他展开，例如球谐函数，其对于较大区域是有用的，尤其在地表的曲率很大的情况下。

广义上，按照以下方式确定展开系数。在位置r，已知所测量的势场(重力或重力梯度)，但是对于特定值对k_m和k_n，

值也是已知的。因而在傅立叶展开(例如，参照后面的方程(3))，有效地单独看待每个数据点，如上述第一实施例，产生可以用于对傅立叶展开系数进行求解的联立方程组。然后可以采用这些系数来确定(通过前向计算)任何二维表面上的重力或重力梯度，即提供地图。将认识到，使用该(以及前一)实施例，因为有效地断开了数据点，所以可以对数据进行随机采样，所需要的只是测量值和关联(x，y，z)位置的集合，以便能够确定针对任何其他点的集合(尤其是表面或平面)的重力或重力梯度场，也就是地图。此外，可以从中看出，通过传统方法，不需要填充所测量的数据以构建矩形。概念上，可以将数据点视为沿x和y方向由针对傅立叶展开的最大波长所确定的最大尺寸的二维区域的任意形状的部分(但是，如可以看出的，不需要所测量的数据位于二维表面上)。因而可以看出，可以独立于测量尺寸来选择k_m和k_n的值(这里k是波数，即2π/λ，而利用传统测量，最大波长必须与矩形测量数据的长度或宽度相对应)。

在另一方面，本发明提供了一种处理来自地表的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的数据处理系统，所述系统包括：数据存储器，用于存储所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；程序存储器，用于存储处理器控制代码；以及与所述数据存储器和所述程序存储器耦合的处理器，用于加载和实现所述控制代码，所述代码包括用于控制所述处理器执行以下动作的代码：输入所述测量后势场数据；确定所述势场测量和所述位置之间的多个关系，每个所述关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；以及针对所述多个关系来确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据。

在另一相关方面，本发明提供了一种用于处理来自航空重力测量的测量数据以提供针对重力地图的数据的方法，所述测量数据包括多个重力势场测量，每个重力势场测量均具有三维的关联测量位置，所述方法包括：使用所述测量数据来估计乘以了谐波展开中所述三维位置的函数的系数，从而确定针对所述重力地图的所述测量区域中的重力场的表示。

本发明还提供了一种用于处理来自航空重力测量的测量数据以提供针对重力地图的数据的方法，所述测量数据包括多个重力势场测量，每个重力势场测量均具有三维的关联测量位置，所述方法包括：使用所述测量数据来估计多个质量元素，每个质量元素乘以了所述重力场的等效源表示中的所述三维位置的函数，从而提供了针对所述重力地图的所述数据。

本发明还提供了一种绘制重力场的地图的方法，所述方法包括：输入来自航空重力测量的测量数据，所述测量数据包括多个重力势场测量，每个重力势场测量均具有三维的关联测量位置；确定所述势场测量与所述位置之间的多个关系，每个所述关系将所述势场测量与乘以了场地图绘制参数的三维的所述关联位置的函数相关；以及通过针对所述多个关系确定所述场地图绘制参数的大致自洽集合来绘制所述重力场的地图。

上述方法的实施例还包括：通过沿并不局限于平行或定义了矩形网格图形的路径组飞行，来测量该场，以获得测量的势场数据。

如上所述，如所期望的，可以利用有效随机采样的数据来使用上述方法的实施例。尽管随机采样可能实际上并不方便，但是将会理解，由于所描述的方法通常是不受限的，因而尤其由于出于一般的安全考虑和航空器能力而限制的尽可能低的飞行，可以选择飞行调查路径以获得更好的初始数据。因而，例如，可以通过另外沿山谷沿线飞行和/或通过使一条或多条飞行路径弯曲以拟合山脉进行修改，补充广义上传统的路径组。

尽管不需要路径大致在相同高度上交叉(“相交”)，但是相交对于降低噪声和低频漂移是有用的，因为相交提供了类似区域的两个密切相关的测量。然而，不同于传统测量，不需要路径在相同的高度上交叉。此外，例如在山谷中许多交叉图形也是可用的，两个分段线性或蜿蜒的路径可以是沿相反方向呈Z字形扭曲，从而提供沿路径的连续交叉。

根据本发明的另一方面，提供了一种用于进行航空势场测量的方法，所述方法包括使航空器沿路径组飞行，并测量所述路径上的点处的势场数据，其中所述路径组具有以下特征的一个或多个：两个路径在相差至少50米的高度上交叉；在所述测量的区域中，沿相同总方向的路径不平行超过5度；所述路径包括曲线路径；所述路径组中的所述路径总地并不大致位于一个表面上；所述路径组中的所述路径总地定义了表面，其中所述路径中的至少一条路径定义了所述表面中的两个正交方向之一，从而所述表面具有与另一正交方向的距离在高度上大于百分之五的改变率。

在上述方法的实施例中，两条路径在相距大于50、100、150或200米的高度交叉。这些路径可以是线性的、分段线性的或曲线的，并且通常相邻的飞行路径可以非平行多于两度、三度、五度、十度或更多。如先前提及的，这些路径总地并不需要大致在一个表面上，尤其由于没有“大致恒定高度”这一限制。然而，针对这种航空测量来构造路径组的方便方法是构造传统测量，然后对其进行修改，以更加接近于特定区域(如期望有更好的覆盖的山谷)中的地面。

构造传统测量的典型方法是在下层地形上执行要在其上执行测量的表面的二维覆盖。在下层地形的顶点，表面高度典型由航空器的最小允许飞行高度(为了安全)来确定，然后根据航空器确定的攀升/下降速率的限制来降低表面高度，典型为百分之二至三的数量级。由于通过矩形网格来限制路径，所以这些限制作用于两个正交方向。相反，利用上述方法的实施例，这种限制仅需要作用于一个方向(一维而非二维)。有效地，一维覆盖可以用于确定飞行路径。因而在所述方法的实施例中描述了路径组，这些路径组总地定义了定义一个或两个正交方向的所述路径所位于的表面，沿另一方向，允许该表面具有其距离大于所允许的航空器攀升/下降速率的高度改变率，例如，大于百分之三、百分之五、百分之十、百分之二十或更多。

本发明还针对这种航空势场飞行测量路径组提供了载有航空器导航数据的数据载体。

本发明还提供了一种处理器控制代码，用于尤其在诸如盘、CD-或DVD-ROM之类的数据载体、诸如只读存储器之类的编程存储器(固件)、或在诸如光或电信号载体之类的数据载体上实现上述方法。用于实现本发明的实施例的代码(和/或数据)包括诸如C之类的传统编程语言(解译或编译)的源、对象或可执行代码、或汇编码、用于建立或控制ASIC(特定用途集成电路)或FPGA(现场编程门阵列)的代码、或者用于诸如Verilog(注册商标)或VHDL(极高速集成电路硬件描述语言)之类的硬件描述语言的代码。如本领域技术人员将会理解，可以在彼此通信的多个耦合组件之间分布这种代码和/或数据。

本发明的其他方面如下：

一种从空中进行势场测量、从而使航空器配备有包括一个或多个势场测量装置(例如，矢量重力计、重力计、磁力计、磁重力梯度计或其他)的一种地理测量设备的方法，其中航空器以覆盖了测量区域的方式(尤其利用针对飞行的每条线路距离航空器的最低可能表面，并利用在测量区域上近似统一的覆盖，其中这种线路符合航空器的安全操作)飞行了断开的、非水平的、非直线的飞行线路的无规律集合。优选地，在该方法中，使用沿测量线路采集的、与所测量的势场测量装置数据一起起作用的算法来确定感兴趣的地表势场量，而无需对数据进行分层、使数据在共同水平面上或将数据分格，这种数据包括势场量自身的测量、以及(瞬间)位置和可选的势场测量装置的状态。

以上描述了从空中进行势场测量的方法，但是其飞行图形随实现每个单元区域内可接受的“相交”数量的目标而轻微改变，这里“相交”是不同朝向的航空器所飞行的线路象征性在空间的相同点处相交的点。

在以上描述的从空中进行势场测量的方法中，故意使来自实质上平行但不一定是直线的一组线路的一些线路曲折，从而强制与平行飞行的组中的其他线路大量相交。

在以上描述的从空中进行势场测量的方法中，以在总测量区域(给定对总路径长度/飞行时间的限制)上强制最大(可实现)数量的相交为目标，飞行所有线路。因而没有测量线路必须与任何其他测量线路平行，没有线路必须是直线的，并且通常大多数线路将不会从一侧至另一侧横跨整个测量。

在以上描述的从空中进行势场测量的方法中，航空器飞行了更加传统的测量图形，如2-D覆盖或恒定高度测量图形或包括线路和以在每个相交处交叉为意图而飞行的连结线的一些其他图形。(例如，导航数据可以定义确切的交叉，但是实际中飞行员通常将仅在诸如10m或20m的容限内实现确切的交叉)。

在以上描述的从空中进行势场测量的方法中，在处理之前对数据进行“分级”。这里的分级是覆盖了包括以下的一个或多个的技术的一般性术语：降噪、去除低频漂移、匹配相邻线路的低频内容、将数据与固定高度的航空器相关等。在处理之前也可以将数据分格。

在以上描述的从空中进行势场测量的方法中，分析的第一阶段包括去除或调整实际无法通过被调查的地理结构或测量数据的地形学(地形)产生的数据。地形的影响也可以用于在去除地球的重力势能之前的任何处理阶段修正数据。

在以上描述的从空中进行势场测量的方法中，使用LIDAR(激光雷达)和IMU(惯性测量单元)的组合结合DGPS(差分全球定位系统)来产生精确的DEM(数字正视图模型)，从而针对航空器运动来修正LIDAR数据。DEM和DGPS数据也可以用于针对地形来修正测量的势场数据。类似于航空器加速度，状态、角速率和角加速度数据也可以用于修正势场装置的输出数据。任何车载或远程传感器可以用于提供针对航空器和/或势场装置的位置和运动信息。

因而优选地，航空器配备有大量附加标准航空地理测量装置中的任何装置，如GPS、DGPS、高度计、海拔测量、压力测量、超球面扫描仪、电磁测量(EM)系统、时域电磁系统(TDEM)、矢量磁力计、加速度计、重力计和包括其他势场测量设备的其他设备。

在以上描述的从空中进行势场测量的方法中，例如根据当前的最佳实践，使用固定或可移动基站中的装置，对测量平面上的装置的输出进行修正。这种设备可以包括GPS和磁装置以及高质量陆地重力计。

在以上描述的从空中进行势场测量的方法中，将根据上述方法中的任一方法所采集的数据与任何基于地面或卫星的测量数据组合，以有助于改进分析，这种数据包括地形、谱、磁或其他数据。

附图说明

以下将参照附图，通过示例对本发明的这些和其他方面进行进一步描述，其中：

图1示出了具有飞行测量数据的航空器以及配置用于实现根据本发明方法实施例的数据处理系统的示例；

图2示出了用于处理测量的势场数据以实现根据本发明的方法实施例的过程的流程图；

图3示出了根据本发明方法实施例的用于进行航空势场测量的飞行路径数据的产生过程的流程图；以及

图4示出了针对航空势场测量的飞行路径的示例组，其中的数据可以根据本发明的实施例进行处理。

具体实施方式

背景

首先描述一些背景技术，有助于理解本发明。

势场数据包括但不限于重力计数据、重力梯度计数据、矢量磁力计数据和真磁梯度计数据。通过管理这些量如何根据空间变化、以及不同类型的测量如何相关的一系列关系，在数学上将这些数据特征化。可以通过标量来导出势场的元素和表示。

针对重力，相关势能是重力标量势Φ(r)，定义为

$\mathrm{\Φ}\left(r\right)=\∫\∫\∫\frac{\mathrm{G\ρ}\left(r^{\′}\right)}{|r-r^{\′}|}{d}^3{r}^{\′}$

其中，r，ρ(r’)，G分别是重力场测量的位置、位置r’处的质量密度、以及重力常数。重力(所体验到的重力场)是标量势的空间导数。重力是方向已知(重力向下作用)的矢量。对于任何所选的笛卡尔坐标系，通过三个分量来表示重力：

$g=(g_x,g_y,g_z)=(\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂x},\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂y},\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂z})$

对于这三个分量，每个分量沿三个方向中的每个方向而改变，因而产生了九个量，形成了重力梯度张量：

$G=\left(\begin{array}{ccc}{G}_{\mathrm{xx}}& G_{\mathrm{xy}}& G_{\mathrm{xz}}\\ G_{\mathrm{yx}}& G_{\mathrm{yy}}& G_{\mathrm{yz}}\\ G_{\mathrm{zx}}& G_{\mathrm{zy}}& G_{\mathrm{zz}}\end{array}\right)=\left(\begin{array}{ccc}\frac{\∂}{\∂x}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂x}& \frac{\∂}{\∂x}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂y}& \frac{\∂}{\∂x}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂z}\\ \frac{\∂}{\∂y}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂x}& \frac{\∂}{\∂y}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂y}& \frac{\∂}{\∂y}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂\mathrm{zx}}\\ \frac{\∂}{\∂z}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂x}& \frac{\∂}{\∂z}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂y}& \frac{\∂}{\∂z}\frac{\∂\mathrm{\Φ}\left(r\right)}{\∂z}\end{array}\right)$

势场的数学原理和基本方程以及关系遵循对标量势函数、其导数、其傅立叶变换和其他数学量的属性的分析。

根据格林定理，在已知闭合表面上的标量势的空间导数(包括标量势本身)的情况下，在该表面所包围的体积内的所有点上已知该空间导数值。其推论是，一旦该量对于所有点都是已知的，则通过微分和积分，可以导出包括标量势自身的标量势的所有其他导数。因而当在包围了该体积的表面上仅已知导数之一时，则有效地知道了该体积内的所有点上的标量势和所有导数。这指示了，标量势的任何导数的任何分量的全部测量允许计算标量势的任何导数的任何其他分量。据此得出(至少在理论上)测量哪个量并不重要，装置的选择简单地归结于哪个装置以最大信噪比测量了期望量。

以上重力标量势的微分(在一些与

关联的问题之后，其中r→0)最终产生：

${\▿}^2\mathrm{\Φ}\left(r\right)=\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂x^2}+\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂y^2}+\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂z^2}=-4\mathrm{\πG\ρ}\left(r\right)$

该方程在不重要的区域内变为拉布拉斯方程，得到重力的重要基本关系：

${\▿}^2\mathrm{\Φ}\left(r\right)=\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂x^2}+\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂y^2}+\frac{{\∂}^2\mathrm{\Φ}\left(r\right)}{\∂z^2}=0$

谐波函数满足拉布拉斯方程，具有可以在对从势场测量中采集到的数据进行分析的过程中所使用的许多属性。

可以使用多种技术来分析和处理数据，这些技术与从测量中采集到的数据一起起作用，作为起点，但是之后改变数据和/或数据格式，因而与所测量的量关联的值全部出现在位于水平的、固定海拔的分析平面上的规则2-D网格上，将所使用的过程称为“分级(levelling)”和“分格(gridding)”。

分格是地球物理学数字处理技术，其基本思想是：

·将测量区域分为其侧面通常与进行测量的主要方向对齐的矩形小区，

·利用完全“等效于”测量数据的数据来替换实际测量数据，该测量数据现在分配有每个小区中部的点处的值。

这种数据被称为“分格”数据。存在许多方式以确保代表测量数据的“所产生的”数据确实“等效于”测量数据。也存在许多方式来选择每个小区的尺寸，但是这些方式全部与沿两个正交方向飞行的线路的平均间隔相关。在没有以这里所描述的规则方式沿测量线路飞行时，尽管不是必须，仍可以使用分格的概念，但是对于小区大小和方位的选择存在更大的范围。

一旦数据处于该格式，则在数学上容易处理得多。不要求这些过程中所涉及的任何数据满足拉布拉斯方程。将该数据视为数值和统计集合，将其他方法应用于这些数值以给出针对水平分析平面的势场的最佳估计。

可以将数据变为2-D傅立叶序列，因而除了需要数据位于固定高度之外，还需要每行数据具有2n个数据点，以使快速傅立叶变换方法起作用，并且沿正交方向采集数据。

表示重力梯度的2-D傅立叶序列的一般形式可以表示为以下形式的2-D空间正弦波之和：

$g_{\mathrm{zz}}(x,y)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}{g}_{\mathrm{zz}}(k_m,k_n)e^{{\mathrm{ik}}_{m}x}{e}^{{\mathrm{ik}}_{n}y}$

其中波数k_m，k_n与测量的尺寸相关，x、y方向上的L_x，L_y分别由以下表示：

$k_m=\frac{2\mathrm{\π}}{L_x},$ $k_n=\frac{2\mathrm{\π}}{L_y}$

该表示仅在恒定高度上有效。然而后面将描述了波数和系数g_zz(k_m，k_n)的高度相关性，这共同形成了重力梯度的2-D傅立叶变换。

优选实施例的详细描述

如果利用了重力或磁势场的属性，则不必如上处理数据。以下的讨论将基于重力场，但是可以容易地扩展至磁场。

首先陈述与重力场相关的一些观察：

1.在各个方面可以对主体外的重力场进行建模，仿佛它来自完全位于主体表面难以察觉的薄层的物质，并精确地符合主体的表面。这种层定义了等效源，即产生大致(理论上精确地)与主体本身相同的重力标识的重力源。可以有许多定义等效源的方式，可以是如上所述的方式，或者可以是严格水平的，这些等效源可以全部或部分在地球表面之上或之下、或者其他，但是具有相同的属性；即，这些等效源产生与如地球所产生的重力场相同的重力场。其他信息可以参照R.J.Blakely，“Potential Theory in Gravity and Magnetic Applications”，Cambridge University Press，1995。

2.主体外的重力场也可以唯一地写作修改后的2-D傅立叶基本函数序列，其在主体外的所有位置上都是有效的。该函数明确地是调和的，并在以下再现。其形式完全通用，并且可以表示主体内任何可能的质量分布。

在优选实施例中，使用重力梯度计测量G_zz作为位置r_measure的函数，并且在不需要产生重力梯度张量的其他元素的情况下与其一起起作用。这可以用于产生下层质量分布的表示。这里讨论可以使用的两种技术，等效源技术和修改后的傅立叶基本函数技术。上述提及的格林定理示出了原理上可以根据G_zz导出下层质量分布，但是在所描述的技术上，并不明确地使用格林定理。

等效源方法

在该实施例中，将测量区域的表面分为小块，典型地是一侧50m的量级，这些被称为小板(后面将其称为质量元素)。易于前向计算每个小板的重力(如Blakely(同上)所述)，调整其质量直至获得对测量数据的最佳整体拟合。通过标准最小平方拟合过程来做出该质量确定。通过将真实测量位置处的数据与在相同的真实测量位置处的所提出的等效源产生的重力场匹配来获得该拟合。该过程在数学上是精确的，不会将任何人为调整引入数据以与水平矩形测量一致。

一旦获得了拟合，便认为该拟合是最初的数据集合。所有用于确定地质结构的后续分析比较任何给定地质结构将会产生的重力场与由等效源所产生的重力场之差，并将该差最小化。这种技术的一个实际优点在于，最佳拟合来自于质量分布(尽管是合成的)，因而最佳拟合情况将自动满足拉普拉斯方程。这不同于上述的传统方法，上述传统方法产生数值最佳拟合，但是并不施加该数据必须满足拉普拉斯方程的附加限制，即数据来自真实的质量分布的限制。

等效源方法不必使用与地形共形的表面，可以使用可以覆盖位于恒定高度、在地球真实表面之上或之下的任何表面、可以穿过地球真实表面等的源。在各个小板的质量中，符合该地形的表面的选择可能产生较小的变化，但是原理上整个结果并不会受到表面的任何合理选择的显著影响。

同样，小板可以是任何大小或形状的，甚至不需要是大小相等的，事实上，如果允许小板的大小和几何形状依据每个区域内的地形变化快慢而改变，则有助于提高分析效率。仅由所使用的小板个数以及在测量的每个位置上在分析中使用的这些小板的个数来确定该过程的数学复杂度。该技术的一个优点在于，对于重力或重力梯度的一些分量，可以只使用数据点区域中的小板，这极大地减小了分析的复杂度。很清楚，该分析是空间域中的分析。

一旦产生了等效源，则可以通过直接前向计算来预测任何表面上的重力标量势的任何导数。从分析和可视化角度来看，该过程是有用的。

更具体地，给定每个表面元素的质量，使用直接计算来预测在每个测量点处针对所测量的量将会获得何值，该值是表示重力矢量的分量或重力梯度张力的分量的量。通常，这将会是以下示出的形式之和。这里，使用gg作为所测量的量的符号，如上所述，在一些优选实施例中所测量的量是G_zz。

${\mathrm{gg}}_{\mathrm{calculated}}\left(r_{\mathrm{measure}}\right)=\underset{\mathrm{all}-\mathrm{masses}}{\mathrm{\Σ}}{m}_{\mathrm{mass}-\mathrm{element}}F(r_{\mathrm{measure}}-r_{\mathrm{mass}-\mathrm{element}})$

方程(1)

在以上的方程中，F被称为格林函数(同上，Blakely，在第185页，在此一并引入作为参考)以及r_mass-element定义了质量元素的位置(例如，重力中心或一些其他定义点)。函数F是标准函数，并且是本领域技术人员已知的(可知针对不同的质量元素几何形状推出)。例如，同上的Blakely在第187页定义了针对矩形棱镜的格林函数F，将其一并引入作为参考；该函数具有8项，每项与棱镜的顶点相对应。

对于该关系的重要事物是，函数“F”是已知的，并且是每个元素的质量。通常，所计算的值不会与所测量的值相同，在一种方法中，构造以下和

$S=\underset{\mathrm{all}-\mathrm{measurements}}{\mathrm{\Σ}}\left[{\left({\mathrm{gg}}_{\mathrm{calculated}}\right)\left(r_{\mathrm{measure}}\right)-{\mathrm{gg}}_{\mathrm{measured}}\left(r_{\mathrm{measure}}\right))}^2\right]$

方程(2)

有效地，量S定义了总平方误差S(噪声和建模误差)，并且通过将其最小化，可以找到最佳拟合gg_calculated。这可以通过标准过程，例如使用Matlab(RTM)来实现。

通过对上述自变量的展开，针对量S，唯一未知的是每个质量元素的实际质量，以及与每个点处的测量相关联的测量噪声。通过要求S相对于每个元素的质量的变化是最小的，获得最佳拟合(即这些质量元素的值的最佳估计)。当获得修正系数时，量S仅是与每个测量相关联的平方噪声项之和，这是可以通过S获得的最小值。通过相对于每个质量元素的质量对S进行微分，针对每个质量元素，相对于m_mass-element的S的一阶导数必须为零。这最终导致了产生联立方程组，其中唯一的变量是每个表面质量元素的质量，以及等式个数等于拟合过程中所涉及的质量元素的个数。方便地，该联立方程组可以通过矩阵表示，并针对表面质量元素值进行求解。由于测量数量通常远大于联立方程的个数(即，质量元素的个数)，例如至少大五或十倍，因而该方程组是过约束的。这是有用的，因为如果不是这种情况，则表面质量元素的值将与噪声拟合；作为替代，使用例如最小平方拟合的优化过程来确定表面质量元素的值，因而减小噪声。此外，由于使用了物理模型重力场(实质量元素)，所以该方程符合拉普拉斯方程，从而意在拟合的是重力场而非噪声。

实际中，例如可以通过根据距离对质量元素加权、并将超过阈值距离的权重设置为零，来有效地忽略位于大于所定义的距离处的质量元素。该阈值距离典型为几千米，例如在1至10千米的范围内，但是在一定程度上取决于地形(例如，也许需要扩展该距离以包括大的附近的山脉)。对质量元素进行加权以忽略位于大于阈值距离的那些质量元素的结果是，表示所计算的测量值之差的矩阵成为稀疏矩阵。这是有用的，因为这种矩阵可以包括例如100K乘500K个数据点，如果这些大部分为零，则可以显著减小处理负担。

本领域技术人员将会理解，上述定义的量S仅是将所计算的数据与所测量的数据拟合的一种方式，存在许多其他方式。一旦定义了针对质量元素的等效表面的平面，自由参数便是质量元素的质量(由于距离函数已知)；本领域技术人员将会理解，可以用以拟合这些自由参数的许多数值技术包括例如迭代数值搜索。在实践中发现有用的一种技术是使用没有平方函数、而采用求模的方程2的版本。

一旦发现了质量元素，前向计算(即对质量元素的作用求和)便可以导出其他标量势分量，并且通过微分可以导出G的其他分量。G的导出值可以与地质模型进行比较(被称为“解译”)，以确定下层的地质结构。

修改后的傅立叶基本函数方法

在所述方法的第二实施例中，使用修改后的傅立叶基本函数技术。在该方法中，使用二维傅立叶展开来表示所测量的数据，但是将重力场的2-D傅立叶变换的每个分量

的每个二维系数A_mn重写为位置无关系数的乘积，(C_mn)乘以高度相关系数

即，明确地，可以将任何质量分布的重力场写为以下形式：

$\mathrm{gg}(x,y,z)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}{C}_{\mathrm{mn}}{e}^{{\mathrm{ik}}_{n}x}{e}^{{\mathrm{ik}}_{n}y}{e}^{-\left|k_{\mathrm{mn}}\right|z},$ $\left|k_{\mathrm{mn}}\right|=\sqrt{({k_m}^2+{k_n}^2)}$

在傅立叶变换语言中，A_mn是重力场的2-D傅立叶变换。该表示完全通用，并且表示重力场的展开的每个分量明确地符合拉普拉斯方程(当微分时，高度相关系数的乘数与2D分量的乘数抵消)。现在明确的是每个傅立叶系数的高度相关性，还清楚了通过知道C_mn系数所提供的重力场表示可以完全确定在没有质量的空间中任何点处的重力场。可以依据该表示的期望精度，在一点对该修改后的傅立叶序列所暗含的展开进行截短。

在位置r处，已知测量的势场(重力或重力梯度)，针对被称为波数的特定值k_m和k_n对，每个修改后的傅立叶系数基本函数的值

也是已知的。每个波数描述了傅立叶基本函数的谐波空间变化，例如，给定波数k_m＝2π/λ_m的每个系数在空间内以时段λ_m表现周期性。针对每个独立的方向，傅立叶方法的另一属性是，每个波数将会是与广义上等同于测量尺寸的波长相对应的最小波数的整数倍。

可以采用其他展开，例如球谐函数(对于较大区域有用，尤其在地表的曲率很大的情况下)。

广义上，以与上述等效源方法实际相同的方式来确定展开系数。一旦选择了每个方向所感兴趣的最大波长，以及用于测量的坐标系的方位和原点的最大波长，便可以在空间中的每个点处确定基本函数值。一旦这样做了，便针对所测量的量计算相同的差方函数S，并通过将这种函数相对于C_mn中的变化最小化来确定针对系数C_mn的最佳值。因此，如上述第一实施例一样，产生用于求解傅立叶展开系数的联立方程组。然后可以采用这些系数来确定(通过前向计算)在任何二维表面上的重力或重力梯度，即提供地图。将会认识到，使用该(以及先前的)实施例，可以对数据进行随机采样，因为有效地断开了数据点，所以可以对数据进行随机采样，所需要的仅仅是测量值和关联(x，y，z)位置集合，以便能够确定针对任何其他点集合或实际平面上的点集合的重力或重力梯度场，这就是地图。此外，可以从中看出，不需要填充所测量的数据以构建矩形。概念上，可以将数据点视为沿x和y方向由针对傅立叶展开的最大波长所确定的最大尺寸的二维区域的任意形状的部分(但是，如可以看出的，不需要所测量的数据位于二维表面上)。因而可以看出，可以独立于测量尺寸来选择k_m和k_n的值(这里k是波数，即2π/λ，而利用传统测量，最大波长必须与矩形测量数据的长度或宽度相对应)。

更具体地，傅立叶序列路径之后的概念是在傅立叶域而非空间域中求解方程(1)的问题。在下文中，将会假设g_zz是所测量的量(如上也称为G_zz)。

可以将任何主体或质量分布中的重力场写为以下的完全通式：

$g_{\mathrm{zz}}(x,y,z)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}{g_{\mathrm{zz}}}^{\′}(k_m,k_n)e^{{\mathrm{ik}}_{m}x}{e}^{{\mathrm{ik}}_{n}y}{e}^{-\left|k_{\mathrm{mn}}\right|z}$ 方程(3)

其中，

$\left|k_{\mathrm{mn}}\right|=\sqrt{{k_m}^2+{k_n}^2}$

在方程(3)中，系数g_zz′(k_m，k_n)在形式上与高度和波数无关，该系数g_zz′(k_m，k_n)仅取决于测量下的地质。应注意，方程(3)的形式与以上给出的2-D傅立叶序列表示的通式相似。仅在

$g_{\mathrm{zz}}(k_m,k_n)={g_{\mathrm{zz}}}^{\′}(k_m,k_n)e^{-\left|k_{\mathrm{mn}}\right|z}$

时，表达式相同。

可以看出针对高度和标准傅立叶变换的波数(即系数g_zz(k_m，k_n))的强相关。这是无法实施在可变高度测量上采集的数据的标准傅立叶变换的原因。也示出了实际上无法根据点来“向上”延拓点的数据；向上延拓仅在已知量的精确频率内容时才起作用。还可以看出，方程(3)展开的每一分量具有不同的高度相关性。

然而，所描述的修改后的傅立叶基本函数方法本身并不试图执行傅立叶变换，如上所述，该方法不是在数学上合理的方式。

该过程是，类似地处理每个数据点，可以向量

分配值，该值在已知测量位置时是精确的。这是由于如后面所述选择k值，并且已知r(x，y，z)。例如，通过概念上将矩形放置于测量区域周围(或略微在测量区域内)(近似表示测量区域)来选择k_x和k_y的最大和最小值，矩形的尺寸定义了最大波长(k值)。例如，可以依据测量的近似平均高度来选择最小值，例如，0.5至3之间的多个平均测量高度。这是由于变化的长度规模(length scale)取决于测量高度，并且广义上，选择最小值来以大于阈值的重要性表示在平均测量高度上期望或看到的变化(考虑了计算能力)。为了示出该概念，z＝测量高度处的exp(-kz)项对于波长等于测量高度的情况衰减了大约e-6。

然后，使基于方程(3)的方程组等于所使用的数据点个数，从而已知每个数据测量点处的重力梯度的测量值，包括噪声和量

根据这一点，唯一未知的是系数g_zz′(k_m，k_n)，然后该拟合过程清楚地成为获得该系数的最佳估计的过程。

参照方程(3)，其给出了计算后的g_zz，通过与方程(1)类比，可以导出方程(2)的一个版本。这提供了等于未知数个数的方程组，在这种情况下，未知数是系数g_zz’的值。

按照与等效源方法相同的方式来实现上述过程，即通过最小平方拟合过程来确定系数。重要问题在于，以这种形式，系数g_zz′(k_m，k_n)明确地不具有高度或波数相关性。

该技术具有与等效源分布类似的优点，即该解决方案是调和的，因为方程(3)的展开中的每一项明确为调和的。装置噪声不是调和的，因而甚至单独地，该技术也具有去除该过程中大部分噪声的可能。

类似于等效源方法，该技术不需要规则网格数据。该技术的一个附加明显优势在于，不同于傅立叶序列方法本身，该技术并不依赖于可用于分析的波长。这是取决于针对可能感兴趣的区域上的完全通项中的修改后的傅立叶序列方法的有效性的巧妙结果。该过程包括选择沿测量平面的两个正交方向的所感兴趣的最大波长，这导致了沿每个方向的所感兴趣的最大波长λ_xmax，λ_ymax。通过下式，该波长表明其自身是在每个方向上的最小波数：

$k_{x\min }=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_{x\max }},$ $k_{y\min }=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_{y\max }}$

方程(3)的展开中的所有波数是它们的整数倍。即：

k_m＝m·k_xmin，k_n＝n·k_ymin

此外有利地，该方法不必以2的幂作为要在每行的分析中使用的数据点的个数，而这是快速傅立叶变换方法中所需要的。这避免了填充数据以实现这种情况的需要。

不同于等效源方法在空间域中进行分析的情况，在频域中进行傅立叶序列分析。频率表示作为给定分量在每个距离测量单位上所具有的周期数的测量的空间频率，这通过以上引入的波数来表示，波数越多，每个距离单位的周期数越多，重复距离越短。所提出的技术并没有在任一域中明确地起作用，该技术仅产生了一方程组，从中可以计算未知数的最佳估计。

在实际项中，傅立叶方法通常需要矩形区域上的测量，最小波数是可以用于表示数据的最长波长，最长波长明确地是每个方向中的测量尺寸。在测量并不符合该标准的情况下，通常创建数据，使得修正后的数据集在被称为填充的过程中满足该标准。

修改后的傅立叶基本函数方法与传统的傅立叶分析相比具有许多优点。例如：

1.根本不必执行任何傅立叶变换，已知将会是何种标准格式，方程(3)中的展开系数要经由最佳拟合过程而非明确计算傅立叶变换来确定。

2.要使用的波数不必是与测量的大小和形状相关的整数。这是由于，已知傅立叶变换的形式，经由刚才描述的最佳拟合过程来确定系数。

也许最显著的结果在于在非恒定高度上采集测量数据。传统上必须在恒定高度上执行2-D傅立叶变换，并且在方程(3)中引入的傅立叶变换是强高度函数。在方程(3)中明确地表示了该功能，系数g_zz′(k_m，k_n)不再是高度相关的。

总结

所描述的方法实施例允许完全利用实际测量数据，并且在数学上是精确的；所述方法不会通过使数据拟合于规则恒定高度网格而破坏数据，所产生的解决方案是调和的，数据并不必以2的幂作为两个正交方向上数据点的个数。

通常，该过程是使用来自每个数据点的所有数据，并试图产生等效表面质量分布、或者在测量区域上最佳拟合数据的修改后的傅立叶序列系数集。

这基本上不同于假定数据来自固定高度(即使不是这样)、并根据所测量的数据产生新数据、以在规则2-D网格上产生数据的处理方法。该“分格”数据集将值分配给所测量的数据，这些所测量的数据不是被测量的数据，并被分配给航空器不采集数据的位置，但这是不正确的。

通过这种过程产生的误差可能在给定情况下固定高度近似不好时更加严重。如果试图尽可能靠近地面飞行以将来自相对靠近地球表面的小特征的信号最大化，则误差也是严重的，这些特征导致了修改后的傅立叶序列展开中的短波长分量的产生。最后，在装置噪声很大时这些误差不太严重，因为在外插中涉及到的误差小于与装置相关联的误差。

现在参照图1，示出了进行势场测量以获得数据从而根据如上所述的方法进行处理的航空器10的示例。航空器10包括在其上安装了重力梯度计14的惯性平台12，重力梯度计14将势场测量数据提供给数据采集系统16。惯性平台14安装有惯性测量单元(IMU)18，该惯性测量单元18也将数据提供至数据采集系统16，典型包括状态数据(例如，斜率、滚动和偏航数据)、角速度和角加速度数据、以及航空器加速度数据。航空器也配备有差分GPS系统20和LIDAR系统22或类似装置，以提供航空器在下层地形上的高度的有关数据。航空器10也可以配备其他装置24，如磁力计、TDEM系统和/或超谱成像系统，再次馈入数据采集系统。数据采集系统16也具有来自一般航空器装置26的输入，该一般航空器装置26可以包括例如高度计、航空和/或地面速度数据等。例如，数据采集系统16可以提供一些初始数据预处理，以便针对航空器运动修正LIDAR数据，和/或将来自IMU18和DGPS 20的数据组合。数据采集系统16也具有通信链路16a和/或非易失性存储器16b，从而能够存储所采集的势场和位置数据以用于以后的处理。也可以提供网络接口(未示出)。

通常(但不必定)离线地、有时在不同于采集测量数据的地区的地区执行数据处理，以产生用于势场测量的地图数据。如图所示，数据处理系统50包括处理器52，该处理器52与代码和数据存储器54、输入/输出系统56(例如包括网络和/或存储介质和/或其他通信的接口)、以及例如包括键盘和/或鼠标的用户接口58耦合。可以在可拆卸存储介质60上提供存储在存储器54中的代码和/或数据。在操作中，数据包括从势场测量中采集的数据，代码包括例如根据图2所示的过程(在下面讨论)来处理该数据以产生地图数据的代码。

现在参照图2，图2示出了在数据处理器上实现的过程示例，在实施例中，该数据处理器可以包括通用计算机系统，用于根据先前描述的技术来处理来自飞行测量的数据。因而在步骤S200，该过程输入测量的势场数据(例如重力梯度)和关联的3D位置数据。可选地，在步骤S200a，例如可以进行一些预处理，以去除异常和/或减少(或增加)或选择要处理的数据。

在步骤S202，该过程构造了与势场数据、3D位置数据和场地图绘制参数集相关的联立方程的过约束组。更具体地，该过程构造了将每个势场测量与乘以了场地图绘制参数的关联3D位置(s(r-r’))的函数相关的联立方程。优选地，该过程使用遵循拉普拉斯方程的势场表示(更具体地，位置函数)，例如上述等效源或修改后的傅立叶基本函数。可选地，在该阶段可以包括地形修正数据(S204)。例如，该过程通过传统的最小平方方法对这些联立方程进行求解(S206)，从而确定场地图绘制参数的大致自洽集合。这些有效地定义了被测量的势场的地图，因而构成了该过程的一个优选输出。然后，(S208)可以使用这些参数来执行前向计算，从而产生表面上势场的期望地图，例如，近似地表区域的表面(如所需要和当需要时)。

针对传统的航空重力，不必非常靠近地表飞行，因为这些装置相对不灵敏，并且还被设置为仅对于长波长敏感。然后，可以在差不多恒定的海拔处、或在可以合理假设数据来自恒定海拔的高度上进行测量。

然而，对于使用新的高灵敏度重力梯度计系统的任何测量(也能够分辨短得多的波长)，期望尽可能接近地面飞行，以便将信号(尤其包含在短波长分量中的信号)最大化。这种测量不可能足以满足恒定海拔标准，但是当前的所有分析方法就像这种测量满足标准一样来进行。

然而，可以利用等效源或修改后的傅立叶变换方法来处理数据，并将其与对于势场测量而言完全新颖的飞行方法进行组合。

就最一般的形式而言，该方法飞行完全断开的测量线路组，从该线路采集数据，并使用该数据和获取的绝对位置，获得等效源小板质量或修改后的傅立叶变换方法的分量的最佳拟合参数。

测量线路不需要是直线，不需要是水平的，不需要沿正交方向，沿大致正交方向的线路不需要交叉。理想地，飞行线路将尽可能符合地形，而无需任何其他限制。当然，位于该测量飞行方法和那些包括恒定海拔正交相交飞行线路的方法之间的任何其他飞行方法都是可能的，并且任何测量可以包括任一飞行方法的元素。解析技术的能力和数学修正使得这些飞行过程成为可能。

清楚地，可以在使用上述过程之一之前或之后另外处理使用该飞行方法采集的数据。可以利用针对噪声降低或数据组合的任何现有技术。通过使用高质量DEM可以极大地消除地形影响，实际上可以使用与更传统的数据处理相关的任一过程。

航空器可以是适于低飞的任何地球物理测量航空器。实际上，航空器将配备有与其符合安全和合理飞行时间而可以携带的那样多的设备。所选设备将测量与地质相关的信号和与航空器动作相关联的信号。尤其在重力相关测量中，使用与航空器动作相关联的信号来修正由选择装置所获得的测量以用于测量地形属性。

图3示出了用于产生针对势场测量定义了飞行图形的飞行测量数据的过程的流程图，可以使用上述方法来处理该过程的结果。

在步骤S300，该过程输入地形数据，并且例如通过用户接口，输入定义了要测量的期望区域的数据。然后，在步骤S302，该过程在所选地形部分上覆盖一表面，以确定包括飞行路径组的飞行测量图形。表面的覆盖可以受到表面内任一方向上的最大斜率的限制，并且最初，飞行路径可以定义矩形图形。然而在步骤S306，该过程调整表面和/或飞行路径，以使路径适合表面。这受到最大攀升/下降速率(根据航空器数据而确定)的限制，但是仅需要沿路径施加该限制。因而在其他方向上，例如与路径正交的方向上，表面的斜率可以超过由该限制所定义的值。此外或可选地，在步骤S308，该过程调整表面和/或飞行路径，以使路径更加接近于下层地形的表面。这可以通过允许路径是以下的一个或多个路径来执行：曲线、非平行、在交叉点处高度上明显不同(例如，多于50米、100米或200米)、位于其斜率大于航空器的最大攀升/下降速率的表面内、或者不在单个表面内。一旦执行了该适应，该过程便输出定义了所确定的飞行路径组的飞行测量数据(S310)。该数据可以由飞行员用以根据所确定的测量图形使航空器飞行。

将会理解，不同于传统的分格系统，这种技术使得能够靠近谷底、沿山谷定义飞行路径，在靠近谷底的情况下，其他技术将会由于网格图形与航空器的最大攀升/下降速率的组合而限制航空器与地面的接近。类似地，将会理解，由于所描述的技术使得能够对来自这种飞行测量的数据进行处理，因而通常可以定义紧贴地形的路径组。其优点在于，根据这种技术获得的数据由于具有在地球表面以上高度的势场扰动的振幅的下降中的显著速率而要准确得多。

图4示出了飞行测量路径的示例，可以通过先前所描述的技术对来自该路径的数据进行处理。图4示出了能够实施这些技术的路径类型的示例。

在平坦区域内，所有飞行策略趋于一致，但是在丘陵或真正的山脉区域，所提出的方法产生了显著的优势。在这种区域内，相对独立地分析测量数据的每个部分，并形成了关于如何在该区域飞行的策略。

在这种情况下，针对地形分析测量区域，并基于所感兴趣的地形目标以及如何可以在该区域上飞行以最佳地识别那些地形目标，来形成飞行测量策略。如果，如通常的情况，需要在所有测量点非常靠近地面，则测量线路将以相对少的连接性而断开。航空器仅可以以给定速率攀升或下降，需要足够的间隔来安全地飞过山脊。向“下”或“上”飞过山谷的线路将不会与“横越”山谷的线路交叉、或者将在较少交叉处相交。通常认为交叉对于除去低频噪声(通常称为漂移)非常重要，这是针对设计测量的另一考虑。

标准测量所期望的设备包括：重力梯度计；重力计；惯性测量单元(IMU)；GPS：LIDAR；磁重力计；磁力计；高度计；雷达高度计；气压计；摄谱仪；超谱成像器；EM；TDEM；导航S/W；数据采集设备。

一旦形成了飞行策略，便进行实际测量，并对数据进行采集。该分析使用来自势场测量系统的真实数据，对于与测量相关联的任何已知人为误差(例如，航空器倾斜、航空器加速度、压力作用等)来修正该数据，然后如上所指示来分析数据。低飞的原因在于，所测量信号(尤其是那些具有短波长、即大波数的信号)的快速高度相关性。

毫无疑问，本领域技术人员将会获得其他有效的可选方式。将会理解，本发明不限于所描述的实施例，包括在所附权利要求的精神和范围内的本领域技术人员所显而易见的修改。

Claims (16)

1.一种处理来自地表的测量区域的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的方法，所述方法包括：

输入所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；

确定所述势场测量和所述位置之间的多个关系，所述多个关系中的每个关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；以及

针对所述多个关系确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据；

其中，所述场通过多个表面元素来表示，每个表面元素具有关联的源强度、大小和几何形状，

其中，每个所述场地图绘制参数包括所述源强度，以及

其中，当所述测量区域具有变化的地形时，选择所述多个表面元素中每个表面元素的所述大小和所述几何形状以与所述变化的地形匹配。

2.如权利要求1所述的方法，其中所述地图数据包括所述场地图绘制参数集合。

3.如权利要求1所述的方法，其中所述地图数据包括绘制在表面上的所述场的地图的数据，所述方法还包括根据所述场地图绘制参数集合来确定所述表面上的所述场。

4.如权利要求1所述的方法，其中所述关系包括联立方程组，其中所述场地图绘制参数的确定包括对所述联立方程进行求解。

5.如权利要求4所述的方法，其中所述联立方程组是过约束的，以及所述求解包括通过降噪过程进行求解。

6.如权利要求1所述的方法，其中通过稀疏矩阵定义所述关系，以及所述场地图绘制参数的确定采用数值技术。

7.如权利要求1所述的方法，其中所述位置函数大致满足拉普拉斯方程。

8.如权利要求1所述的方法，其中所述测量的势场数据包括重力计数据和重力梯度计数据中的一个或多个，以及所述场表示重力场。

9.一种具有计算机处理器控制代码的非易失性计算机可读介质，用于实现权利要求1的方法。

10.一种处理来自地表的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的数据处理系统，所述系统包括：

数据存储器，用于存储所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；

程序存储器，用于存储处理器控制代码；以及

与所述数据存储器和所述程序存储器耦合的处理器，用于加载和执行所述控制代码，所述代码包括用于控制所述处理器进行以下动作的代码：

输入所述测量的势场数据；

确定所述势场测量和所述位置之间的多个关系，所述多个关系中的每个关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；

针对所述多个关系确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据，

其中，所述处理器还被配置为：

通过多个表面元素来表示所述场，每个表面元素具有关联的源强度、大小和几何形状，其中，每个所述场地图绘制参数包括所述源强度，以及

当所述测量区域具有变化的地形时，选择所述多个表面元素中每个表面元素的所述大小和所述几何形状以与所述变化的地形匹配。

11.一种用于处理来自航空重力测量的测量数据以提供针对重力地图的数据的方法，所述测量数据包括多个重力势场测量，每个重力势场测量均具有三维的关联测量位置，所述方法包括：

使用以下表达式，确定针对所述重力地图的所述测量区域中的重力场的表示：

$g_{\mathrm{zz}}(x,y,z)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}{g_{\mathrm{zz}}}^{\′}(k_m,k_n)e^{{\mathrm{ik}}_{m}x}{e}^{{\mathrm{ik}}_{n}y}{e}^{-\left|k_{\mathrm{mn}}\right|z}$

其中，

g_zz′(k_m，k_n)是系数，

x、y和z定义了所述三维中的所述位置，

k_m和k_n为正波数，并且

以及

使用所述测量数据来估计所述系数。

12.如权利要求11所述的方法，其中所述系数形式上与高度无关。

13.一种处理来自地表的势场测量的测量的势场数据以确定用于绘制所述场的地图的地图数据的方法，所述方法包括：

输入所述测量的势场数据，所述测量的势场数据包括定义了多个势场测量和关联位置的数据，每个所述位置定义了三维的所述势场测量的位置；

确定所述势场测量和所述位置之间的多个关系，每个所述关系将所述势场测量与乘以场地图绘制参数的三维的所述关联位置的函数相关；以及

针对所述多个关系确定所述场地图绘制参数的大致自洽集合，从而确定所述地图数据；

所述方法还包括：使用航空器测量所述场以捕获所述测量的数据从而用于确定所述多个关系，以及所述测量包括沿具有以下特征中的一个或多个的路径组飞行：

两个路径在相差至少50米的高度上交叉；

在所述测量的区域中，沿相同总方向的路径不平行超过5度；

所述路径组中的所述路径并不都基本位于一表面上；

所述路径组中的所述路径总地定义了一表面，其中所述路径中的至少一条路径定义了所述表面中的两个正交方向之一，从而所述表面具有与另一正交方向的距离在高度上大于百分之五的改变率。

14.一种具有计算机处理器控制代码的非易失性计算机可读介质，用于实现权利要求13所述的方法。

15.一种使用航空势场测量来测量势场数据的方法，所述方法包括沿路径组使航空器飞行，并测量所述路径上的点处的势场数据，其中所述路径组具有以下特征的一个或多个：

两个路径在相差至少50米的高度上交叉；

在所述测量的区域中，沿相同总方向的路径不平行超过5度；

所述路径组中的所述路径并不都大致位于一表面上；

所述路径组中的所述路径总地定义了一表面，其中所述路径中的至少一条路径定义了所述表面中的两个正交方向之一，从而所述表面具有与另一正交方向的距离在高度上大于百分之五的改变率。

16.一种具有导航数据的非易失性计算机可读介质，所述导航数据定义了如权利要求15所述的航空势场飞行测量路径组。

Images