MATLAB在重力基点网平差中的应用

王生林，李凤廷，梁 晰，霍成胜，白国龙

（青海省第三地质矿产勘查院，青海 西宁 810008）

0 前言

MATLAB［1－5］是 Matrix和 Laboratory两个英文单词的前三个字母的组合，是美国 Math－Works公司最早在二十世纪七十年代推出的产品，最初用FORTRAN语言编写，主要功能是实现程序库的接口功能。进入九十年代以来，随着软件的不断升级和更新，MATLAB发展成为国际公认的标准计算软件，在数值计算及可视化方面的功能不断增强，成为当今最优秀的科技应用软件之一，具有强大的科学计算能力、可视化功能、开放式可扩展环境，所附带的工具箱支持三十多个领域的计算、仿真等应用，因此在许多科学领域中，MATLAB成为计算机辅助设计和分析、算法研究及其应用开发的基本工具和首选平台。同时，MATLAB具有其它高级语言难以比拟的一些优点即M文件编写简单、效率高、易学易懂，在信号处理、通信、自动控制及科学计算等领域中被广泛应用，特别是在矩阵运算、数组处理等方面具有很大的优势，是工程师、科研工作者最易上手的编程语言、最好的工具和环境。因此，MATLAB语言也被通俗地称为演算纸式的科学算法语言。

重力勘探方法较为成熟，基点网平差计算也不例外，在2003年冯治汉［6］将以前单点平差的计算方法推导成矩阵的运算，同时用MATLAB编程验证了《区域重力调查规范》中的例子，这给基点网平差计算带来了一次跨越式的发展，重力基点网由笔算进入了初步的计算机时代。在实际的应用当中，矩阵的运算只有用MATLAB计算方便，但MATLAB在编译生成独立运行文件时存在一定的局限性（MATLAB自身的函数大多数不能进行编译）。作者通过编写的M文件，将指定文件中基点平差的数据读入，再进行计算得到准确的结果。使用者只需将M文件和对应的Excel文件放在同一路径下，打开MATLAB编译环境界面，改变路径，输入M文件名再回车，会提示使用者在Excle相应的位置输入相关参数，然后就可以进行计算了。

1 重力基点网平差原理［7－14］

1.1 平差计算

重力勘探中基点网平差，多采用条件式平差，具体方法和精度的评价如下。

设基点网有n个边段，m个未知基点，并由r（r＝n－m）个闭合圈组成。各圈的闭合差为W；各边段的改正数为V；改正数条件方程系数为A；A′为A的转置；各边段的独立增量数构成对角矩阵P；L为各边段的平差前的重力增量值；X为各边段平差后的重力增量值；各基点在各边段的方向矩阵为F；平差后各基点的重力值为G；联系数为K，则改正数条件方程和联系数法方程分别为：

联立式（1）、式（2）得

1.2 精度评价

单位权中误差μ为式（6）。

转换系数Q满足方程式（7）。

某基点平差值函数的权倒数G为式（8）。

将式（7）代入式（8）得到式（9）。

某基点的平差值函数中误差MG为

基点网的精度，即为整个网内最弱点中的误差MGmax。

2 程序的部分MATLAB源代码

该程序采用函数sum＝xlsread（′文件名.xls′），读取指定格式和文件名下的平差基本数据，再按照基点网平差原理进行计算，得到最终结果。以下是该程序的部分程序源代码，供同行参考。

% A系数矩阵；P权矩阵；L各边增量矩阵；W闭合差矩阵；F各边段的向量矩阵

% 各基点的重力值的均方误差，选择最大值作为本网的均方误差

Mg1＝sqrt（miu＊pg）；

3 M文件的使用及应用效果

3.1 M文件的使用

使用《区域重力调查规范》中的例子，如图1所示。

图1 某重力基点网分布示意图Fig.1 Sketch map of gravitational base point net

校正数条件方程式为：

条件式系数为：

将编写的M文件（名为JDWPC）和一个名为“jidianwangpingcha.xls”的文件放在同一目录下，在MATLAB编译环境中，输入“JDWPC”再回车，就会提示在Excel中所要输入基点网平差基本参数如图2所示。

图2 M文件执行时界面图Fig.2 Interface map of execution the m file

在Excel文件中输入的数据格式如图3所示。

编辑好Excel文件中的各项数据后保存，在运行界面输入“0＋Enter”就得到计算的结果（见图4）。由于数值的四舍五入等原因，造成闭合差不为零，这时将不符值分在不与邻环接界的权较小的边上，本例中将V3舍去尾数取为“8”，V4不是舍去尾数取为“12”而是进位取为“13”，通过运行界面的提示重新输入（见图5）。继续计算各个基点的相对已知点的重力差值，如果输入已知基点的绝对重力值就可以计算出各个未知基点的绝对重力值，当输入“0”时，就显示相应基点相对于已知点的重力差，最后显示出基点网各未知基点的平差精度（见图6）。

3.2 应用效果

通过以上M文件的运行计算结果为：

四舍五入得到结果为：

闭合差不为零，通过人为的调整使得个闭合圈闭合差为零后，将最终的结果重新输入为：

计算得到平差后各边段的增量值。

各基点相对重力值的均方误差：

取其最大均方误差值为该基点网的精度为±24×10－8m／s2，计算结果与规范一致。

图3 Excel文件中各项数据的格式示意图Fig.3 Sketch map of data format in excel file

4 结论

MATLAB是一个很优秀的计算编程软件，对数据的处理及基本的绘图有很好地利用前景。M文件的编写简单，对基本软件Excel内的数据读写简便，深受科研工作者的青睐。

M文件的使用简单方便，使用者只需安装MATLAB编译环境，将名为“JDWPC”的 M文件和一个名为“jidianwangpingcha.xls”文件都放在运行界面的窗口目录下，在Excel的sheet1下，输入改正数条件方程式系数；在sheet2中的A列依次输入各边段的增量数（权数）P，B列输入各边段平差前的增量值L，C列输入闭合圈的闭合差；在sheet3中输入方向矩阵F，不需改变原程序中的任何内容就可以进行计算。特别是边段较多的网，如果在原程序中输入相应的数字就非常麻烦而且也容易出错，对于不懂MATLAB语言的使用者来说，一旦源程序出错将无法改正，本原程序只需在基本工具Excel表格中进行编辑，使用方便不容易出错，是重力基点网平差的较为人性化且好用的程序。

［1］周建兴，岂兴明，矫津毅，等.MATLAB从入门到精通［M］.北京：人民邮电出版社，2008.

［2］刘维.精通Matlab与c／c＋＋混合程序设计（第2版）［M］.北京：北京航空航天大学出版社2008.

［3］汪洋兵，马玄龙.Excel在重力基点网平差中的应用［J］.资源环境与工程，2010（6）：701－706.

［4］郭良辉，孟小红，石磊.基于 Matlab的重力基点网平差实验教学法［J］.科技信息：科学教研，2008（18）：24.

［5］叶景艳，钱美平，周锡明，等.利用VB编程完成基点网联测中的各项计算［J］.物探化探计算技术，2004（1）：71－77.

［6］冯治汉.MATLAB及其在重力基点网平差中的应用［J］.物探化探计算技术，2003，25（4）：336－339.

［7］罗孝宽，郭绍雍.应用地球物理教程—重力、磁法［M］.北京：地质出版社，1991.

［8］中华人名共和国地质行业标准.大比例尺重力勘查规范DZ／T0171－1997［S］.北京：中国标准出版社，1997.

［9］同济大学应用数学系，工程数学.线性代数［M］.北京：高等教育出版社，2003.

［10］刘天佑.应用地球物理数据采集与处理［M］.武汉：中国地质大学出版社，2004.

［11］王宝仁，程新文.一种简易快速的重力基点网平差方法［J］.石油物探，1988（02）：91－99.

［12］朱松涛.重力基点网的广义逆矩阵平差法［J］.物探与化探，1983（01）：26－36.

［13］朱松涛.水准网（重力重点网）的广义逆矩阵平差法［J］.长安大学学报：地球科学版，1982（02）：107－117.

［14］俞炯霞.用条件观测平差法进行重力基点网的平差［J］.物探化探电子计算技术，1982（1）：62－69.