自定義坐標系的建立及其坐標變換實現(xiàn)

周朝憲，董少波，和志軍，姜高珍，燕丹晨

(1.中國地質(zhì)大學(xué)，北京 100083;2.南凱礦業(yè)有限公司，北京 100023;

3.中色地科礦產(chǎn)勘查股份有限公司，北京 100012;4.中國地質(zhì)科學(xué)院，北京 100037;5.國家海洋環(huán)境預(yù)報中心，北京 100081)

1 引言

20世紀90年代以來，地質(zhì)勘查工作的一大進展就是地理信息系統(tǒng)(GIS)越來越多地滲入并支撐著當今的勘查工作?，F(xiàn)今的地質(zhì)勘查工作要求地質(zhì)勘查人員掌握投影、大地基準(datum)和數(shù)學(xué)橢球體(ellipsoid)的含義和關(guān)系，以免造成錯誤和不必要的問題。

所謂投影，包括兩部分內(nèi)容:一是對地球這個橢球體參數(shù)最佳近似的理論化，即構(gòu)建大地基準;二是把這一理論化的地球大地基準體上的坐標點(大地坐標系的經(jīng)緯度)轉(zhuǎn)到平面坐標(即方里網(wǎng))上(周朝憲等，2013)。大地基準實際也包括兩部分內(nèi)容:首先把實際地球先簡單化為數(shù)學(xué)橢球體。這個數(shù)學(xué)橢球體由地球半徑長軸和扁平率(flattening，扁平率為橢圓率(ellipticity)或者偏心率(eccentricity)的函數(shù))來定義。由于實際地球和數(shù)學(xué)橢球體還存在比較大的差異(主要是地球質(zhì)心和數(shù)學(xué)橢球體并不嚴格重合)，難以滿足測繪和各種工程所需，所以大地基準還包括對數(shù)學(xué)橢球體進行地球質(zhì)心的偏移補償。ellipsoid、datum和實際地球表面的關(guān)系見圖1。

地球上任一點的經(jīng)緯度都是建立在datum基礎(chǔ)上，而不是建立在ellipsoid之上。任一點的經(jīng)緯度在不同datum體系下的數(shù)值基本不相同，通常相差幾十米乃至數(shù)百米。了解這一點十分重要，有些地質(zhì)人員在地質(zhì)圖中沒有給出datum，這會給使用圖件的人帶來很大麻煩。

目前，GPS全球定位系統(tǒng)所采用的UTM(Universal Transverse Mercator，通用橫軸墨卡托)投影其datum是WGS(World Geodetic System)84橢球體。理論上講，投影方法和datum可以自由組合，從而產(chǎn)生眾多的坐標系統(tǒng)。另外，由于一些歷史原因，在不同的歷史時期，許多國家和地區(qū)都建立有各自的基準datum和ellipsoid。如20世紀50年代以來，我國使用過與北京54坐標系對應(yīng)的Krassovsky橢球體和西安80坐標系的IAG1975橢球體。

圖1 地球?qū)嶋H地形、大地基準和數(shù)學(xué)橢球體的關(guān)系Fig.1 The relationship among the Earth actual surface，Datum/Geoid and Ellipsoid

實際地質(zhì)勘查中，如果地質(zhì)人員不了解這些坐標知識，會給工作帶來一些麻煩。如某CXB地質(zhì)隊在埃塞俄比亞工作時，發(fā)現(xiàn)手持式GPS給出的坐標和當?shù)氐牡匦?地質(zhì)圖不相吻合，便按照國內(nèi)處理北京54坐標系和西安80坐標系的辦法，根據(jù)當?shù)氐匦螆D上的控制點，建立了使用WGS 84為橢球體(手持式GPS接收機坐標所依據(jù)的橢球體)，并自行修改手持GPS的D x、D y和D z。這實際是建立了獨特的TM投影，而地形圖上標注著其使用的datum和橢球體分別是Adindan和Clarke_1880。在已知圖件的投影和datum的情況下，正確的做法有2種:(1)把地形圖通過坐標變換到WGS 84坐標下，以便和手持式GPS接收機數(shù)據(jù)套合。因為有些手持式GPS接收機往往沒有Adindan和Clarke_1880體系;(2)把手持式GPS坐標下的地質(zhì)成果通過GIS軟件轉(zhuǎn)換成Adindan和Clarke_1880坐標系，從而把不同坐標體系下的數(shù)據(jù)進行套合。

為了使用該地質(zhì)隊自定義坐標體系下的地質(zhì)成果，可在GIS軟件中根據(jù)其datum采用WGS 84，以及D x、D y和D z建立了該自定義坐標系，從而利用軟件快捷地進行不同坐標體系的轉(zhuǎn)換，解決不同坐標體系下地質(zhì)數(shù)據(jù)的套合使用。

周朝憲等(2013)給出了使用MapInfo?對UTM投影和Gauss-Krüger投影下方里網(wǎng)坐標和其經(jīng)緯度的轉(zhuǎn)換過程。其實，對于國內(nèi)北京54坐標系和西安80坐標系(都采用Gauss-Krüger投影，UTM投影和Gauss-Krüger投影都屬于TM橫軸墨卡托投影)，由于該兩種坐標系只有數(shù)學(xué)橢球體而沒有datum，必須先在GIS軟件中建立自定義坐標系，即按照自定義坐標系的方式處理即可。否則，基本無法利用諸多通用GIS軟件對北京54坐標系以及西安80坐標系下的地質(zhì)數(shù)據(jù)進行處理。

即使在國內(nèi)，地質(zhì)勘查人員也經(jīng)常面臨坐標變換的問題。如北京54坐標系和西安80坐標系之間的相互轉(zhuǎn)換以及和其他投影坐標系(如CGCS2000以及UTM等等)間的變換。因而，建立自定義坐標系以及實現(xiàn)不同坐標系下數(shù)據(jù)的套合使用已經(jīng)迫切擺在地質(zhì)勘查人員面前。

熊忠招(2010)、陳士銀(1997)、暢開獅(2008)等從測繪的角度探討了如何建立個別的自定義坐標系，但其方法不適合于地勘工作;楊啟和(1986)、夏蘭芳等(2007)、沈本忠(1986)、Osborne(2008)、Dozier(1980)和Kawase(2011)等就坐標變換給出了嚴格的數(shù)學(xué)換算，但是這些運算對絕大多數(shù)地質(zhì)勘查人員而言，不具有可操作性;還有人提出了簡易的算法公式，但由于簡易公式所得結(jié)果的精度太差，在實際工作中是不可以接受的(見周朝憲等，2013);劉健等(2005)、王寶軍等(2010)發(fā)展出一些專門的軟件程序用以坐標變換，但是這些軟件一則普及性很低，廣大地勘人員一般無法使用，二則也面臨構(gòu)建datum和投影等問題，并不適用于廣大地勘人員使用。目前，對廣大地質(zhì)人員而言，比較可行的辦法就是使用現(xiàn)成的比較通用的GIS軟件來進行轉(zhuǎn)換(周朝憲等，2013)，而坐標轉(zhuǎn)換的前提往往就是在這些軟件中建立起自定義坐標系。

常用的GIS軟件包括 MapInfo?、Geosoft?Oasis MontajTM、ArcGIS?以及國內(nèi)較為通用的 MapGIS?。MapGIS?6.7版本(現(xiàn)在最常用的版本)及其以前版本的軟件由于其顯示的坐標僅僅是圖面坐標，還算不上是完整的勘查GIS軟件，盡管其描圖功能較為強大。完整的勘查GIS軟件至少應(yīng)具有三個特點:(1)圖面顯示坐標為該點在所采用的GIS坐標體系中的真實坐標數(shù)值。當然，其不隨著圖面比例尺變化而變化。(2)含有目前常見的坐標體系，并可以進行坐標相互轉(zhuǎn)換。(3)基本可以在其中建立自定義坐標系。MapInfo?、Geosoft?Oasis MontajTM和ArcGIS?都具備這三個特點，這些軟件在非洲等國家地勘人員中的普及程度也遠遠高于我國，建議我國地質(zhì)勘查人員盡快熟悉。

本文使用西方地質(zhì)勘查應(yīng)用最廣的MapInfo?和Geosoft?Oasis MontajTM進行展示，如何建立自定義坐標，即如何在這些軟件中建立自定義datum，并進行坐標變換。實際上，ArcGIS?也可以定義自定義坐標系，也可以使用其處理我國的北京54坐標系和西安80坐標系下的GIS數(shù)據(jù)。步驟與MapInfo?和Geosoft?Oasis MontajTM類似，其投影信息存儲在數(shù)據(jù)集的PRJ文本文件中，也是通過直接修改該文本文件來建立自定義坐標系。陳悟天(2010)和周朝憲等(2013)認為，ArcGIS?沒有國內(nèi)現(xiàn)用的 Gauss-Krüger投影坐標的說法是不準確的，可以在該軟件中建立國內(nèi)各個地區(qū)的datum，然后進行坐標變換。

2 使用MapInfo建立自定義坐標系及其坐標變換

MapInfo?的坐標文件主要是mapinfow.prj，mapinfow.prj是個文本文件，可以直接編輯，其結(jié)構(gòu)參見表1。

MapInfo?自定義投影參數(shù)中，對于橫軸墨卡托投影坐標(UTM投影和Gauss-Krüger投影都屬于橫軸墨卡托投影)而言，主要有如下參數(shù):

投影代號(Type)，基準面(datum)，單位(U-nit)，原點經(jīng)度(Origin Longitude，對 UTM 和 Gauss-Krüger投影而言，其原點經(jīng)度皆為該帶的中央經(jīng)線的經(jīng)度值)，原點緯度(Origin Latitude，為赤道緯度值0)，比例因子(Scale Factor，對UTM投影和Gauss-Krüger投影而言，比例因子分別為0.9996和1)，東偽偏移(False Easting)，北偽偏移(False Northing)。關(guān)于UTM和Gauss-Krüger投影的上述參數(shù)更多內(nèi)容具體可參見周朝憲等(2013)。

在MapInfo?中，橫軸墨卡托投影的代號為8，其余各個投影的代號參見MapInfo?的說明書。各個基準面datum的MapInfo?代號也參見MapInfo?說明書。如，WGS 84 datum的代號是104。

例如:對東經(jīng)34°和北緯8°的點，其 UTM投影選取WGS 84 world作為大地基準(其數(shù)學(xué)橢球體是WGS 84)的MapInfo?參數(shù)為:

"UTM Zone 36，Northern Hemisphere(WGS 84)p32636"，8，104，7，33，0，0.9996，500000，0

即:UTM的36帶，投影方法為橫軸墨卡托投影，代號為8;datum為WGS 84，代號為104;坐標單位為m，代號為7;該帶屬于UTM的36帶，其起始經(jīng)線(即中央經(jīng)線)為33°E;其起始緯度為0°;UTM投影的比例因子是0.9996;東偽偏移為500，000m;北偽偏移為0m。

如上文所述，北京54坐標系和西安80坐標系沒有datum，這就要求我們在datum這一項中予以定義。在MapInfo?中，基準的定義主要有數(shù)學(xué)橢球體與相應(yīng)的參數(shù)組成(即三參數(shù)或七參數(shù))。北京54坐標系采用的是Krassovsky橢球體，其在MapInfo?中的代號是3，WGS 84數(shù)學(xué)橢球體的MapInfo?代號是28。

例如，河南新縣北部某地的大地坐標(WGS 84)為 31°46’06.9’’N 和 114°36’41.6’’E，其北京 54坐標系該3度帶(帶號38，中央經(jīng)線114°E，三參數(shù))的平移校正D x= －9.0，D y= －113和D z= －39，那么其datum在MapInfo?中格式定義為:(999，3，－9.0，－113，－39)。其中 999 代表自定義 datum，3 代表Krassovsky橢球體。

則可以寫出下列語句:

"－－－ Xinxian(D x:－9.0，D y:－113，D z:－39;Beijing_54_Krasovsky)－－－"，8，999，3，－9，－113，－39，7，114，0，1，500000，0

……………… 語句(1)

其中:

8代表投影類型為橫軸墨卡托投影，這里即為Gauss－Krüger投影;

999表示自定義datum，其后的3代表 Krassovsky橢球體，－9.0，－113，－39分別為 D x，D y和D z;

7代表單位為m;

114代表起始經(jīng)線(對橫軸墨卡托投影即為該帶的中央經(jīng)線)為114°E;

0代表起始緯線為赤道;

1代表比例因子為1;

500000代表東偽偏移500，000m;

0代表北偽偏移0m。

表1 MapInfo?中坐標系統(tǒng)各個參數(shù)定義格式Table 1 The parameter form of coordinate system in Map Info?

然后，打開MapInfo?軟件所在的目錄下的mapinfow.prj文件，把上述語句(1)寫入，保存，使用MapInfo?軟件調(diào)用該投影即可。關(guān)于不同坐標系的轉(zhuǎn)換過程步驟可參見周朝憲等(2013)。

對上文提到的某CXB地質(zhì)隊在埃塞俄比亞西部Abobo的地區(qū)，其大地坐標(WGS 84)為7°39’N和34°39’E，采用的橢球體是WGS 84，屬于6度帶，該區(qū)采用三參數(shù)格式:D x=－100，D y=2和D z=－1，該6度帶中央經(jīng)線為33°E。按照上文所述在mapinfow.prj中加入如下語句即可建立該坐標系。

"－－－ Abobo_ChuanXiBei(D x:－100，D y:2，D z:－1;Local_WGS 84)－－－"

"CXB(WGS 84)_CM33_Local_WGS 84"，8，999，28，－100，2，－1，7，33，0，0.9996，500000，0 …………語句 (2)

其中:

8代表橫軸墨卡托投影;

999表示自定義datum，其后的28代表MapInfo?中 WGS 84 數(shù)學(xué)橢球體，－100，2，－1 分別為 D x，D y和D z;

7代表方里網(wǎng)單位為m;

33代表起始經(jīng)線(對橫軸墨卡托投影即為中央經(jīng)線)為 33°E;

0代表起始緯線為赤道;

0.999 6代表比例因子為0.999 6;

500 000代表東偽偏移500，000m;

0代表北偽偏移0m。

以上是以三參數(shù)為例，對七參數(shù)的自定義坐標系，其 datum在 MapInfo?中格式定義為:(9999，3，D x，D y，D z，R x，R y，R z，K，0)。其中 9999 代表自定義datum，3代表Krassovsky橢球體，R x、R y和 R z為旋轉(zhuǎn)校正，K為比例校正。這里需要說明的是:在mapinfow.prj文件中添加自定義坐標語句時，語句應(yīng)盡可能短，否則MapInfo?調(diào)用該語句時經(jīng)常會出現(xiàn)錯行，從而引起混亂。

因而，使用MapInfo?軟件設(shè)立自定義坐標系的步驟歸納為:

(1)確立投影類型及其投影的相關(guān)參數(shù)(比例因子、起始經(jīng)緯度和經(jīng)向和緯向偽偏移量)。投影類型的MapInfo?代號參見MapInfo?Professional User Guide①。

(2)確定datum。確立datum就要確立數(shù)學(xué)橢球體和datum的相關(guān)參數(shù)。橢球體的MapInfo?代號參見Geosoft?Oasis MontajTM軟件下的ellipsoid.csv文件。

建立自定義坐標系后，便可以使用MapInfo?實現(xiàn)坐標變換，也可以使用建立自定義坐標系后的Geosoft?Oasis MontajTM來實現(xiàn)坐標變換，后者用戶界面更為友好一些。

3 使用Geosoft?Oasis Montaj TM建立自定義坐標系

地質(zhì)勘查工作一般不會建立全新的投影方法，通常是建立適合于各個地區(qū)的datum。建立自定義坐標系的核心就是建立適合于工作區(qū)的datum。下面簡單敘述如何使用Geosoft?Oasis MontajTM軟件建立埃塞俄比亞西部Abobo地區(qū)的自定義datum。

首先，修改 Geosoft?Oasis MontajTM目錄下的“datumtrf.csv”，在其中加入如表2內(nèi)容:

然后，修改 Geosoft?Oasis MontajTM目錄下的“l(fā)datumtrf.csv”，在其中加入表3內(nèi)容。

因為無論MapInfo?還是Geosoft?Oasis MontajTM對datum的命名都遵循EPSG(European Petroleum Survey Group)/POSC(Petrotechnical Open Software Corporation)命名規(guī)則。自定義的datum不屬于其范圍，故此前面加上*以示區(qū)別，如埃塞俄比亞西部Abobo地區(qū)的自定義datum為“*Abobo special”，在表2和表3中該名稱必須完全一致。

表2 Geosoft?Oasis Montaj TM中datumtrf.csv各個參數(shù)定義格式②Table 2 The parameter ofcustom coordinate system in datumtrf.csv of Geosoft? Oasis Montaj TM②

表3 Geosoft?Oasis Montaj TM中自定義坐標ldatum.csv中各個參數(shù)定義格式②Table 3 The parameter of custom coordinate system in ldatum.csv of Geosoft? Oasis Montaj TM②

圖 2 Geosoft? Oasis Montaj TM界面Fig.2 Coordinate dialogue box of Geosoft?Oasis Montaj TM

建立自定義坐標系后，可以把坐標數(shù)據(jù)導(dǎo)入Geosoft?Oasis MontajTM。第一步:通過圖2菜單上的“坐標”選擇子菜單“坐標系”，選擇我們上文新建的“*Abobo special”，把這些數(shù)據(jù)設(shè)置成以X_CXB和Y_CXB分別為經(jīng)向和緯向坐標的“*Abobo special”坐標系。第二步:通過圖2的菜單上的“坐標”選擇“新建投影坐標系”，以X_CXB和Y_CXB分別為當前的經(jīng)向和緯向坐標(大地基準為“*Abobo special”)，得出 WGS 84 經(jīng)緯度(即 Long_WGS 84 和Lat_WGS 84)。圖2中的經(jīng)緯度單位是度分秒，如圖 2中 Point 19行對應(yīng)的經(jīng)度值34.46.15.8375898670 為 34°46’15.8375898670’’。第三步:同第二步類似，以 Long_WGS 84和 Lat_WGS 84為當前坐標數(shù)值(大地基準為WGS 84)，選擇“新建投影坐標系”，新建UTM WGS 84坐標系，即得出圖2中UTM WGS 84經(jīng)向(X_UTM)和緯向(Y_UTM)坐標。

4 結(jié)論

投影坐標包括建立大地基準(datum)和把地球橢球體上的坐標點轉(zhuǎn)到平面坐標(即方里網(wǎng))上。大地基準(datum)包括數(shù)學(xué)橢球體(ellipsoid)和根據(jù)地球質(zhì)量對該數(shù)學(xué)橢球體的補償。數(shù)學(xué)橢球體、大地基準和投影方式三者可以組合成多種投影坐標系。地質(zhì)勘查人員經(jīng)常會遇到坐標轉(zhuǎn)換問題和建立自定義坐標問題。我國北京54坐標系和西安80坐標系地質(zhì)數(shù)據(jù)的GIS化，基本上需要建立自定義坐標。因為其只有數(shù)學(xué)橢球體，而缺乏大地基準。建立自定義坐標靠手工不可能實現(xiàn)，必須借助軟年實現(xiàn)。GIS軟件已經(jīng)成為地質(zhì)勘查的基本工具，而地質(zhì)數(shù)據(jù)的GIS化的第一步就是建立坐標系統(tǒng)。

本文舉例說明如何使用 MapInfo?和 Geosoft?Oasis MontajTM兩種常見GIS勘查軟件建立自定義坐標系，主要是如何在該兩個軟件中定義自定義大地基準，并且給出了如何進行坐標變換的步驟，根據(jù)這些步驟可以得到可靠的坐標數(shù)值。

致謝 論文的撰寫過程中得到中國科學(xué)院遙感與數(shù)字地球研究所徐新超博士的寶貴建議，在此表示感謝。

［注釋］

① Pitney Bowes Software Inc.2010.MapInfo professional user guide(edition v.10.5)

② Geosoft Inc.2010.Oasis Montaj v.7.2 mapping and processing system tutorials

Chen Wu-tian.2010.Construction survey in countries using UTM projection［J］.Science and Technology Information，(8):36－37(in Chinese)

Chen Shi-yin.1997.Method to establishment local independent coordinate system［J］.Bulletin of Survey and Mapping，(10):5－7(in Chinese with English abstract)

Chang Kai-shi.2008.Discussion on problems related with the independent coordinate system establishment［J］.City Survey，(1):86－89(in Chinese with English abstract).

Dozier J.1980.Improved algorithm for calculation of UTM and geodetic coordinates［R］.NOAA Technical Report NEES 81:1－19

Kawase K.2011.A general formula for calculating meridian arc length and its application to coordinate conversion in the Gauss-Krüger pro-jection［J］.Bulletin of the Geospatial Information Authority ofJapan，59:1－13

Liu Jian，Liu Gao-Feng.2005.Algorithm of Coordinates Conversion in Gauss-Kruger Projection［J］.Computer Emulation:22(10):119－121(in Chinese with English abstract)

Osborne P.2008.The Mercator Projection［M］.Edinburg，1－189

Shen Ben-zhong.1986.A new calculation for Mercator transverse projection ellipsoidal coordination［J］.Journal of Xi’an College of Geology，8(1):71－86(in Chinese with English abstract)

Xia Lan-fang，Hu Peng，Huang Meng-long.2007.A numerical implementation of analytical transformation using map projection［J］.Science of Surveying and Mapping，32(3):69－71(in Chinese with English abstract)

Xiong Zhong-zhao.2010.Establishment of independent coordinate systems on the UTM projection［J］.Geospatial Information，8(2):41－43(in Chinese with English abstract)

Wang Bao-jun，Song Guo-min，Luo Fen-yong，Mao Yu-zhu，Chen Lingyu.2010.Map Projection Transformation Interface Technology［J］.Journal of Geomatics Science and Technology，27(6):463－466

Yang Qi-he.1981.The Gauss-Krüger projection and the Transverse Mercator projection［J］.Bulletin of Survey and Mapping，(6):34－37(in Chinese)

Yang Qi-he.1986.The research of theory and application of map projection transformation［J］.Journal of PLA Institute of Surveying and Mapping，(1):65－72(in Chinese)

Zhou Chao-xian，F(xiàn)ang Zhi-feng，Yu Cai-hong，Zhang Yun-guo，Gao Ying-bo，Yan Dan-chen，Yang Qiang.2013.UTM projection and Gauss-Krüger projection and their conversion ［J］.Geology and Exploration，49(5):882－889(in Chinese with English abstract)

［附中文參考文獻］

暢開獅.2008.建立城市獨立坐標系相關(guān)問題的探討［J］.城市勘測，(1):86－89

陳士銀.1997.建立地方獨立坐標系的方法［J］.測繪通報，(10):5－7

陳悟天.2010.使用UTM投影坐標系國家的施工測量［J］.科技資訊，(8):36－37

劉 健，劉高峰.2005.高斯－克呂格投影下的坐標變換算法研究［J］.計算機仿真，22(10):119－121

沈本忠.1986.橢球面橫墨卡托投影坐標計算新公式［J］.西安地質(zhì)學(xué)院學(xué)報，8(1):71－86

王寶軍，宋國民，羅奮勇，毛玉柱，陳令羽.2010.地圖投影變換接口技術(shù)［J］.測繪科學(xué)技術(shù)學(xué)報，27(6):463－466

夏蘭芳，胡 鵬，黃夢龍.2007.地圖投影解析變換的數(shù)值實現(xiàn)方法［J］.測繪科學(xué)，32(3):69－71

熊忠招.2010.淺談UTM投影下獨立坐標系統(tǒng)建立［J］.地理空間信息，8(2):41－43

許婭婭，黃文元.2008.山區(qū)公路測量坐標系的選擇方法研究［J］.測繪通報，(7):26－28

楊啟和.1981.高斯－克呂格投影和橫墨卡托投影［J］.測繪通報，(6):34－37

楊啟和.1986.地圖投影變換理論和應(yīng)用研究［J］.解放軍測繪學(xué)院學(xué)報，(1):65－72

周朝憲，房志峰，于彩虹，張云國，高應(yīng)波，燕丹晨，楊強.2013.UTM投影和Gauss－Krüger投影及其變換實現(xiàn)［J］.地質(zhì)與勘探，49(5):882－889