趙 輝 王建偉 李思卓 樊康堯

(1.自然資源部大地測(cè)量數(shù)據(jù)處理中心, 陜西 西安 710054; 2.自然資源部第二地形測(cè)量隊(duì), 陜西 西安 710054)

0 引言

面積計(jì)算是國土調(diào)查的一項(xiàng)重要工作,對(duì)于不規(guī)則圖斑的平面面積計(jì)算已有簡單易行的公式,但高斯-克呂格投影是等角橫切橢圓柱投影,存在面積變形,離中央經(jīng)線越遠(yuǎn)變形越大,而橢球面積可避免此問題。已有許多學(xué)者對(duì)平面面積和橢球面積及不同坐標(biāo)系的變化做出深入研究[1],王解先等[2]討論了高斯投影中央子午線變化和投影面高程變化引起的面積變形,分析了橢球面積和平面面積的差異;黨亞民等[3]研究了不同坐標(biāo)系坐標(biāo)轉(zhuǎn)換對(duì)國土資源調(diào)查面積量算的影響;楊潤書等[4]計(jì)算了高原地區(qū)不同坐標(biāo)系及投影面引起的面積誤差。

地球橢球面是一個(gè)曲面,橢球面積的計(jì)算相對(duì)復(fù)雜。在計(jì)算由子午圈和平行圈圍成的橢球面梯形時(shí),對(duì)微分面積進(jìn)行積分,并將被積函數(shù)展開級(jí)數(shù)分項(xiàng)求積后再相加[5];施一民等[6]利用測(cè)地坐標(biāo)推導(dǎo)了橢球面上測(cè)地格網(wǎng)曲邊矩形的面積;林綠等[7]將橢球面上區(qū)域投影到赤道平面,利用二重積分的方式計(jì)算曲面面積。對(duì)于不規(guī)則圖斑,第三次全國國土調(diào)查技術(shù)規(guī)程中在多邊形邊上內(nèi)插一定的加密點(diǎn),將任意封閉區(qū)域分割成有限個(gè)梯形圖塊,求出各小梯形圖塊后累加得到總面積[8];史守正等[9]改進(jìn)了圖斑面積的定積分近似計(jì)算方法,簡化了橢球面積的計(jì)算過程。也可以通過求取平面面積和橢球面積的轉(zhuǎn)換關(guān)系,將平面面積轉(zhuǎn)換為橢球面積。劉洋等[10]顧及地球曲率,對(duì)高斯投影面積變形進(jìn)行拋物線擬合,然后修正高斯投影面積;茹仕高等[11]從高斯投影面積變形原理出發(fā),引入格網(wǎng)改正法建立了高斯投影面積與橢球面積的修正系數(shù)。

本文從橢球面積計(jì)算嚴(yán)密公式出發(fā),對(duì)大地坐標(biāo)形式的不規(guī)則圖斑直接采用數(shù)值積分方法計(jì)算橢球面積,并通過算例分析驗(yàn)證方法的可行性與精度。

1 橢球面積計(jì)算方法

1.1 橢球面梯形圖幅面積

由兩條子午線和兩條平行圈圍成的橢球面梯形,微分面積dP等于子午線微分弧長dx和平行圈微分弧長dy的乘積:

dx=MdBdy=NcosBdL

(1)

則

(2)

求積后得

(3)

由于公式(3)計(jì)算復(fù)雜,對(duì)式(2)中被積函數(shù)進(jìn)行級(jí)數(shù)展開,再分項(xiàng)積分得:

(4)

1.2 任意圖斑橢球面積計(jì)算的傳統(tǒng)方法

本文將文獻(xiàn)[8]的任意圖斑橢球面積計(jì)算稱為傳統(tǒng)方法,對(duì)于不規(guī)則多邊形圖斑可以分割成有限個(gè)梯形圖塊,對(duì)每個(gè)梯形圖塊進(jìn)行疊加得到總面積。如圖1所示,主要思路為:指定一條經(jīng)線L0,將多邊形P1P2邊的兩端點(diǎn)投影到經(jīng)線L0上，得到P′1P′2,圍成梯形圖塊P′1P′2P2P1,然后對(duì)P1P2邊插入一些加密點(diǎn)Pi,沿緯線再將梯形圖塊分割為多個(gè)小梯形圖塊P′iP′i+1Pi+1Pi,分隔到一定程度,小梯形面積可根據(jù)公式(4)進(jìn)行計(jì)算,累加小梯形面積得到梯形圖塊P′1P′2P2P1的面積,累加所有多邊形邊的梯形圖塊面積得到多邊形圖斑面積。

圖1 傳統(tǒng)方法任意圖斑橢球面積計(jì)算

2 數(shù)值積分橢球面積計(jì)算方法

圖2 數(shù)值積分方法任意圖斑橢球面積計(jì)算

(5)

B的積分下限為B0,積分上限P1P2為大地線,根據(jù)大地線微分方程,可建立大地線上經(jīng)緯度的關(guān)系

B=u(L,BP1,LP1,BP2,LP2)

(6)

對(duì)B積分

(7)

(8)

(1)矩形法:

(9)

(2)梯形法:

(10)

(3)Simpson法:

(11)

3 算例分析

以1∶10 000標(biāo)準(zhǔn)分幅大小進(jìn)行實(shí)驗(yàn)分析,選取I49G024001圖幅,四角點(diǎn)坐標(biāo)見表1。參數(shù)使用2000國家大地坐標(biāo)系橢球參數(shù),長半軸a=6 378 137 m,扁率f=1/298.257 222 101。

表1 標(biāo)準(zhǔn)分幅I49G024001四角點(diǎn)大地坐標(biāo)

如表2所示，根據(jù)公式(3)嚴(yán)密公式計(jì)算的橢球面積為26 367 263.059 3 m2,根據(jù)公式(4)計(jì)算的橢球面積同樣為26 367 263.059 3 m2,再利用公式(9)～(11)驗(yàn)證數(shù)值積分方法的正確性。由公式(3)可知,橢球面積與緯度、經(jīng)差有關(guān),而與經(jīng)度無關(guān),因此,對(duì)于經(jīng)線緯線圍成的標(biāo)準(zhǔn)分幅圖斑不需要?jiǎng)澐肿訁^(qū)間,1個(gè)區(qū)間即可直接計(jì)算。

表2 數(shù)值積分方法計(jì)算的梯形圖斑橢球面積 單位：m2

為了模擬任意不規(guī)則圖斑,選擇表1中ABC、ACD構(gòu)成兩個(gè)三角形。這時(shí)已經(jīng)不能用嚴(yán)密公式計(jì)算其準(zhǔn)確的橢球面積,采用1.2中傳統(tǒng)方法和本文提出的數(shù)值積分方法計(jì)算,進(jìn)行比較分析。

在區(qū)間劃分方面,采用二分的方式,第二次劃分是前一次的2倍,加快迭代過程。表3～6限于篇幅,ABC和ACD顯示面積減去了13 180 000 m2,兩個(gè)三角形的面積和減去了26 360 000 m2。

表3列出了傳統(tǒng)方法計(jì)算的三角形結(jié)果,可以看出當(dāng)劃分1個(gè)區(qū)間時(shí),計(jì)算的兩個(gè)三角形面積相等,ABC補(bǔ)的面積大于割的面積,ACD補(bǔ)的面積小于割的面積,最終面積都等于沿緯線中點(diǎn)對(duì)ABCD圖斑切割后的分割圖塊面積。隨著劃分個(gè)數(shù)的增加,ABC和ACD面積都逐漸趨于穩(wěn)定,但受變量的數(shù)值精度影響,兩個(gè)三角形的面積和發(fā)生變化。

表3 傳統(tǒng)方法計(jì)算的不規(guī)則圖斑橢球面積

表4～6列出了三種數(shù)值積分方法的計(jì)算結(jié)果顯示,相對(duì)于傳統(tǒng)方法,數(shù)值積分方法的面積和較為穩(wěn)定。矩形法和梯形法計(jì)算結(jié)果接近,劃分210個(gè)區(qū)間面積變化小于0.01 m2,而Simpson法表現(xiàn)出了優(yōu)越性,劃分24個(gè)區(qū)間面積變化小于0.01 m2,劃分28個(gè)區(qū)間面積變化小于0.001 m2,劃分210個(gè)區(qū)間面積變化小于0.000 1 m2。

表4 矩形法數(shù)值積分計(jì)算的不規(guī)則圖斑橢球面積

表5 梯形法數(shù)值積分計(jì)算的不規(guī)則圖斑橢球面積

表6 Simpson法數(shù)值積分計(jì)算的不規(guī)則圖斑橢球面積

同時(shí),結(jié)果也顯示分割的兩個(gè)三角形ACD比ABC橢球面積大10 081.51 m2,傳統(tǒng)方法按一個(gè)區(qū)間計(jì)算,則ABC面積多補(bǔ)了5 040.75 m2,ACD面積多割了5 040.75 m2。

4 結(jié)束語

針對(duì)不規(guī)則圖斑的橢球面積計(jì)算,傳統(tǒng)方法是在不規(guī)則圖斑的邊界上插入加密點(diǎn),構(gòu)造近似橢球面梯形圖塊,按規(guī)則梯形橢球面積計(jì)算公式求解后累積求和得到總面積。本文提出了基于數(shù)值積分思想的不規(guī)則圖斑橢球面積計(jì)算方法，利用實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較分析,得出以下結(jié)論:

(1)傳統(tǒng)方法為了提高計(jì)算精度需要插入一定密度的加密點(diǎn),隨著密度的增加，受數(shù)值精度影響變大,而數(shù)值積分方法影響較小;

(2)數(shù)值積分方法中Simpson法比矩形法和梯形法效率更高,能夠快速求出不規(guī)則圖斑的橢球面積,且精度可靠。