Gsolve：一个Python计算机程序，用于相对重力测量转换为绝对重力值和重力异常

SoftwareX 7（2018）129原始软件出版物Gsolve，一个Python计算机程序，具有图形用户界面，用于将相对重力测量转换为绝对重力值和重力异常Jack McCubbinea，*，Fabio Caratori Tontinic，Vaughan Stagpoolec，Euan Smithb，格兰特·a澳大利亚澳大利亚地球科学局社区安全和地球监测司大地测量科，澳大利亚首都地区2601，堪培拉378号邮政信箱b惠灵顿维多利亚大学Kelburn，Wellington 6012，新西兰cGNS Science，1 Fairway Dr，Avalon，Lower Hutt 5011，New Zealandar t i cl e i nf o文章历史记录：2016年9月1日收到2017年10月4日收到修订版，2018年关键词：相对重力Python最小二乘地形改正a b st ra ct开发了一个带有图形用户界面的Python程序（Gsolve），以协助相对重力测量的常规Gsolve计算相对重力测量的每个测量点的重力，该测量点参考至少一个已知的重力值。在测量数据的处理过程中，考虑了太阳和月亮的潮汐效应该计算基于最小二乘法公式，其中相对于相对重力观测值和一些提供的重力参考点值，最小化每个测量位置处的绝对重力与重力仪动态相关参数之间的差异。该程序还允许用户计算自由空气重力异常，相对于GRS 80和GRS 67参考椭球，从确定的重力值，并计算地形校正在每个调查网站使用棱镜公式和用户提供的数字高程模型。本文评述了将相对重力仪观测值换算成重力值的数学框架然后，它继续详细说明如何使用软件实现处理步骤版权所有©2018作者.由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据现行守则第3版用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-16-00075法律代码许可证GNU使用的代码版本系统Gsolve使用Python 2.7编写，GUI使用PyQt4。Gdal，地理空间数据抽象库是用Python、C和C++编写的。所有使用的库都是在开源软件许可证（如GNU）编译要求，操作环境依赖性我们建议安装最新的Anaconda Python2.7 64位发行版，因为它需要注意安装依赖项，除了gdal，然后可以通过命令提示符使用'否则Gsolve使用以下外部 python 库 ： PyQt4 ， Tkinter ， tkfileditram ， numpy ， xlrd ， xlwt ， xlutils ，matplotlib，ctypes和gdal。GUI是为Windows设置的如果可用，链接到开发人员文档/手册问题支持电子邮件jack. ga.gov.au通讯作者。电子邮件地址：Jack. ga.gov.au（J. McCubbine）。https://doi.org/10.1016/j.softx.2018.04.0031. 介绍地球Reilly [1]2352-7110/©2018作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx*130J. McCubbine等人/SoftwareX 7（2018）129G阿夫林b⏐ ⏐⏐ ⏐⏐ˆˆ• − −˜ ˜ˆ•×−•×•−ˆ−•×•×−•ˆ×将相对重力测量值转换为重力值的鲁棒最小二乘公式该公式依赖于相对测量，包括在绝对重力基准点的至少一次测量，并考虑了潮汐重力效应和重力仪漂移。最小二乘公式以前已实现方法1：正态最小二乘法。线性方程组由Eqs. （2）和（3）可以写成由等式（1）给出的矩阵形式（四）、CUP−1-Tx−zWoodward和Carman [2]在一个名为DSOLVE的FORTRAN计算机程序中。我们首先回顾Reilly重力测量观测值与重力值之间的关系。然后我们给两个Q0 0a=其中，H 中文（简体）不同地将相对重力观测方程约束到现有重力值的替代方法。对于这些公式的实际实施，我们已经开发了一个新的Python计算机程序（Gsolve）与图形用户界面。它进一步解释了赖利[ 1 ]原始方法中没有包括的相对重力测量和重力仪动力学的额外实际例如，观测数据的不确定性和重力参考值的不同不确定性水平Gsolve首先使用用户提供的校准表和/或根据观测结果确定的重 力 仪 校 准 因 子 将 相 对 重 力 仪 现 场 观 测 值 转 换 为 mGal（105m/s2）然后使用朗曼公式[3]校正潮汐重力效应的观测值，最后通过最小二乘法求解观测点的重力值，同时考虑重力仪漂移。详细给出了执行这些处理步骤的程序的操作2. 将相对重力观测值转换为重力值的方法框架Reilly [1]给出了在位置j和时间ti的相对重力仪读数xi，j，xi，j=gj−a−b（ti−t1）+zi，j−ei（1）其中，gj是位置j处的真实重力值，b（t1，t1）是对应于计量器已经漂移了时间t1的值，z i，j 是位置j和时间的潮汐引力效应ti（这是一个可以通过朗曼公式计算的已知量），a是重力仪零点（即在没有任何漂移或潮汐效应的情况下，相对重力仪刻度盘读数为零时的绝对重力值），以及Ei是与时间t1处的仪表读数相关联的（小）误差项。收集已知量和右侧的误差项，给出以下形式的观测方程g j− a − b（t i− t1）= x i，j − z i，j − e i.（二）在所有使用相对重力仪的测量中，存在M个位置处的N个观测值的子集，其中存在先前记录的重力值gj。相对于真实重力值gj，值gj还必须包括（小）误差项j。这导致第二个方程的形式，gj=gj+j。（三）为了将相对重力观测值与绝对重力值联系起来，必须至少有一个这种类型的方程即必须在至少一个已知重力的位置进行相对重力观测。线性方程组由Eqs.（2）和（3）可以用来确定最佳拟合参数gj，a和b通过最小广场。g是gjP是一个N M矩阵，对于位置j和时间t i的观测，j列和i行为1，其他位置为01是N1个向量T是N1向量的时间t it1，Q是M M矩阵，对于对应于等式（1）的位置，在行j的对角线上具有1。（3）其他地方都是0x z是每次观测的仪表读数减去潮汐效应的向量，h是具有元素gi的M1向量，其中矩阵Q的行中具有1，否则为零。然后可以寻求最小二乘解，其最小化向量e和ε中的平方误差项的和。这将给出绝对重力的最佳拟合参数估计，测点g、重力仪基线a和漂移系数b.该方法不将参考点与观测方程组解耦，并假设误差项不为零。这使得参考地点绝对重力值的最小二乘解与先前记录的值不同。这在参考绝对重力值不确定时是有利的，例如，如果这些值是有日期的或者是用不准确的仪器测量的。方法2：解耦最小二乘法，遵循Reilly [1]。Reilly（2）和（3）通过将参考重力值与观测方程解耦解耦的参考重力值是通过修改方程。（二）、对于在参考重力站点进行的观测，从观测方程的右侧减去参考重力值gj，并且丢弃该特定位置的值gj这给出了一个新的观测方程的形式，-a−b（ti−t1）=xi，j−zi，j−g<$j−ei（5）（4）这相当于对P和x z。将P的列j中的值g j设置为等于零，并且从向量x z中减去值gj，使得观测方程具有（5）的形式。该方法固定参考重力值，使得它们在最小二乘解中不会改变，从而迫使残差项精确等于零。然而，参考地点的估计重力值参数方差是非零的，尽管解精确地拟合到给定的先前记录的重力值。相对于所述的其他两种方法，这相当于参考重力值的中等置信水平。方法3：约束最小二乘法。这种方法解决了一组方程的形式方程。（1）具有类似于等式（1）的无误差形式的约束方程。（2）采用约束最小二乘法。最小二乘矩阵方程的公式如下。考虑以下等式Xβ=Y−e（6）·····J. McCubbine等人/SoftwareX 7（2018）129131Gg一一××+⏐ ⏐ ⏐ ⏐ ⏐⏐⏐ˆ˜b˜bλ布拉奇⏐ ⏐ ⏐ ⏐其中，X= [P，−1，−T]，β=，Y=[x−z]，其中P，T，x，z，并且e如方法1中所述。约束由矩阵方程给出，Kβ= c（7）其中K具有维数k（M2）对于某个k.这里k是参考重力站点的数量，因此K具有与Q相同的形式，但删除了零行。c是一个k1向量，并且包含在矩阵K中对应于其中具有1的列的位置的每行中的指定重力值。该方程类似于方程的无误差形式。（三）、与方法2类似，位置j处的参考重力值是固定的。这里的目标是最小化，S=（Y−Xβ）′（Y−Xβ）+λ′（Kβ−c）2）（8）其中λ是拉格朗日乘子的k×1向量S具有最小值，其中（1）≤S =0和（2）S =0。这是给3.1. 相对重力数据处理收集调查数据。相对重力测量、位置标识、时间和日期（UTC）以及每次观测的位置必须首先导入或手动输入程序。这些输入数据将用于构造线性方程组图1显示了执行此操作的图形用户界面;这是用户运行程序时显示的第一个页面用户可以从Excel电子表格中导入数据（按（1）），并使用示例测量数据文件example_survey_data.xls提供的指定格式，也可以手动逐个输入可以使用按钮（4）删除数据的列表，并且可以通过单击（2）将测量数据从Gsolve导出为电子表格。利用相对重力仪进行的重力测量通常被视为一系列每日观测（或“循环”），其中包括每天在第一个站点进行的一次收集通过，∂β ∂λ通过这种方式进行的原始观测，可以在不同的时间内确定重力仪漂移率和重力仪零点2X′（Y−Xλ）−K′λ=0（9）或（对于λ=λ/2）X′（Y−Xβ）−K′λ=0（10）这意味着，X′Xβ−K′λ=X′Y（11）和Kβ= c（12）这可以写成矩阵形式如下，期这有助于解释数据中的“零”，其中重力仪的零点由于撞击和颠簸而用户可以为观测数据的各个部分分配单独的循环编号。当求解线性方程组时，由方程组给出。（2）和（3），Gsolve计算每个回路的单独重力仪零点值a和漂移系数b。这是有利的，因为它可以用于解释重力仪零点值的偏差和数据的分离日之间的可变漂移率。绝对重力参考点和重力仪校准表存储在Gsolve中。Gsolve存储2个数据库，第一个包含用于Eq.（3）第二X′XK ′βX ′YK0（十三）包含制造商提供的重力仪校准表，用于将重力仪刻度盘读数转换为重力单位x i，j.图图2显示了可以使用的图形用户界面⏐β˜⏐<$X′XK′<$−1<$X′Y<$与数据库进行交互。新型重力仪校准表可以通过单击按钮（5）并按照屏幕上的说明（文件格式）误差项α被认为是零，方差也为零。当参考重力值被认为没有误差时，这种方法是最合适的，并且最适合参考重力值的高置信度本文给出了三种将相对重力观测值转换为重力值的方法尽管理论上存在差异，但一般来说，每种方法都应返回几乎相同的重力值。然而，由于参考重力值的精度的理论考虑，所导出的重力值方差3. 格索尔Gsolve是一个Python计算机程序，旨在通过使用直观的图形用户界面，使相对重力测量数据的常规处理用户友好。该程序将相对重力观测值转换为重力值，并可进一步将导出的重力值转换为重力异常。由软件执行的相对重力转换过程校正重力仪读数的漂移、潮汐重力效应和重力仪零点的跳跃（以下称为“零点”）以及附加重力仪校准因子。此外，重力异常减少工具可以计算自由空气重力异常（相对于GRS 80和GRS 67参考椭球），并且可以利用用户提供的数字高程模型计算重力地形校正该计划的细节在以下小节中给出。由calibration_table_example. xls给出），数据库可以通过从中选择一个下拉选项来查看(6)并且可以通过单击（7）从数据库中删除表。通过单击“导入”（8）并导入相同格式的文件，可以将新的绝对重力值添加到数据库中 作为example_absolute_gravity.xls，或单击建立线性方程组。用最小二乘法求解线性方程组的图形用户界面如图所示。3.第三章。首先，必须使用下拉菜单（12）指定接下来，用户必须指定参考绝对重力场，以构建Eq.（3）点击（13）。该程序将提示用户通过将调查的站点名称与存储的站点名称进行匹配来在数据库中查找它们，或者手动输入它们，任何错误添加的参考站点都可以通过选择该行并单击（14）来删除。然后，用户必须选择是否为中确定的每个环计算单独的重力仪零点a和漂移率b。点击（15）查看调查数据通常，用多个回路处理数据会导致较小的gj参数标准误差估计值，因为单独回路之间的重力仪零点中的误差会得到更好的解释。软件允许用户选择残差置信区间。这可以通过编辑（16）中的值来改变，其中最佳拟合参数由下式找到.（132J. McCubbine等人/SoftwareX 7（2018）129˜ ˜G一⏐⏐⏐ ⏐βFig. 1. 图形用户界面，将相对重力测量数据导入Gsolve。图二. （a）显示了图形用户界面，用于从Gsolve数据库中查看、导入和删除重力仪校准表。(b)显示了用于从Gsolve数据库中查看、导入、修改、保存和删除绝对重力参考站点的图形用户界面。默认值为100%。当此值小于100%时，Gsolve将生成初始解，识别指定置信区间之外的离群值，然后计算删除离群值的新解这对于消除错误测量确定的参数gj、aj和bj。附加重力仪校准因子。Reilly [1]确定了一种方法来确定相对重力仪读数的进一步校准比例因子β。这在重力仪校准表不可用或不精确的情况下很有用为了通过最小二乘法计算该值，Eq.（4）必须修改为Eq。（十五）、指定的校准表（和β因子）。然后计算潮汐重力效应，并使用Longman [3]公式，然后整理线性方程组，求出参数g、a、b、β的最小二乘解决心已定表中列出了每个站点的最小二乘绝对重力测定值以及测定值的标准误差和所用观测值的数量通过点击（21）并按照屏幕上的说明，可以将表中的数据P−1−T xx-z⏐e⏐详细输出。 最小二乘解的详细输出为2000年0月0日0b=h+（十五）在“详细结果和绘图”选项卡中给出这在图4中示出。对于测量中的每个观测，站点名称、位置、刻度盘和校准刻度盘读数、日期和时间、潮汐效应、回路为了确定β值，必须参考测量数据至少有两个绝对重力场β值可以通过点击（18）由Gsolve确定。或者，它可以由用户在（17）中使用先前导出的值来指定，并且用于在确定最小二乘解之前校准重力仪读数。默认选项是不计算β并将其设置为零，即它被排除在计算之外解线性方程组。然后，用户可以指定要使用的最小二乘法（即，正常最小二乘法、解耦最小二乘法或约束最小二乘法），以找到每个测量位置处的绝对重力。最小二乘计算最后通过按下“求解”（20）启动重力仪读数首先使用同一性和残差在表中给出每个回路的每小时漂移率由（22）给出，β值（无论是指定值还是计算值）由（23）给出。任何单个环路的漂移函数和漂移函数周围的观测残差可以使用（25）绘制，类似地，观测残差的累积分布函数可以使用（24）绘制这些图有助于识别离群值或使用置信区间是否合适（即，是否出现来自漂移函数的大的异常值或累积分布图中的长尾详细的输出数据可以通过点击导出到xls文件（26）。J. McCubbine等人/SoftwareX 7（2018）129133−√图三. 指定要使用的重力仪校准表、绝对重力参考点、最小二乘解的设置并计算绝对重力值。3.2. 计算重力异常自由空气异常一旦在每个观测点获得绝对重力值，Gsolve也可用于确定自由空气重力异常。自由空气异常ΔgFA定义为观测重力g与在观测高度和纬度处估算的地球近似椭球体重力之间的差。即gFA=g−γ（φ）−δgFA（ 16）其中，γ（φ）是椭球重力，δg FA是使用线性公式0估算的自由空气效应。3086HmGal/m，其中H为观测点高度，单位为m。γ（φ）可根据GRS 80或GRS 67参考椭球参数计算。这通过从下拉菜单中选择相应的椭圆体并点击更新按钮（27）来完成。用于查看和导出（使用按钮（28））重力校正值和重力异常的图形用户界面在图2B中示出。五、椭球重力是使用Moritz [4]公式计算的，公式如下：（十七）、1+ksin2（φ）和GRS 67参考椭球体为• γa= 978031。84558（mGal）• k= 0。001931663383• e2 = 0。00669460532856地形校正。布格重力异常是指对观测点周围地形的重力影响进行校正后的自由空气重力异常地形的重力影响通常在两个部分中进行校正，布格板δgbs和地形校正δgTC，因此，BA=BA−δgbs+ δgTC。（ 十八）Bouguer平板校正δbs（单位：mGal）由方程（1）给出（十九）、δgbs=−0。0419ρH（ 19）式中，H为观测高度，单位为米，ρ为岩石密度，单位为g/cm3。通常使用2.67 g/cm3的值Gsolve可用于从数字高程模型（DEM）确定地形校正δgTC该过程使用Naggy [5]（）=22（十七）棱镜公式用户可以读取数字高程模型（DEM）γφγa1−e sin（φ）其中φ是椭球纬度。用于GRS 80椭球体的参数是，• γa= 978032。67715（mGal）• k= 0。001931851353• e2 = 0。00669438002290Python库GDAL支持的50多种格式中的任何一种（2016），这里使用简单的xyz格式作为示例，并通过文件“NLDEM.txt”作为示例数据集导入DEM的用户界面如图所示。 六、按钮(29) 允许用户选择他们的DEM文件，格式被自动破译，然后绘制，分辨率在文本框（30）中指定为了继续计算，用户必须134J. McCubbine等人/SoftwareX 7（2018）129见图4。 详细的结果输出图形用户界面。图五. 测点重力异常推导的图形用户界面。在这里可以选择GRS80或GRS67参考椭球体，数据可以导出到xls文件。J. McCubbine等人/SoftwareX 7（2018）129135见图6。 用于导入数字高程模型文件的图形用户界面。首先导入位置数据文件。其中包含重力场的东向和北向。图形用户界面这样做是显示在图。7.第一次会议。此输入文件与读入的相对重力数据处理文件类型相同，但在xls文件的位置表上必须有两个附加通道，即“东距”和“北距”。东距和北距必须以米为单位，并且与DEM的投影相同。这些附加列显示在示例文件“example_survey_data.xls”中当读取重力观测点时（通过点击（31）），程序提示用户指定用于校正的重力观测点的高程类型。提示提供了三个选项，1. Gsolve可以从DEM中确定高度。这可以通过单击2. 用户可以将测量点的DEM高度保存在测量数据输入文件中，Gsolve将读取它们（它们存储在输入数据位置表的第7列中）。这是通过对第一个选项单击“否"并对第二个选项单击”是“来执行的3. 或者，用户可以将现场观测高度（例如GPS）保存在测量输入数据文件中，Gsolve将读取它们（它们存储在位置表的第4列中）。这是通过对第一个和第二个选项单击“否"来执行的如果重力场位于DEM边界内，则通常使用选项1。如果重力场位于DEM边界之外，但之前已计算DEM高度，则可将其放置在勘测数据文件的第7example_survey_data.xls）的2应该使用。在某些情况下，使用者可能会发现最好使用在野外建立的高度或某些其他高度进行地形校正，在这种情况下，可以使用选项3。然后，用户必须选择 校正仅考虑具有指定内半径和外半径的环形物内的地形）和岩石密度（34）（g/cm 3）（使用图1所示的图形用户界面）。8.第八条。为了启动计算，用户需要点击按钮（35）。当计算完成后，图上的“x”标记变为白色，并且地形校正计算的最终结果在图形用户界面的“结果”选项卡下给出结果表包含重力场名称、位置数据和地形校正。还给出了计算的地形校正的参数，并且可以使用按钮（36）导出数据该软件使用户能够分几步计算地形校正。例如，如果要在陆上和海上进行地形校正，用户可以首先使用密度值为2.67 g/cm3的陆上DEM数据，然后使用密度值为1.67 g/cm3的测深数据。或者，类似地，用户可能发现使用不同分辨率的几个DEM来加速计算是有利的例如，可以用2m分辨率DEM计算从重力场向外到170 m半径的内部地形校正，可以用20 m DEM计算从170 m向外到2160 m的中间地形校正，并且可以用20 m DEM计算向外到21900 m或更远的外部校正。136J. McCubbine等人/SoftwareX 7（2018）129见图7。 图形用户界面导入重力场数据进行地形校正。图8.第八条。 图形用户界面设置参数和计算地形校正。然后，可以对各个地形校正进行求和，以确定每个重力场的总校正。4. 结论本文评述了Reilly此外，两个替代的最小二乘配方的问题已经给出。每种方法都将相对重力测量值与参考重力值联系起来，但精确度不同介绍了一个新的Python计算机程序Gsolve该计划的目的是协助例行用最小二乘法将相对重力测量数据归约为重力值。它有一个直观的图形用户界面，为用户unfamiliar与命令行界面，使它可以更容易，比以前的基于FORTRAN的程序更容易采用。该程序采用的处理链考虑了太阳和月球的潮汐重力效应、重力仪漂移和测量数据中的重力误差。当使用两个或多个参考绝对重力站点时，它还确定重力仪的进一步校准因子。除了将相对重力观测值转换为重力值外，还能够计算相对于GRS 80和GRS 67参考椭球体的J. McCubbine等人/SoftwareX 7（2018）129137图9.第九条。 图形用户界面设置参数和计算地形校正。并根据用户提供的数字高程模型和指定的密度参数计算地形校正这为用户提供了所有必要的工具，以稳健和一致的方式计算布格重力引用[1] Reilly WI重力仪观测的调整。新西兰J GeoelGeophys1970;73：697-702.[2] 伍德沃德DJ，卡门AF。减少精密重力观测的计算机程序，地球物理部，技术说明第93号，1984年。[3] 朗曼IM.太阳和月亮引起的潮汐加速度的计算公式。《地球物理研究杂志》，1959年;64：2351-5。[4] 莫里茨·H 1980年大地测量参考系统。Bull Géodésique 1980;54（3）：395http://dx.doi.org/10.1007/BF02521480网站。[5] 纳吉·D直角棱镜的万有引力。 Geophysics1966;31（2）：362-71.