紀志剛，余 銳，宣 偉，郭劍琴

（1.天津市市政工程設計研究院，天津300051 2.武漢大學 測繪學院，湖北 武漢430079）

在許多大型線狀工程中，地球曲率引起的長度變形對項目影響顯著，需要采用高斯投影的方式對數(shù)據(jù)進行處理?，F(xiàn)有程序多基于克拉索夫斯基橢球或1975國際橢球[1,2]，不利于解算。筆者將CGCS2000橢球參數(shù)代入原始方程推出新公式，并將公式編制成程序，利用計算機進行解算，滿足了需求。

1 2000國家大地坐標系

2000國家大地坐標系的原點為包括海洋和大氣的整個地球的質量中心。坐標系的Z軸由原點指向歷元2000.0的地球參考極方向，X軸指向格林尼治參考子午線與地球赤道面（歷元2000.0）的交點，Y軸與Z軸、X軸構成右手正交坐標系[3]。采用廣義相對論意義下的尺度，2000國家大地坐標系采用的地球橢球為CGCS2000橢球[4],其基本參數(shù)如表1所示。

表1 CGCS2000橢球基本參數(shù)表[5]

2 程序實現(xiàn)

2.1 程序過程概況

該程序的實現(xiàn)過程如圖1所示。

2.2 程序概述

圖1 程序運算流程圖

該程序基于VS2008平臺用C#語言編寫，主要模塊為6個函數(shù)，分別為solve_x0()、conversion_degree()、conversion_Text()、Gauss_ZhengSuan()、Gauss_FanSuan()和Gauss_HuanDai()。

1）Double solve_x0(doubleB)，該函數(shù)的功能是利用給定的參數(shù)B，去求解赤道到緯度B所對應平行圈的子午線弧長：

式中，M為子午圈曲率半徑。將M按級數(shù)展開，并進行積分得弧長計算公式為：

將CGCS2000橢球參數(shù)代入得各系數(shù)數(shù)值，如表2所示。

表2 弧長計算公式系數(shù)表

將參數(shù)代入式（2）得弧長計算公式為：

2）conversion_degree()，該函數(shù)的功能是將rad轉換為（°）。由于坐標反算的結果為大地坐標，坐標換帶時也需要將投影坐標首先反算得到大地坐標。

3）conversion_Text()，該函數(shù)的功能是將大地坐標（只考慮經緯度，不考慮高程）轉換為rad。程序中大地坐標的輸入形式用下例說明：如大地緯度為95°23′18.3″，輸入形式為95.231 83。

4）Gauss_ZhengSuan()，該函數(shù)的功能是高斯正算，即將已知大地坐標轉換為目標帶的高斯投影坐標。Gauss_ZhengSuan()所解算出來的大地坐標，其y坐標在原來的基礎上增加了500 km?；贑GCS2000橢球元素，適合電算的高斯正算公式為[6]：

式中,

將CGCS2000坐標系參數(shù)代入式（7），然后匯編入Gauss_ZhengSuan()函數(shù)中。

5）Gauss_FanSuan()，該函數(shù)的功能是將高斯投影坐標反算為大地坐標。其所接收的大地坐標的y坐標也是在原來的基礎上增加了500 km。Gauss_FanSuan()函數(shù)的核心為高斯反算公式，基于CGCS2000橢球元素的高斯反算公式為：

接下來求解Bf。用迭代法，首先令

式中，X為已知投影點X坐標，6 367 449.145 71為式（5）的第1個系數(shù)。每次按式（11）迭代計算，直至

6）Gauss_HuanDai()，該函數(shù)的功能是高斯換帶，即將已知帶投影坐標換算為目標帶坐標。換帶類型不同，所得結果會有所不同。常規(guī)的分帶有6°帶和3°帶2種，函數(shù)根據(jù)分帶的不同，又細分為4種類型：6°帶換算為 6°帶、6°帶換算為 3°帶、3°帶換算為 3°帶、3°帶換算為6°帶[7]。其中6°帶和3°帶之間的相互轉換具有相對復雜性，下面重點介紹二者的相互轉換過程。

① 6°帶換算為3°帶。以圖2為例，首先在6°帶內利用已知點A的坐標，采用Gauss_FanSuan()函數(shù)計算出A點的緯度B以及經度與中央子午線的經差l。若l<1.5°，在3°帶內，A點位于以117°為中央子午線的投影帶內，此時l不變，再利用Gauss_ZhengSuan()函數(shù)計算出A點在該帶內的投影坐標；若l>1.5°，在3°帶內，A點位于以120°為中央子午線的投影帶內，此時l=6°-l，再利用Gauss_ZhengSuan()函數(shù)計算出A點在該帶內的投影坐標。

② 3°帶換算為6°帶。首先利用Gauss_FanSuan函數(shù)計算出A點的緯度B和經度與中央子午線的經差l。判斷A點所在3°帶中央子午線是否滿足L0=6n-3。若滿足，則A點坐標不變；若不滿足，l=1.5°+l，再利用Gauss_ZhengSuan函數(shù)計算出A點在該帶內的投影坐標。程序實例圖如圖3、圖4所示。

圖3 程序實例圖（1）

圖4 程序實例圖（2）

3 算例分析

為檢驗轉換結果的準確性，下面均采用Powercoor軟件轉換結果作參考值，同類型的轉換結果進行比較，以獲取轉換精度。為了同軟件中的顯示格式相同，各表中經緯度均采用如下表示形式：30°30′15.2″表示為30.301 52。以下表中初始數(shù)據(jù)均為模擬數(shù)據(jù)。

表3 正算結果分析表/m

3.1 正反算數(shù)據(jù)表及結果分析

坐標分量中誤差計算公式分別為：

點位中誤差為：

對表3各點分別作中誤差分析，得出最大點位中誤差為0.011 49 m。

表3和表4中，原始大地坐標經過正反算后又得到新的大地坐標，將二者進行比較，結果如表5所示。

表4 反算結果分析表

表5 新舊大地坐標對照表

表3～表5中數(shù)據(jù)顯示，在X方向和B方向轉換結果較參考值均偏大，Y方向和L方向轉換結果較參考值均偏小，存在一定的系統(tǒng)誤差。

3.2 換帶計算分析表及結果分析

將坐標（5 635 146.875 0，629 151.817 9）分別以6°→ 6°，6°→ 3°，3°→ 6°，3°→ 3°4 種方式進行換帶計算，并將結果與參考值進行比較，見表6。對表6結果同樣利用式（11）進行點位中誤差分析，經計算得各點轉換后最大中誤差為1.047 cm。

表6 換帶計算表

4 結 語

綜上所述，坐標正算結果中誤差最大差值可達到1.149 cm,一般情況其精度在mm級；反算得到的大地坐標精度可以達到10-4（″）量級；坐標換帶精度一般也可以達到mm級，基本滿足了工程建設的需要。

[1]林曉靜,張小紅.ITRF2005與CGCS2000坐標轉換方法與精度分析[J].大地測量與地球動力學,2010,30(2):117-119

[2]唐運海,何誠.坐標轉換及換帶計算的研究與實驗分析[J].測繪與空間信息,2011,34（2）:1-5

[3]李征航,黃勁松.GPS測量與數(shù)據(jù)處理［M］.武漢:武漢大學出版社,2005

[4]程鵬飛,文漢江.2000國家大地坐標系橢球參數(shù)與GRS80和WGS84的比較[J].測繪學報,2009,38(3):89-94

[5]孔祥元,郭際明,劉宗泉.大地測量學基礎[M].武漢:武漢大學出版社,2007

[6]劉飛,周琳琳.GPS大地坐標向地方坐標轉換的實用方法研究[J].華東師范大學學報,2005(1):73-77

[7]胡圣武.對高斯投影分帶的研究[J].地理空間信息, 2012,10(1):54-56