基于最小二乘和窮舉算法的高程抵償面優(yōu)化

鄧標(biāo), 李鵬, 李德亮, 王延霞

(滁州學(xué)院 地理信息與旅游學(xué)院,安徽 滁州 239000)

0 引 言

大地測量計(jì)算的基準(zhǔn)面是參考橢球面,而橢球面是一個(gè)數(shù)學(xué)曲面,在其面上進(jìn)行計(jì)算是相當(dāng)復(fù)雜的,若能把橢球面上的元素歸算到平面上來,并在平面上進(jìn)行計(jì)算,問題就簡單得多[1]．當(dāng)前解決由橢球面到平面的投影方法有多種,但無論采用何種方式都會(huì)產(chǎn)生非線性變形問題．因高斯投影具有保角映射的特點(diǎn),目前國際上普遍采用此方法進(jìn)行計(jì)算[2],當(dāng)然高斯投影方法也不可避免地出現(xiàn)投影變形問題．眾所周知,在距中央子午線較遠(yuǎn)或海拔較高的區(qū)域,國家控制網(wǎng)的邊長與地面實(shí)際邊長之間產(chǎn)生的差值較大[3],不能滿足文獻(xiàn)[4]對(duì)控制網(wǎng)長度變形小于2.5 cm/km的要求．因此,如何消除或減弱邊長變形問題,是科研工作者和工程技術(shù)人員必須面臨的現(xiàn)實(shí)問題．

1 高斯投影變形分析與處理方法

1.1 高斯投影變形分析

高斯投影的實(shí)質(zhì)是將橢球面按一定的數(shù)學(xué)法則投影為平面．在高斯投影計(jì)算過程中,決定高斯投影的計(jì)算參數(shù)有多項(xiàng),但起決定作用的參數(shù)只有中央子午線和高程面兩項(xiàng)[5]．高斯投影時(shí)需將地面觀測長度歸化到國家統(tǒng)一的參考橢球面上,由此產(chǎn)生的變形稱為高程歸化改正,其計(jì)算公式為[12]

(1)

式中:D為實(shí)測邊長度，m;ΔS1為測距邊長歸化到參考橢球面上的改正數(shù)，m;RA為測距邊參考橢球面曲率半徑，m;Hm為測距邊高出大地水準(zhǔn)面的平均高程，m．ΔS1始終為負(fù)值,說明將地面觀測值歸化到參考橢球面時(shí),長度總是縮小的,并且ΔS1和Hm成正比．

將參考橢球面上的邊長歸算到高斯平面上發(fā)生的變形叫高斯投影變形,其產(chǎn)生的變形值計(jì)算公式為[13]:

(2)

式中:S為參考橢球面上邊長，m;ym為測距邊兩端點(diǎn)近似橫坐標(biāo)的平均值，m;Rm為參考橢球面在測距邊中點(diǎn)的平均曲率半徑，m．ΔS2為正值,表明將參考橢球面的長度歸化到高斯平面上時(shí),長度總是增長的,并且ΔS2與ym的平方成正比．由于D值和S值對(duì)兩個(gè)改正數(shù)影響較小,并用Rm替代RA則兩項(xiàng)歸算綜合變形可表示為[14]

(3)

由式(3)可知:兩次投影變換之間存在一定程度的抵消,如何使抵消程度最大化,是研究投影變形的關(guān)鍵問題所在[15]．

1.2 處理方法

當(dāng)前,在工程實(shí)踐中,應(yīng)用比較廣泛的處理方法主要有以下幾種,現(xiàn)歸納總結(jié)如下:

1)抵償坐標(biāo)系:它是根據(jù)高斯投影改正選擇一個(gè)合適的高程參考面作為整個(gè)或部分測區(qū)的投影面,使得由高程歸化產(chǎn)生的改正抵償高斯投影變形改正．簡言之,確定抵償坐標(biāo)系的思想就是用適當(dāng)?shù)母叱袒鶞?zhǔn)面以抵償分帶投影變形[16]．工程中常選擇測區(qū)范圍內(nèi)最大橫坐標(biāo)處為抵償點(diǎn),即y0=ymax,使得該處的高程歸化改正ΔS1與高斯投影變形ΔS2相等,那么該處就得到了完全的補(bǔ)償,得到的抵償投影面為

(4)

3)抵償任意帶坐標(biāo)系:它是通過既改變Hm,又改變ym,來共同抵償兩項(xiàng)歸算改正變形[17]．

2 基于最小二乘法確定最佳抵償投影面

由前文分析可知,現(xiàn)有各種算法都有一定的局限性,因此,本文采用最小二乘法結(jié)合窮舉算法求取測區(qū)最優(yōu)抵償高程面和最佳中央子午線．盡管目前流行的“一點(diǎn)一方位”法能夠解決部分工程的實(shí)際問題,但是在相對(duì)高差較大時(shí),該方法也易于造成變形超限的問題．因此,討論最佳抵償投影面仍具有現(xiàn)實(shí)的指導(dǎo)意義和實(shí)用價(jià)值．

最小二乘法推導(dǎo)最佳中央子午線和最優(yōu)抵償高程面方法如下:假定最優(yōu)抵償高程面為H0,最佳中央子午線在國家標(biāo)準(zhǔn)3°高斯投影面中值為Y0．那么由式(3)可知

(5)

(6)

進(jìn)而對(duì)于抵償高程面和移動(dòng)后中央子午線的要求為

(7)

H0)].

(8)

(9)

式(8)、式(9)兩項(xiàng)求和并根據(jù)極值定理可得:

[(Yi-Y0)2-2Rm(Hi-H0)]=0.

(10)

式(10)經(jīng)整理后得到:

(11)

即

Y0)-Rm]=0,

(12)

因?yàn)閇(Yi-Y0)-Rm]≠0,所以由式(12)可得:

(13)

等式兩邊同除以n則式(13)可化為

(14)

Hi∈[Hmin,Hmax]時(shí),有:

(15)

同理當(dāng)Yi∈[Ymin,Ymax],可得

(16)

(17)

式中:Ymax為測區(qū)控制點(diǎn)東坐標(biāo)最大值；Ymin為測區(qū)控制點(diǎn)東坐標(biāo)最小值;Hmax為測區(qū)控制網(wǎng)高程最大值;Hmin為測區(qū)控制網(wǎng)高程最小值．根據(jù)式(13)～(17)可得:

(18)

(19)

(20)

同理可得:

(21)

由式(21)可得出,H0的最大值為當(dāng)

(22)

式(22)為H0的極限,這和式(18)中根號(hào)里面值大于零的限制是有同樣要求的．而Y0離項(xiàng)目區(qū)中點(diǎn)越遠(yuǎn)時(shí),H0的值則越來越小．由式(21)可得,若Y0=0,既中央子午線為國家標(biāo)準(zhǔn)3°帶時(shí),可得最優(yōu)抵償高程面為

(23)

若在地勢平坦且海拔較低的的平原區(qū)域,取H0=0,可得最佳中央子午線值:

(24)

由式(24)可知,最優(yōu)抵償高程面是和最佳中央子午線相互關(guān)聯(lián)的,只有取得最佳中央子午線后才能獲得最優(yōu)抵償高程面,同理在確定了抵償高程面H0時(shí),可得中央子午線值．但在式(24)中存在兩個(gè)變量H0和Y0,且不存在具有實(shí)際意義的自由極值,所以難以求解．對(duì)于求算極值的問題,文獻(xiàn)[18-19]中首次分別研究多變量極值的搜索法,但是相對(duì)模型簡單,適用于一般的定參數(shù)對(duì)象的線性化問題,且受算法穩(wěn)定性的要求,收斂速度較慢．為此,本文引用人工智能中的窮舉算法進(jìn)行求解最優(yōu)值．

3 窮舉算法確定最優(yōu)值及算例

3.1 窮舉算法原理

根據(jù)以上原理,以GPS控制點(diǎn)所處位置的投影變形代表整個(gè)控制區(qū)的投影變形,設(shè)計(jì)計(jì)算參數(shù)如下:

1)初始投影參數(shù):抵償高程面為H0=Hmin,中央子午線由式(18)確定;

2)計(jì)算范圍:抵償高程面范圍為

3)計(jì)算步長:抵償高程面步長為1 m,中央經(jīng)線步長為0.1 km,按照式(19)和式(20)進(jìn)行計(jì)算后,取值選擇至0.1 km;

需要說明的是,H0和Y0值的搜索范圍可根據(jù)控制區(qū)域的情況人為選擇,依據(jù)測區(qū)的實(shí)際情況加以確定,但高程H0不要超越理論高程的極值

在使用過程中,根據(jù)測區(qū)情況,由式(22)可以較快地設(shè)定高程的取值范圍,在MATLAB中可以較快地搜索出最佳中央子午線和最優(yōu)抵償高程面的取值．

3.2 算例分析

為說明以上原理的實(shí)用性及有效性,下面通過相應(yīng)實(shí)例來驗(yàn)證．實(shí)例數(shù)據(jù)取自文獻(xiàn)[20]:藏區(qū)某流域位于東經(jīng)98°04′～98°36′,北緯31°11′～31°51′之間,中央子午線經(jīng)度為99°,測區(qū)平均海拔超過3 300 m,且測區(qū)高差相差較大,各測區(qū)實(shí)測數(shù)據(jù)如表1所示．

表1 某流域各測區(qū)實(shí)測數(shù)據(jù)

在進(jìn)行控制網(wǎng)平差前,必須先估算高程歸化改正和高斯投影變形,保證兩項(xiàng)綜合長度變形值不大于2.5 cm/km,根據(jù)式(1)、(2)計(jì)算各測區(qū)綜合變形,結(jié)果如表2所示．

表2 國家坐標(biāo)系下各測區(qū)綜合變形cm/km

測區(qū)高程歸化改正高斯投影改正綜合投影改正A-59.110.12-58.99B-58.530.18-58.35C-57.370.22-57.15D-58.030.30-57.73E-56.620.36-56.26F-55.120.59-54.54G-58.3580.67-57.69

由表2可知,藏區(qū)某流域長度變形以高程歸化改正為主,采用國家標(biāo)準(zhǔn)3°帶高斯投影,各測區(qū)綜合長度變形值遠(yuǎn)大于文獻(xiàn)[4]的要求,故不能采用國家標(biāo)準(zhǔn)3°帶高斯投影．

文章首先采用抵償高程面法進(jìn)行求取,而不改變中央子午線．由式(21)可得H0=3616 m,根據(jù)式(3)計(jì)算各測區(qū)的綜合變形．為便于對(duì)比分析,本文與常規(guī)方法進(jìn)行對(duì)比,取測區(qū)平均高程面(H0=3669 m),國家標(biāo)注中央子午線Y0=0,計(jì)算結(jié)果如表3所示．

表3 抵償高程面計(jì)算結(jié)果比較cm/km

測區(qū)傳統(tǒng)常規(guī)方法高程歸化改正高斯投影改正綜合投影改正本文方法高程歸化改正高斯投影改正綜合投影改正A-1.5 0.1-1.4-2.4 0.2-2.2B-0.9 0.2-0.8 -1.84 0.2-1.6C 0.2 0.2 0.4-0.6 0.2-0.4D-0.4 0.3-0.1-1.3 0.3-1.0E 1.0 0.4 1.3 0.1 0.4 0.5F 2.5 0.6 3.1/超限 1.6 0.6 2.2G-0.7 0.66-0.1-1.6 0.66-0.9

由表3可得:傳統(tǒng)方法和本文方法相比,本文方法計(jì)算出來的綜合投影改正較傳統(tǒng)方法誤差分布更均勻．傳統(tǒng)方法對(duì)于高差較大的地區(qū),易造成變形超限情況的發(fā)生．

在驗(yàn)證本文推導(dǎo)的抵償高程面法后,選取抵償高程面及任意帶高斯投影建立該流域獨(dú)立坐標(biāo)系,利用式(18)和H0的限制條件,采用窮舉算法得測區(qū)最優(yōu)值H0=3589 m,Y0=41.7 km．即抵償投影面高程為3 589 m,將該流域99°中央子午線向東平移測區(qū)值41.7 km,計(jì)算獨(dú)立坐標(biāo)系下各測區(qū)的綜合長度變形,計(jì)算結(jié)果如表4所示．

表4 最小二乘法各測區(qū)綜合變形cm/km

測區(qū)高程歸化改正高斯投影改正綜合投影改正A-2.781.25-1.53B-2.201.09-1.10C-1.040.99-0.04D-1.700.84-0.85E-0.280.75 0.46F 1.210.49 1.69G-2.020.42-1.60

由表4可知,在獨(dú)立坐標(biāo)系下,投影到抵償高程面的高程歸化改正,很好地抵償了長度變形問題,各測區(qū)綜合變形均小于2.5 cm/km,滿足整個(gè)流域投影變形需求．

表5 測區(qū)綜合變形比較 cm/km

4 結(jié) 論

高海拔區(qū)域長度投影變形以高程歸化改正為主,采用國家標(biāo)準(zhǔn)3°帶高斯投影,各測區(qū)綜合長度變形值遠(yuǎn)大于現(xiàn)行工程測量規(guī)范的限差要求．利用最小二乘法在獲取最優(yōu)抵償高程面和最佳中央子午線時(shí),考慮了高程和中央子午線距離等參數(shù),非常符合地理位置不同、地形不同、取值不同的實(shí)際情況,對(duì)限制投影變形有較為明顯的優(yōu)勢．

最小二乘法和窮舉算法結(jié)合能快而準(zhǔn)確地找到最優(yōu)抵償高程面和最佳中央子午線,在滿足相關(guān)規(guī)范的前提下,本文方法計(jì)算的誤差分布較為均勻．采用最小二乘法推導(dǎo)的H0和Y0間的關(guān)系及H0取值范圍的約束,很大程度上降低了窮舉算法的計(jì)算量,避免了全局搜索的缺點(diǎn)．

最小二乘法結(jié)合窮舉算法求取最優(yōu)投影面和最佳中央子午線法,可作為工程測量中控制投影變形的一種參考和借鑒．在實(shí)際工程應(yīng)用中,使用本文方法時(shí),Y0值的搜索范圍可以根據(jù)工程實(shí)際選擇,但應(yīng)注意曲面代替水平面的適用范圍．