陳 帥,王 鵬

(天津市普迅電力信息技術(shù)有限公司,天津 300000)

0 引 言

目前,高等級水準(zhǔn)測量仍然是礦區(qū)進行開采沉陷觀測的主要方法,而水準(zhǔn)網(wǎng)的數(shù)據(jù)處理流程繁瑣且量大,大多數(shù)礦區(qū)當(dāng)前使用的水準(zhǔn)網(wǎng)數(shù)據(jù)處理程序或軟件都不同程度的存在不足,如:程序不能自動進行閉合條件的搜索和閉合差的計算,單純依靠人工計算,費時費力;一般的平差程序是以固定基準(zhǔn)模型設(shè)計的,不能根據(jù)礦區(qū)的實際情況及工程本身特點選擇合適的平差基準(zhǔn)模型來處理數(shù)據(jù)[1-2];不能進行抗差處理,消除觀測粗差的影響等。

VB6.0是一種可視化、面向?qū)ο蟮拈_發(fā)語言[3]，文章基于VB6.0的語言環(huán)境,進行了礦區(qū)水準(zhǔn)網(wǎng)綜合數(shù)據(jù)處理程序的編寫,實現(xiàn)了對不同平差基準(zhǔn)的水準(zhǔn)網(wǎng)進行最小二乘估計和抗差估計,以及自動進行閉合環(huán)的搜索和閉合差的計算等功能,完善了礦區(qū)水準(zhǔn)網(wǎng)數(shù)據(jù)處理流程,減少了作業(yè)人員的工作量,提高了工作效率。

1 數(shù)學(xué)模型

1.1 平差基準(zhǔn)模型

l=L-(AX0+d),

(1)

根據(jù)最小二乘原理可知:

(2)

單位權(quán)中誤差為[5]:

(3)

式中,P為觀測值的權(quán)陣。

R(A)=m

給定約束條件為

由最小二乘原理[7]可知:

(4)

式中,

Qp=(ATPA+PxSSTPx)-1

=(N+PxSSTPx)-1.

(5)

取S為

ST=(1 1 … 1).

當(dāng)采用秩虧基準(zhǔn)時

Px=E.

(6)

當(dāng)采用擬穩(wěn)基準(zhǔn),

(7)

式中:μ2>d,μ1+μ2=μ.

1.2 抗差模型

(8)

式中:Wi的計算采用IGG3函數(shù),以固定基準(zhǔn)為例,其抗差解為

(9)

當(dāng)參數(shù)兩次估值之差的絕對值中的最大值小于迭代精度時,迭代停止,單位權(quán)中誤差為

(10)

式中:n0為權(quán)因子,且等于0的觀測值個數(shù)。

2 程序設(shè)計

2.1 功能設(shè)計

程序可實現(xiàn)的設(shè)計功能有:閉合差計算、近似高程計算、平差計算等,以上功能分別基于固定基準(zhǔn)、秩虧基準(zhǔn)及擬穩(wěn)基準(zhǔn)，將采集的水準(zhǔn)數(shù)據(jù)按照特定格式整理好,選擇基準(zhǔn)類型,從導(dǎo)入模塊進行導(dǎo)入,一次進行操作,注意在平差前需進行抗差參數(shù)設(shè)置。

2.2 搜索閉合條件并計算閉合差

外業(yè)水準(zhǔn)測量結(jié)束后,為保證測量結(jié)果的質(zhì)量,一般都要進行閉合差的計算,是否符合相應(yīng)等級水準(zhǔn)測量要求，由于環(huán)線和附合路線的閉合差計算都涉及到搜索最短路徑的,因此首先進行最短路徑搜索程序的設(shè)計。

限于篇幅,以下只給出最短路徑搜索的關(guān)鍵步驟代碼:建立兩個數(shù)組ljd和S,分別記錄每個點的鄰接點號及其到目標(biāo)點的路線長度;設(shè)m點為目標(biāo)點點號。

For j=1 To tn

ljd(j)=-1

S (j)=1E+30

Next i

S(m)=0

ljd(m)=m

Dim p1%, p2%

p1=begin(j)

p2=end(j)

If S(p2)> S(p1)+S12 Then

ljd(p2)=p1

S(p2)=S(p1)+S12

ElseIf S(p1)> S(p2)+S12 Then

Ljd(p1)=p2

S(p1)=S(p2)+S12

將上述編碼編寫成過程FSP(),在進行閉合差搜索與最小獨立閉合環(huán)搜索編碼時調(diào)用此過程。

閉合差限差的計算公式如下:

(11)

式中,sigma0為給定的先驗單位權(quán)中誤差。

2.3 求解未知點的近似高程

在水準(zhǔn)網(wǎng)進行平差之前需求解未知點的近似高程值,程序設(shè)計思路如下:給所有未知點的近似高程絕對值很大的負數(shù),對網(wǎng)中所有觀測值進行搜索,當(dāng)觀測值出現(xiàn)起點和終點一個大于該負數(shù),一個小于該負數(shù),說明該觀測值的一端近似高程已知,另一端未知,再根據(jù)已知點高程值和高差推算出未知點的近似高程值,那么在之后的程序判別過程中該點就作為已知點來處理,以此類推,直到求解出所有未知點的近似高程時程序結(jié)束。

計算近似高程的程序代碼如下:

Dim Elevation()As Double '存放高程和高程近似值

Dim begin(), end()'編號數(shù)組,存儲的是點的序

Dim h()' 高差觀測值數(shù)組

For j=1 To su 'su為高程未知點個數(shù)

Elevation (sn+j)=-9999 'sn為高程已知點個數(shù)

Next j

n=0

Do While k < su

For j=1 To sh 'sh為觀測高差數(shù)

If Elevation (begin(j))>-9999And Elevation (end(j))<-9999Then '起始點已知,終點未知

Elevation (end(j))= Elevation (begin (j))+h(j)

n=n+1

ElseIf Elevation (end(j))>-9999 And Elevation (begin(j))<-9999 Then '起始點未知,終點已知

Elevation (begin(j))=Elevation end(j))-h(j)

n=n+1

End If

Next j

Loop

2.4 誤差方程的形成

程序運行時應(yīng)生成誤差方程的B陣、P陣和l陣，P陣和l陣的形成相對較為簡單，限于篇幅，下面給出B陣形成的代碼。對于固定基準(zhǔn)，代碼如下：

For j = 1 To sh

If begin(j)> sn Then

b(j, begin (j)- sn)= -1

End If

If end(j)> sn Then

b(j, end(j)- sn)= 1

End If

Next i

2.5 抗差迭代計算

抗差估迭代計算過程對兩次求解的參數(shù)進行比較，當(dāng)它們的之差的絕對值的最大值小于迭代條件的時候，迭代停止。下面給出固定基準(zhǔn)中迭代實現(xiàn)的代碼：

nn=1

Do While nn > eps

For i = 1 To su

b(j)= a(j, 1)

Next j

Call Ca_UintV(sh, sdd, d, dv, M)'利用中位數(shù)計算單位權(quán)方差因子M以及標(biāo)準(zhǔn)化殘差數(shù)組dv的計算

For j = 1 To sh

W(j)= Wj("IGG3", dv(j), n0, n1)'權(quán)因子計算函數(shù)

p(j, ij)= p1(j, j)* W(i)

Next j

Call zxec(e, r, s, d, l, q, PVV)'最小二乘求解過程

For j = 1 To su

b(j)= Abs(a(j, 1)- b(j))

Next j

nn = max(b)'max為編寫的最大值求解函數(shù)

Loop

3 算例分析

采用某地實測閉合水準(zhǔn)路線數(shù)據(jù)進行程序驗證，網(wǎng)中共有25個高差觀測值，B0為已知高程點，具體數(shù)據(jù)如表1所示。

表1 觀測數(shù)據(jù)和已知點高程

現(xiàn)用以上數(shù)據(jù)導(dǎo)入程序，以固定基準(zhǔn)為例進行如下操作：

1)閉合條件搜索和閉合差計算

環(huán)閉合差：B02-B03-B05-B04-B02,W=0.000 1(限差：0.006 2).

環(huán)閉合差：B04-B05-B06-W17-B04,W=0.001 1(限差：0.005 4).

環(huán)閉合差：W17-B06-J11-B0-W13-W17,W=0.001 1(限差：0.006 4).

環(huán)閉合差：B02-B03-B01-B02,W=0.004 5(限差：0.006 2).

環(huán)閉合差：B03-B09-B08-B03,W=-0.000 6(限差：0.005 3).

環(huán)閉合差：B05-B03-B09-J33-B05, W=0.000 8(限差：0.004 9).

環(huán)閉合差：B0-J11-J39-CP6-B0,W=-0.002 1(限差：0.006 9).

環(huán)閉合差：J33-B05-B06-J11-J33,W=0.001 1(限差：0.005 4).

環(huán)閉合差：J39-B09-J33-J39,W=0.001 1(限差：0.004 3).

環(huán)閉合差：B08-B03-B01-B08,W=0.002 6(限差：0.005 9).

路線閉合差計算:無附合路線。

2)固定基準(zhǔn)的最小二乘估計

對表1中的數(shù)據(jù)作固定基準(zhǔn)最小二乘估計，求得的高程平差值與高程中誤差如表2所示。

表2 固定基準(zhǔn)最小二乘估計

3)秩虧基準(zhǔn)的最小二乘估計

對表1中的數(shù)據(jù)作秩虧基準(zhǔn)最小二乘估計，求得的高程平差值與高程中誤差如表3所示。

表3 秩虧基準(zhǔn)最小二乘估計

由表2和表3可以看出，由固定基準(zhǔn)和秩虧基準(zhǔn)的最小二乘結(jié)算的高程平差值基本一致。

4)在B05-B04觀測高差中加入-0.04 m模擬粗差，J33-B09觀測高差中加入+0.02 m模擬粗差，對含粗差的數(shù)據(jù)分別進行最小二乘估計和抗差估計，將結(jié)果對比，如表4所示。

表4 最小二乘估計與抗差估計對比表

通過表4可以看出，在數(shù)據(jù)中加入粗差后進行最小二乘估計，其單位權(quán)中誤差顯著增大，為6.9 mm，高程平差值與表2中高程平差值相差較大。而同樣的數(shù)據(jù)進行抗差估計的結(jié)果與表2中高程平差值基本一致，從而驗證了抗差估計程序的準(zhǔn)確性與可靠性。

4 結(jié)束語

基于VB6.0語言環(huán)境設(shè)計了礦區(qū)水準(zhǔn)網(wǎng)綜合數(shù)據(jù)處理的程序，實現(xiàn)了對三種不同平差基準(zhǔn)的水準(zhǔn)網(wǎng)數(shù)據(jù)進行最小二乘估計和抗差估計，同時實現(xiàn)了程序自動進行閉合環(huán)的搜索和閉合差的計算等功能，并對程序?qū)崿F(xiàn)中的關(guān)鍵性問題進行了比較詳細的分析并給出了部分程序?qū)崿F(xiàn)的代碼。通過實際算例，分析驗證了程序的可靠性與準(zhǔn)確性，同時也彌補了現(xiàn)有礦區(qū)水準(zhǔn)網(wǎng)平差程序或軟件中存在的不足，減少了作業(yè)人員的工作量，提高了工作效率。

[1] 付江缺，高萬柱，張襲遠，等.水準(zhǔn)網(wǎng)自動處理模型及程序?qū)崿F(xiàn)[J].測繪工程，2013，22(5)：69-70.

[2] 劉忠林，薛維剛.精密跨河水準(zhǔn)測量方法與數(shù)據(jù)處理[J].測繪空間與地理信息，2009,32(4):204-206.

[3] 潘地林.Visual Basic 程序設(shè)計[M].北京：高等教育出版社，2006.

[4] 黃維彬.近代平差理論及其應(yīng)用[M].北京：解放軍出版社，1992.

[5] 李 巍，趙 亮,張占偉，等.基于VC的水準(zhǔn)網(wǎng)平差程序的設(shè)計與研究[J].測繪通報，2013,(5)：20-22.

[6] 查 明，段五杏，徐必林，等.全國GPS二級網(wǎng)數(shù)據(jù)處理和精度分析[J].全球定位系統(tǒng)，2001,26(1)：14-18.

[7] 崔希璋，於宗儔，陶本藻,等.廣義測量平差[M].2版.武漢：武漢大學(xué)出版社,2009.