基于秩虧法處理地表形變監測數據平差的設計

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

引言

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

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

1 秩虧自由網平差原理以及數學模型

1.1 間接平差

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

(1)

(2)

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

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

(3)

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

(4)

式中:V——觀測值改正數。

假設權陣為P,則平差的準則見式(5):

VTPV=min

(5)

(6)

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

BTPV=0

(7)

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

(8)

(9)

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

1.2 秩虧自由網平差法與附加條件法

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

(10)

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

x——參數改正值;

u——參數個數。

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

(11)

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

(12)

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

(13)

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

(14)

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

(15)

2 程序需求分析

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

3 程序設計與集成

3.1 程序結構

軟件主要結構為:數據讀取與導入模塊、數據分類以及調用模塊、數據運算模塊以及保存和清除模塊,軟件主要實現步驟如圖1所示。

圖1 軟件實現步驟

3.2 數據結構

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

3.3 關鍵算法

3.3.1提取誤差方程系數矩陣

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

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

圖2 系數矩陣提取流程

3.3.2提取權陣

水準路線定權和路線長度有關,可依據此來定權。

3.3.3提取附加條件系數矩陣

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

3.3.4提取誤差方程常數項

式(3)中B已經提取出,X0是各個點的參數估計值,如果給出,則可以直接計算。當導入數據后,發現未給出參數估計值,可以通過給出假定起算點高程值來進行參數估計值賦值,假設參數估計值矩陣為X。

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

圖3 第一類遍歷流程

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

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

圖4 第二類遍歷流程

4 程序測試

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

點擊導入數據按鈕,然后點擊計算,就可以得到各種需要的數據,點擊保存可以得到結果表。為了驗證平差程序的正確性,以文獻[3]的數據來進行驗算,見表1。水準網示意如圖5所示,軟件運行平差結果如圖6所示,對比文獻[3]的平差結果,見表2。

表1 測試數據

表2 對比結果

圖5 水準網示意圖

圖6 平差結果

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

5 結語

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