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.