王生林，李鳳廷，梁 晰，霍成勝，白國龍

（青海省第三地質(zhì)礦產(chǎn)勘查院，青海 西寧 810008）

0 前言

MATLAB［1－5］是 Matrix和 Laboratory兩個(gè)英文單詞的前三個(gè)字母的組合，是美國 Math－Works公司最早在二十世紀(jì)七十年代推出的產(chǎn)品，最初用FORTRAN語言編寫，主要功能是實(shí)現(xiàn)程序庫的接口功能。進(jìn)入九十年代以來，隨著軟件的不斷升級和更新，MATLAB發(fā)展成為國際公認(rèn)的標(biāo)準(zhǔn)計(jì)算軟件，在數(shù)值計(jì)算及可視化方面的功能不斷增強(qiáng)，成為當(dāng)今最優(yōu)秀的科技應(yīng)用軟件之一，具有強(qiáng)大的科學(xué)計(jì)算能力、可視化功能、開放式可擴(kuò)展環(huán)境，所附帶的工具箱支持三十多個(gè)領(lǐng)域的計(jì)算、仿真等應(yīng)用，因此在許多科學(xué)領(lǐng)域中，MATLAB成為計(jì)算機(jī)輔助設(shè)計(jì)和分析、算法研究及其應(yīng)用開發(fā)的基本工具和首選平臺。同時(shí)，MATLAB具有其它高級語言難以比擬的一些優(yōu)點(diǎn)即M文件編寫簡單、效率高、易學(xué)易懂，在信號處理、通信、自動(dòng)控制及科學(xué)計(jì)算等領(lǐng)域中被廣泛應(yīng)用，特別是在矩陣運(yùn)算、數(shù)組處理等方面具有很大的優(yōu)勢，是工程師、科研工作者最易上手的編程語言、最好的工具和環(huán)境。因此，MATLAB語言也被通俗地稱為演算紙式的科學(xué)算法語言。

重力勘探方法較為成熟，基點(diǎn)網(wǎng)平差計(jì)算也不例外，在2003年馮治漢［6］將以前單點(diǎn)平差的計(jì)算方法推導(dǎo)成矩陣的運(yùn)算，同時(shí)用MATLAB編程驗(yàn)證了《區(qū)域重力調(diào)查規(guī)范》中的例子，這給基點(diǎn)網(wǎng)平差計(jì)算帶來了一次跨越式的發(fā)展，重力基點(diǎn)網(wǎng)由筆算進(jìn)入了初步的計(jì)算機(jī)時(shí)代。在實(shí)際的應(yīng)用當(dāng)中，矩陣的運(yùn)算只有用MATLAB計(jì)算方便，但MATLAB在編譯生成獨(dú)立運(yùn)行文件時(shí)存在一定的局限性（MATLAB自身的函數(shù)大多數(shù)不能進(jìn)行編譯）。作者通過編寫的M文件，將指定文件中基點(diǎn)平差的數(shù)據(jù)讀入，再進(jìn)行計(jì)算得到準(zhǔn)確的結(jié)果。使用者只需將M文件和對應(yīng)的Excel文件放在同一路徑下，打開MATLAB編譯環(huán)境界面，改變路徑，輸入M文件名再回車，會提示使用者在Excle相應(yīng)的位置輸入相關(guān)參數(shù)，然后就可以進(jìn)行計(jì)算了。

1 重力基點(diǎn)網(wǎng)平差原理［7－14］

1.1 平差計(jì)算

重力勘探中基點(diǎn)網(wǎng)平差，多采用條件式平差，具體方法和精度的評價(jià)如下。

設(shè)基點(diǎn)網(wǎng)有n個(gè)邊段，m個(gè)未知基點(diǎn)，并由r（r＝n－m）個(gè)閉合圈組成。各圈的閉合差為W；各邊段的改正數(shù)為V；改正數(shù)條件方程系數(shù)為A；A′為A的轉(zhuǎn)置；各邊段的獨(dú)立增量數(shù)構(gòu)成對角矩陣P；L為各邊段的平差前的重力增量值；X為各邊段平差后的重力增量值；各基點(diǎn)在各邊段的方向矩陣為F；平差后各基點(diǎn)的重力值為G；聯(lián)系數(shù)為K，則改正數(shù)條件方程和聯(lián)系數(shù)法方程分別為：

聯(lián)立式（1）、式（2）得

1.2 精度評價(jià)

單位權(quán)中誤差μ為式（6）。

轉(zhuǎn)換系數(shù)Q滿足方程式（7）。

某基點(diǎn)平差值函數(shù)的權(quán)倒數(shù)G為式（8）。

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

某基點(diǎn)的平差值函數(shù)中誤差MG為

基點(diǎn)網(wǎng)的精度，即為整個(gè)網(wǎng)內(nèi)最弱點(diǎn)中的誤差MGmax。

2 程序的部分MATLAB源代碼

該程序采用函數(shù)sum＝xlsread（′文件名.xls′），讀取指定格式和文件名下的平差基本數(shù)據(jù)，再按照基點(diǎn)網(wǎng)平差原理進(jìn)行計(jì)算，得到最終結(jié)果。以下是該程序的部分程序源代碼，供同行參考。

% A系數(shù)矩陣；P權(quán)矩陣；L各邊增量矩陣；W閉合差矩陣；F各邊段的向量矩陣

% 各基點(diǎn)的重力值的均方誤差，選擇最大值作為本網(wǎng)的均方誤差

Mg1＝sqrt（miu＊pg）；

3 M文件的使用及應(yīng)用效果

3.1 M文件的使用

使用《區(qū)域重力調(diào)查規(guī)范》中的例子，如圖1所示。

圖1 某重力基點(diǎn)網(wǎng)分布示意圖Fig.1 Sketch map of gravitational base point net

校正數(shù)條件方程式為：

條件式系數(shù)為：

將編寫的M文件（名為JDWPC）和一個(gè)名為“jidianwangpingcha.xls”的文件放在同一目錄下，在MATLAB編譯環(huán)境中，輸入“JDWPC”再回車，就會提示在Excel中所要輸入基點(diǎn)網(wǎng)平差基本參數(shù)如圖2所示。

圖2 M文件執(zhí)行時(shí)界面圖Fig.2 Interface map of execution the m file

在Excel文件中輸入的數(shù)據(jù)格式如圖3所示。

編輯好Excel文件中的各項(xiàng)數(shù)據(jù)后保存，在運(yùn)行界面輸入“0＋Enter”就得到計(jì)算的結(jié)果（見圖4）。由于數(shù)值的四舍五入等原因，造成閉合差不為零，這時(shí)將不符值分在不與鄰環(huán)接界的權(quán)較小的邊上，本例中將V3舍去尾數(shù)取為“8”，V4不是舍去尾數(shù)取為“12”而是進(jìn)位取為“13”，通過運(yùn)行界面的提示重新輸入（見圖5）。繼續(xù)計(jì)算各個(gè)基點(diǎn)的相對已知點(diǎn)的重力差值，如果輸入已知基點(diǎn)的絕對重力值就可以計(jì)算出各個(gè)未知基點(diǎn)的絕對重力值，當(dāng)輸入“0”時(shí)，就顯示相應(yīng)基點(diǎn)相對于已知點(diǎn)的重力差，最后顯示出基點(diǎn)網(wǎng)各未知基點(diǎn)的平差精度（見圖6）。

3.2 應(yīng)用效果

通過以上M文件的運(yùn)行計(jì)算結(jié)果為：

四舍五入得到結(jié)果為：

閉合差不為零，通過人為的調(diào)整使得個(gè)閉合圈閉合差為零后，將最終的結(jié)果重新輸入為：

計(jì)算得到平差后各邊段的增量值。

各基點(diǎn)相對重力值的均方誤差：

取其最大均方誤差值為該基點(diǎn)網(wǎng)的精度為±24×10－8m／s2，計(jì)算結(jié)果與規(guī)范一致。

圖3 Excel文件中各項(xiàng)數(shù)據(jù)的格式示意圖Fig.3 Sketch map of data format in excel file

4 結(jié)論

MATLAB是一個(gè)很優(yōu)秀的計(jì)算編程軟件，對數(shù)據(jù)的處理及基本的繪圖有很好地利用前景。M文件的編寫簡單，對基本軟件Excel內(nèi)的數(shù)據(jù)讀寫簡便，深受科研工作者的青睞。

M文件的使用簡單方便，使用者只需安裝MATLAB編譯環(huán)境，將名為“JDWPC”的 M文件和一個(gè)名為“jidianwangpingcha.xls”文件都放在運(yùn)行界面的窗口目錄下，在Excel的sheet1下，輸入改正數(shù)條件方程式系數(shù)；在sheet2中的A列依次輸入各邊段的增量數(shù)（權(quán)數(shù)）P，B列輸入各邊段平差前的增量值L，C列輸入閉合圈的閉合差；在sheet3中輸入方向矩陣F，不需改變原程序中的任何內(nèi)容就可以進(jìn)行計(jì)算。特別是邊段較多的網(wǎng)，如果在原程序中輸入相應(yīng)的數(shù)字就非常麻煩而且也容易出錯(cuò)，對于不懂MATLAB語言的使用者來說，一旦源程序出錯(cuò)將無法改正，本原程序只需在基本工具Excel表格中進(jìn)行編輯，使用方便不容易出錯(cuò)，是重力基點(diǎn)網(wǎng)平差的較為人性化且好用的程序。

［1］周建興，豈興明，矯津毅，等.MATLAB從入門到精通［M］.北京：人民郵電出版社，2008.

［2］劉維.精通Matlab與c／c＋＋混合程序設(shè)計(jì)（第2版）［M］.北京：北京航空航天大學(xué)出版社2008.

［3］汪洋兵，馬玄龍.Excel在重力基點(diǎn)網(wǎng)平差中的應(yīng)用［J］.資源環(huán)境與工程，2010（6）：701－706.

［4］郭良輝，孟小紅，石磊.基于 Matlab的重力基點(diǎn)網(wǎng)平差實(shí)驗(yàn)教學(xué)法［J］.科技信息：科學(xué)教研，2008（18）：24.

［5］葉景艷，錢美平，周錫明，等.利用VB編程完成基點(diǎn)網(wǎng)聯(lián)測中的各項(xiàng)計(jì)算［J］.物探化探計(jì)算技術(shù)，2004（1）：71－77.

［6］馮治漢.MATLAB及其在重力基點(diǎn)網(wǎng)平差中的應(yīng)用［J］.物探化探計(jì)算技術(shù)，2003，25（4）：336－339.

［7］羅孝寬，郭紹雍.應(yīng)用地球物理教程—重力、磁法［M］.北京：地質(zhì)出版社，1991.

［8］中華人名共和國地質(zhì)行業(yè)標(biāo)準(zhǔn).大比例尺重力勘查規(guī)范DZ／T0171－1997［S］.北京：中國標(biāo)準(zhǔn)出版社，1997.

［9］同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系，工程數(shù)學(xué).線性代數(shù)［M］.北京：高等教育出版社，2003.

［10］劉天佑.應(yīng)用地球物理數(shù)據(jù)采集與處理［M］.武漢：中國地質(zhì)大學(xué)出版社，2004.

［11］王寶仁，程新文.一種簡易快速的重力基點(diǎn)網(wǎng)平差方法［J］.石油物探，1988（02）：91－99.

［12］朱松濤.重力基點(diǎn)網(wǎng)的廣義逆矩陣平差法［J］.物探與化探，1983（01）：26－36.

［13］朱松濤.水準(zhǔn)網(wǎng)（重力重點(diǎn)網(wǎng)）的廣義逆矩陣平差法［J］.長安大學(xué)學(xué)報(bào)：地球科學(xué)版，1982（02）：107－117.

［14］俞炯霞.用條件觀測平差法進(jìn)行重力基點(diǎn)網(wǎng)的平差［J］.物探化探電子計(jì)算技術(shù)，1982（1）：62－69.