趙 輝 王建偉 李思卓 樊康堯

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

0 引言

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

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

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

1 橢球面積計算方法

1.1 橢球面梯形圖幅面積

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

dx=MdBdy=NcosBdL

(1)

則

(2)

求積后得

(3)

由于公式(3)計算復雜,對式(2)中被積函數進行級數展開,再分項積分得:

(4)

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

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

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

2 數值積分橢球面積計算方法

圖2 數值積分方法任意圖斑橢球面積計算

(5)

B的積分下限為B0,積分上限P1P2為大地線,根據大地線微分方程,可建立大地線上經緯度的關系

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

(6)

對B積分

(7)

(8)

(1)矩形法:

(9)

(2)梯形法:

(10)

(3)Simpson法:

(11)

3 算例分析

以1∶10 000標準分幅大小進行實驗分析,選取I49G024001圖幅,四角點坐標見表1。參數使用2000國家大地坐標系橢球參數,長半軸a=6 378 137 m,扁率f=1/298.257 222 101。

表1 標準分幅I49G024001四角點大地坐標

如表2所示，根據公式(3)嚴密公式計算的橢球面積為26 367 263.059 3 m2,根據公式(4)計算的橢球面積同樣為26 367 263.059 3 m2,再利用公式(9)～(11)驗證數值積分方法的正確性。由公式(3)可知,橢球面積與緯度、經差有關,而與經度無關,因此,對于經線緯線圍成的標準分幅圖斑不需要劃分子區(qū)間,1個區(qū)間即可直接計算。

表2 數值積分方法計算的梯形圖斑橢球面積 單位：m2

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

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

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

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

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

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

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

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

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

4 結束語

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

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

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