□□ 陳俊祺,上官凱林,劉會芳 (武警山西總隊參謀部,山西 太原 030012)

引言

變形監(jiān)測外業(yè)工作時會遇到各種復(fù)雜的情況。由于被監(jiān)測地物不斷變形的特殊性質(zhì),在其附近較難找到穩(wěn)定的工作點。要在被監(jiān)測地物附近建立基準(zhǔn)點和工作點組成首級網(wǎng)[1],然而這個基準(zhǔn)點可能和監(jiān)測物一樣存在變形,且可靠性尚未可知。最終數(shù)據(jù)是點位高程而非觀測值高差,可將點位高程設(shè)置為參數(shù),采用間接平差或者秩虧法平差法進(jìn)行內(nèi)業(yè)工作。但秩虧自由網(wǎng)法和間接平差各有利弊。在只有一期觀測值時,且基準(zhǔn)任意的自由網(wǎng),采用秩虧自由網(wǎng)法比間接平差法工作量大,且在一定條件下二者平差結(jié)果相同。但當(dāng)具有類似變形監(jiān)測多期觀測值時,仍然采用間接平差法,結(jié)果受基準(zhǔn)點影響較大,可采用秩虧法進(jìn)行平差[2-4]。

現(xiàn)實生活中的外業(yè)工作得到的數(shù)據(jù)往往是龐大的,目前多采用誤差平均分配法進(jìn)行內(nèi)業(yè)工作。其主要原因是由于采用加權(quán)法平差工作量過大。為了提高內(nèi)業(yè)工作質(zhì)量,迫切需要一種可以處理龐大數(shù)據(jù)的秩虧平差軟件。而目前關(guān)于秩虧水準(zhǔn)網(wǎng)的平差軟件少之又少,大多為算法理論研究。根據(jù)文獻(xiàn)[5]和文獻(xiàn)[6]可知,秩虧法平差大致可以分三類:一類是利用廣義逆矩陣?yán)碚摰腗ittermayer偽逆解法;第二類是Mittermayer附加條件法、Pelzer偽觀測值法和Perelmuter消去條件法;第三類是通過一定技術(shù)將將秩虧自由網(wǎng)平差轉(zhuǎn)化為間接平差的Wolf直接解法。下文將基于第二類附加條件法進(jìn)行軟件的開發(fā)。

1 秩虧自由網(wǎng)平差原理以及數(shù)學(xué)模型

1.1 間接平差

設(shè)t為必要觀測值,n為觀測值數(shù)目。間接平差的函數(shù)模型見式(1):

(1)

(2)

將式(2)代入式(1)見式(3):

l=L-(BX0+d)=L-L0

(3)

式中L0為觀測值的近似值矩陣,l為誤差方程常數(shù)項,也就是觀測值與其近似值之差,因而可以推導(dǎo)出誤差方程矩陣見式(4):

(4)

式中:V——觀測值改正數(shù)。

假設(shè)權(quán)陣為P,則平差的準(zhǔn)則見式(5):

VTPV=min

(5)

(6)

將最后一項轉(zhuǎn)置后見式(7):

BTPV=0

(7)

將式(4)帶入式(7)見式(8):

(8)

(9)

式(9)中系數(shù)矩陣NBB為滿秩的非奇異矩陣,可以求得唯一解,這就是間接平差函數(shù)模型求解方法[7]。

1.2 秩虧自由網(wǎng)平差法與附加條件法

附加條件法的基本思路是:用間接平差法平差自由網(wǎng)時,由于網(wǎng)中沒有足夠的起算數(shù)據(jù),所以法方程系數(shù)矩陣NBB秩虧。設(shè)秩虧數(shù)為m,添加m個參數(shù)條件的限制方程,可按照附有限制條件的間接平差解算方法進(jìn)行解算。條件方程見式(10):

(10)

式中:S——參數(shù)限制條件方程系數(shù)矩陣;

x——參數(shù)改正值;

u——參數(shù)個數(shù)。

限制條件方程應(yīng)與法方程線性無關(guān),這一要求等價于關(guān)系式(11):

(11)

限制條件方程中的方程式也與法方程線性無關(guān)。所以S矩陣的秩為m。利用拉格朗日乘數(shù)法求偏導(dǎo)數(shù),然后得到法方程(通過聯(lián)立參數(shù)限制條件得到法方程組),最終得到計算秩虧自由網(wǎng)最優(yōu)解,見式(12):

(12)

其單位權(quán)中誤差見式(13):

(13)

式(13)中根號下分母應(yīng)不為0,否則使用軟件計算時容易出錯,即出現(xiàn)單位權(quán)中誤差無窮大的情況,從而導(dǎo)致在Matlab編程時協(xié)因數(shù)陣元素全部為INF(在Matlab中INF為無窮大)。對于秩虧水準(zhǔn)網(wǎng)來說,秩虧數(shù)為1,根據(jù)平差后所有高程點的改正值之和為0,列出條件方程見式(14):

(14)

因而對于秩虧水準(zhǔn)網(wǎng)來說,S的具體形式見式(15)[7]:

(15)

2 程序需求分析

隨著我國社會的不斷發(fā)展,高層建筑不斷出現(xiàn),民用設(shè)施不斷完善更新,人們對于開采業(yè)安全性的需求不斷提升,變形監(jiān)測的應(yīng)用也越來越廣泛。然而變形監(jiān)測的數(shù)據(jù)是多期觀測的大量數(shù)據(jù),人工解算效率低下且易出錯,已經(jīng)不適應(yīng)當(dāng)代發(fā)展的需求。外業(yè)一般得到的原始數(shù)據(jù)有起點和終點點號、觀測數(shù)據(jù)、路線長度等。根據(jù)式(12),計算改正值所需要的矩陣分別為誤差方程系數(shù)矩陣B、權(quán)陣P、附加條件系數(shù)矩陣S以及誤差方程常數(shù)項l等四項,所以代碼編寫的主要任務(wù)是從外業(yè)的原始數(shù)據(jù)得到解算參數(shù)改正值的所需矩陣。

3 程序設(shè)計與集成

3.1 程序結(jié)構(gòu)

軟件主要結(jié)構(gòu)為:數(shù)據(jù)讀取與導(dǎo)入模塊、數(shù)據(jù)分類以及調(diào)用模塊、數(shù)據(jù)運算模塊以及保存和清除模塊,軟件主要實現(xiàn)步驟如圖1所示。

圖1 軟件實現(xiàn)步驟

3.2 數(shù)據(jù)結(jié)構(gòu)

在Matlab訪問以及將Excel數(shù)據(jù)表格的點號、路線長度以及觀測高差放置在原數(shù)據(jù)矩陣中,以原數(shù)據(jù)矩陣的列向量為單元,對數(shù)據(jù)進(jìn)行分類,分別對起點、終點、觀測高差以及路線長度以不同的變量命名,以方便調(diào)用。

3.3 關(guān)鍵算法

3.3.1提取誤差方程系數(shù)矩陣

誤差方程系數(shù)矩陣B中元素只有0、-1和1,而且矩陣的每一行有且僅有一個“1”和一個“-1”,其他元素均為0。當(dāng)每一行的列號所對應(yīng)的點號是終點時,這個元素應(yīng)為“1”,當(dāng)為起點時,這個元素應(yīng)為“-1”,否則為“0”。這樣就可以由起點和終點提取出誤差方程系數(shù)矩陣B。

先定義系數(shù)矩陣B為一個零矩陣。然后找出每行起點和終點點號對應(yīng)的列號,并分別賦值-1和1,檢索原數(shù)據(jù)矩陣的起點和終點點號。假設(shè):I為原始數(shù)據(jù)矩陣,I[1,1]代表矩陣I的第一行第一列,n為數(shù)據(jù)行數(shù),系數(shù)矩陣流程如圖2所示。

圖2 系數(shù)矩陣提取流程

3.3.2提取權(quán)陣

水準(zhǔn)路線定權(quán)和路線長度有關(guān),可依據(jù)此來定權(quán)。

3.3.3提取附加條件系數(shù)矩陣

由秩虧水準(zhǔn)網(wǎng)平差參數(shù)改正值之和為0,可以推導(dǎo)出附加條件系數(shù)矩陣為式(7),直接定義即可。

3.3.4提取誤差方程常數(shù)項

式(3)中B已經(jīng)提取出,X0是各個點的參數(shù)估計值,如果給出,則可以直接計算。當(dāng)導(dǎo)入數(shù)據(jù)后,發(fā)現(xiàn)未給出參數(shù)估計值,可以通過給出假定起算點高程值來進(jìn)行參數(shù)估計值賦值,假設(shè)參數(shù)估計值矩陣為X。

參數(shù)自動賦值編程是假定起點1的高程值為0,并定義一個n行1列的矩陣,然后從起點到終點第一次逐行遍歷,只要起點為1,就可以根據(jù)高差對終點對應(yīng)的點號賦值;再從終點到起點第二次逐行遍歷,只要終點為1,就可以根據(jù)高差對起點對應(yīng)的點號賦值。這兩次遍歷方法是一樣的,稱為第一類遍歷。第一類遍歷流程如圖3所示。

圖3 第一類遍歷流程

這樣完成了賦值的初步工作。由于并不是所有的點都和點號為1的點有關(guān),所以還有一些點仍然未進(jìn)行賦值,所以需要下一步工作。從起點到終點進(jìn)行第三遍遍歷。依次判斷每一行的起點、終點和高差進(jìn)行判斷,判斷條件為:起點參數(shù)估計值不為0;起點點號不為1;終點參數(shù)估計值為0。若都滿足,根據(jù)兩點高差和起點參數(shù)估計值對終點參數(shù)估計值進(jìn)行賦值。這樣就完成了第三遍遍歷。

在經(jīng)歷多組數(shù)據(jù)測試以后。發(fā)現(xiàn)第三遍遍歷在水準(zhǔn)網(wǎng)較為復(fù)雜時會出錯,具體原因是正向遍歷并不能完成所有點的賦值,還需要進(jìn)行一遍反向遍歷,可將所有點的參數(shù)估計值都完成賦值,所以需要進(jìn)行第四次反向遍歷,判斷條件為:終點參數(shù)估計值不為0;終點點號不為1;起點參數(shù)估計值為0。若都滿足,可以根據(jù)觀測值和終點點號對應(yīng)的參數(shù)估計值對起點進(jìn)行賦值。第三次和第四次遍歷是相似的,稱為第二類遍歷。這樣就完成了賦值工作。當(dāng)然,程序還有許多需要正反遍歷進(jìn)行賦值以及計算工作,方法和X0的賦值基本類似。第二類遍歷流程如圖4所示。

圖4 第二類遍歷流程

4 程序測試

在GUI界面和運算代碼編寫完成之后,利用Matlab自帶的deploytool工具將運算代碼和界面代碼等打包為可以脫離Matlab獨立運行的程序。保存目錄中會有一個exe文件,就是設(shè)計好的軟件界面啟動文件,可以創(chuàng)建其快捷方式并放在桌面。這個保存目錄是已更改的名字,并不是原名,原名為for-test。打開之后可以進(jìn)行數(shù)據(jù)測試工作。

點擊導(dǎo)入數(shù)據(jù)按鈕,然后點擊計算,就可以得到各種需要的數(shù)據(jù),點擊保存可以得到結(jié)果表。為了驗證平差程序的正確性,以文獻(xiàn)[3]的數(shù)據(jù)來進(jìn)行驗算,見表1。水準(zhǔn)網(wǎng)示意如圖5所示,軟件運行平差結(jié)果如圖6所示,對比文獻(xiàn)[3]的平差結(jié)果,見表2。

表1 測試數(shù)據(jù)

表2 對比結(jié)果

圖5 水準(zhǔn)網(wǎng)示意圖

圖6 平差結(jié)果

由表2可知,1號點和4號點平差結(jié)果是有差異的,分別差了0.1 mm,但是總的平差改正值之和是不變的,依然為0。

5 結(jié)語

目前對于地物變形監(jiān)測的主要方法是基于水準(zhǔn)的多期觀測,數(shù)據(jù)必然是大量的,并且多采用誤差平均分配法進(jìn)行內(nèi)業(yè)處理,其主要原因是因為通過定權(quán)分配誤差計算工作量比較大。無論間接平差還是秩虧自由網(wǎng)平差,都需要進(jìn)行大量的參數(shù)估計值賦值工作。利用Matlab開發(fā)的程序進(jìn)行參數(shù)估計值賦值工作,極大地減少了變形監(jiān)測內(nèi)業(yè)計算的繁瑣程度。