自定义坐标系的建立及其坐标变换实现

周朝宪，董少波，和志军，姜高珍，燕丹晨

(1.中国地质大学，北京 100083;2.南凯矿业有限公司，北京 100023;

3.中色地科矿产勘查股份有限公司，北京 100012;4.中国地质科学院，北京 100037;5.国家海洋环境预报中心，北京 100081)

1 引言

20世纪90年代以来，地质勘查工作的一大进展就是地理信息系统(GIS)越来越多地渗入并支撑着当今的勘查工作。现今的地质勘查工作要求地质勘查人员掌握投影、大地基准(datum)和数学椭球体(ellipsoid)的含义和关系，以免造成错误和不必要的问题。

所谓投影，包括两部分内容:一是对地球这个椭球体参数最佳近似的理论化，即构建大地基准;二是把这一理论化的地球大地基准体上的坐标点(大地坐标系的经纬度)转到平面坐标(即方里网)上(周朝宪等，2013)。大地基准实际也包括两部分内容:首先把实际地球先简单化为数学椭球体。这个数学椭球体由地球半径长轴和扁平率(flattening，扁平率为椭圆率(ellipticity)或者偏心率(eccentricity)的函数)来定义。由于实际地球和数学椭球体还存在比较大的差异(主要是地球质心和数学椭球体并不严格重合)，难以满足测绘和各种工程所需，所以大地基准还包括对数学椭球体进行地球质心的偏移补偿。ellipsoid、datum和实际地球表面的关系见图1。

地球上任一点的经纬度都是建立在datum基础上，而不是建立在ellipsoid之上。任一点的经纬度在不同datum体系下的数值基本不相同，通常相差几十米乃至数百米。了解这一点十分重要，有些地质人员在地质图中没有给出datum，这会给使用图件的人带来很大麻烦。

目前，GPS全球定位系统所采用的UTM(Universal Transverse Mercator，通用横轴墨卡托)投影其datum是WGS(World Geodetic System)84椭球体。理论上讲，投影方法和datum可以自由组合，从而产生众多的坐标系统。另外，由于一些历史原因，在不同的历史时期，许多国家和地区都建立有各自的基准datum和ellipsoid。如20世纪50年代以来，我国使用过与北京54坐标系对应的Krassovsky椭球体和西安80坐标系的IAG1975椭球体。

图1 地球实际地形、大地基准和数学椭球体的关系Fig.1 The relationship among the Earth actual surface，Datum/Geoid and Ellipsoid

实际地质勘查中，如果地质人员不了解这些坐标知识，会给工作带来一些麻烦。如某CXB地质队在埃塞俄比亚工作时，发现手持式GPS给出的坐标和当地的地形/地质图不相吻合，便按照国内处理北京54坐标系和西安80坐标系的办法，根据当地地形图上的控制点，建立了使用WGS 84为椭球体(手持式GPS接收机坐标所依据的椭球体)，并自行修改手持GPS的D x、D y和D z。这实际是建立了独特的TM投影，而地形图上标注着其使用的datum和椭球体分别是Adindan和Clarke_1880。在已知图件的投影和datum的情况下，正确的做法有2种:(1)把地形图通过坐标变换到WGS 84坐标下，以便和手持式GPS接收机数据套合。因为有些手持式GPS接收机往往没有Adindan和Clarke_1880体系;(2)把手持式GPS坐标下的地质成果通过GIS软件转换成Adindan和Clarke_1880坐标系，从而把不同坐标体系下的数据进行套合。

为了使用该地质队自定义坐标体系下的地质成果，可在GIS软件中根据其datum采用WGS 84，以及D x、D y和D z建立了该自定义坐标系，从而利用软件快捷地进行不同坐标体系的转换，解决不同坐标体系下地质数据的套合使用。

周朝宪等(2013)给出了使用MapInfoⓇ对UTM投影和Gauss-Krüger投影下方里网坐标和其经纬度的转换过程。其实，对于国内北京54坐标系和西安80坐标系(都采用Gauss-Krüger投影，UTM投影和Gauss-Krüger投影都属于TM横轴墨卡托投影)，由于该两种坐标系只有数学椭球体而没有datum，必须先在GIS软件中建立自定义坐标系，即按照自定义坐标系的方式处理即可。否则，基本无法利用诸多通用GIS软件对北京54坐标系以及西安80坐标系下的地质数据进行处理。

即使在国内，地质勘查人员也经常面临坐标变换的问题。如北京54坐标系和西安80坐标系之间的相互转换以及和其他投影坐标系(如CGCS2000以及UTM等等)间的变换。因而，建立自定义坐标系以及实现不同坐标系下数据的套合使用已经迫切摆在地质勘查人员面前。

熊忠招(2010)、陈士银(1997)、畅开狮(2008)等从测绘的角度探讨了如何建立个别的自定义坐标系，但其方法不适合于地勘工作;杨启和(1986)、夏兰芳等(2007)、沈本忠(1986)、Osborne(2008)、Dozier(1980)和Kawase(2011)等就坐标变换给出了严格的数学换算，但是这些运算对绝大多数地质勘查人员而言，不具有可操作性;还有人提出了简易的算法公式，但由于简易公式所得结果的精度太差，在实际工作中是不可以接受的(见周朝宪等，2013);刘健等(2005)、王宝军等(2010)发展出一些专门的软件程序用以坐标变换，但是这些软件一则普及性很低，广大地勘人员一般无法使用，二则也面临构建datum和投影等问题，并不适用于广大地勘人员使用。目前，对广大地质人员而言，比较可行的办法就是使用现成的比较通用的GIS软件来进行转换(周朝宪等，2013)，而坐标转换的前提往往就是在这些软件中建立起自定义坐标系。

常用的GIS软件包括 MapInfoⓇ、GeosoftⓇOasis MontajTM、ArcGISⓇ以及国内较为通用的 MapGISⓇ。MapGISⓇ6.7版本(现在最常用的版本)及其以前版本的软件由于其显示的坐标仅仅是图面坐标，还算不上是完整的勘查GIS软件，尽管其描图功能较为强大。完整的勘查GIS软件至少应具有三个特点:(1)图面显示坐标为该点在所采用的GIS坐标体系中的真实坐标数值。当然，其不随着图面比例尺变化而变化。(2)含有目前常见的坐标体系，并可以进行坐标相互转换。(3)基本可以在其中建立自定义坐标系。MapInfoⓇ、GeosoftⓇOasis MontajTM和ArcGISⓇ都具备这三个特点，这些软件在非洲等国家地勘人员中的普及程度也远远高于我国，建议我国地质勘查人员尽快熟悉。

本文使用西方地质勘查应用最广的MapInfoⓇ和GeosoftⓇOasis MontajTM进行展示，如何建立自定义坐标，即如何在这些软件中建立自定义datum，并进行坐标变换。实际上，ArcGISⓇ也可以定义自定义坐标系，也可以使用其处理我国的北京54坐标系和西安80坐标系下的GIS数据。步骤与MapInfoⓇ和GeosoftⓇOasis MontajTM类似，其投影信息存储在数据集的PRJ文本文件中，也是通过直接修改该文本文件来建立自定义坐标系。陈悟天(2010)和周朝宪等(2013)认为，ArcGISⓇ没有国内现用的 Gauss-Krüger投影坐标的说法是不准确的，可以在该软件中建立国内各个地区的datum，然后进行坐标变换。

2 使用MapInfo建立自定义坐标系及其坐标变换

MapInfoⓇ的坐标文件主要是mapinfow.prj，mapinfow.prj是个文本文件，可以直接编辑，其结构参见表1。

MapInfoⓇ自定义投影参数中，对于横轴墨卡托投影坐标(UTM投影和Gauss-Krüger投影都属于横轴墨卡托投影)而言，主要有如下参数:

投影代号(Type)，基准面(datum)，单位(U-nit)，原点经度(Origin Longitude，对 UTM 和 Gauss-Krüger投影而言，其原点经度皆为该带的中央经线的经度值)，原点纬度(Origin Latitude，为赤道纬度值0)，比例因子(Scale Factor，对UTM投影和Gauss-Krüger投影而言，比例因子分别为0.9996和1)，东伪偏移(False Easting)，北伪偏移(False Northing)。关于UTM和Gauss-Krüger投影的上述参数更多内容具体可参见周朝宪等(2013)。

在MapInfoⓇ中，横轴墨卡托投影的代号为8，其余各个投影的代号参见MapInfoⓇ的说明书。各个基准面datum的MapInfoⓇ代号也参见MapInfoⓇ说明书。如，WGS 84 datum的代号是104。

例如:对东经34°和北纬8°的点，其 UTM投影选取WGS 84 world作为大地基准(其数学椭球体是WGS 84)的MapInfoⓇ参数为:

"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带，其起始经线(即中央经线)为33°E;其起始纬度为0°;UTM投影的比例因子是0.9996;东伪偏移为500，000m;北伪偏移为0m。

如上文所述，北京54坐标系和西安80坐标系没有datum，这就要求我们在datum这一项中予以定义。在MapInfoⓇ中，基准的定义主要有数学椭球体与相应的参数组成(即三参数或七参数)。北京54坐标系采用的是Krassovsky椭球体，其在MapInfoⓇ中的代号是3，WGS 84数学椭球体的MapInfoⓇ代号是28。

例如，河南新县北部某地的大地坐标(WGS 84)为 31°46’06.9’’N 和 114°36’41.6’’E，其北京 54坐标系该3度带(带号38，中央经线114°E，三参数)的平移校正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代表起始经线(对横轴墨卡托投影即为该带的中央经线)为114°E;

0代表起始纬线为赤道;

1代表比例因子为1;

500000代表东伪偏移500，000m;

0代表北伪偏移0m。

表1 MapInfoⓇ中坐标系统各个参数定义格式Table 1 The parameter form of coordinate system in Map InfoⓇ

然后，打开MapInfoⓇ软件所在的目录下的mapinfow.prj文件，把上述语句(1)写入，保存，使用MapInfoⓇ软件调用该投影即可。关于不同坐标系的转换过程步骤可参见周朝宪等(2013)。

对上文提到的某CXB地质队在埃塞俄比亚西部Abobo的地区，其大地坐标(WGS 84)为7°39’N和34°39’E，采用的椭球体是WGS 84，属于6度带，该区采用三参数格式:D x=－100，D y=2和D z=－1，该6度带中央经线为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 数学椭球体，－100，2，－1 分别为 D x，D y和D z;

7代表方里网单位为m;

33代表起始经线(对横轴墨卡托投影即为中央经线)为 33°E;

0代表起始纬线为赤道;

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

500 000代表东伪偏移500，000m;

0代表北伪偏移0m。

以上是以三参数为例，对七参数的自定义坐标系，其 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为旋转校正，K为比例校正。这里需要说明的是:在mapinfow.prj文件中添加自定义坐标语句时，语句应尽可能短，否则MapInfoⓇ调用该语句时经常会出现错行，从而引起混乱。

因而，使用MapInfoⓇ软件设立自定义坐标系的步骤归纳为:

(1)确立投影类型及其投影的相关参数(比例因子、起始经纬度和经向和纬向伪偏移量)。投影类型的MapInfoⓇ代号参见MapInfoⓇProfessional User Guide①。

(2)确定datum。确立datum就要确立数学椭球体和datum的相关参数。椭球体的MapInfoⓇ代号参见GeosoftⓇOasis MontajTM软件下的ellipsoid.csv文件。

建立自定义坐标系后，便可以使用MapInfoⓇ实现坐标变换，也可以使用建立自定义坐标系后的GeosoftⓇOasis MontajTM来实现坐标变换，后者用户界面更为友好一些。

3 使用GeosoftⓇOasis Montaj TM建立自定义坐标系

地质勘查工作一般不会建立全新的投影方法，通常是建立适合于各个地区的datum。建立自定义坐标系的核心就是建立适合于工作区的datum。下面简单叙述如何使用GeosoftⓇOasis MontajTM软件建立埃塞俄比亚西部Abobo地区的自定义datum。

首先，修改 GeosoftⓇOasis MontajTM目录下的“datumtrf.csv”，在其中加入如表2内容:

然后，修改 GeosoftⓇOasis MontajTM目录下的“ldatumtrf.csv”，在其中加入表3内容。

因为无论MapInfoⓇ还是GeosoftⓇOasis MontajTM对datum的命名都遵循EPSG(European Petroleum Survey Group)/POSC(Petrotechnical Open Software Corporation)命名规则。自定义的datum不属于其范围，故此前面加上*以示区别，如埃塞俄比亚西部Abobo地区的自定义datum为“*Abobo special”，在表2和表3中该名称必须完全一致。

表2 GeosoftⓇOasis Montaj TM中datumtrf.csv各个参数定义格式②Table 2 The parameter ofcustom coordinate system in datumtrf.csv of GeosoftⓇ Oasis Montaj TM②

表3 GeosoftⓇOasis Montaj TM中自定义坐标ldatum.csv中各个参数定义格式②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

建立自定义坐标系后，可以把坐标数据导入GeosoftⓇOasis MontajTM。第一步:通过图2菜单上的“坐标”选择子菜单“坐标系”，选择我们上文新建的“*Abobo special”，把这些数据设置成以X_CXB和Y_CXB分别为经向和纬向坐标的“*Abobo special”坐标系。第二步:通过图2的菜单上的“坐标”选择“新建投影坐标系”，以X_CXB和Y_CXB分别为当前的经向和纬向坐标(大地基准为“*Abobo special”)，得出 WGS 84 经纬度(即 Long_WGS 84 和Lat_WGS 84)。图2中的经纬度单位是度分秒，如图 2中 Point 19行对应的经度值34.46.15.8375898670 为 34°46’15.8375898670’’。第三步:同第二步类似，以 Long_WGS 84和 Lat_WGS 84为当前坐标数值(大地基准为WGS 84)，选择“新建投影坐标系”，新建UTM WGS 84坐标系，即得出图2中UTM WGS 84经向(X_UTM)和纬向(Y_UTM)坐标。

4 结论

投影坐标包括建立大地基准(datum)和把地球椭球体上的坐标点转到平面坐标(即方里网)上。大地基准(datum)包括数学椭球体(ellipsoid)和根据地球质量对该数学椭球体的补偿。数学椭球体、大地基准和投影方式三者可以组合成多种投影坐标系。地质勘查人员经常会遇到坐标转换问题和建立自定义坐标问题。我国北京54坐标系和西安80坐标系地质数据的GIS化，基本上需要建立自定义坐标。因为其只有数学椭球体，而缺乏大地基准。建立自定义坐标靠手工不可能实现，必须借助软年实现。GIS软件已经成为地质勘查的基本工具，而地质数据的GIS化的第一步就是建立坐标系统。

本文举例说明如何使用 MapInfoⓇ和 GeosoftⓇOasis MontajTM两种常见GIS勘查软件建立自定义坐标系，主要是如何在该两个软件中定义自定义大地基准，并且给出了如何进行坐标变换的步骤，根据这些步骤可以得到可靠的坐标数值。

致谢 论文的撰写过程中得到中国科学院遥感与数字地球研究所徐新超博士的宝贵建议，在此表示感谢。

［注释］

① 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，Fang 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.建立城市独立坐标系相关问题的探讨［J］.城市勘测，(1):86－89

陈士银.1997.建立地方独立坐标系的方法［J］.测绘通报，(10):5－7

陈悟天.2010.使用UTM投影坐标系国家的施工测量［J］.科技资讯，(8):36－37

刘 健，刘高峰.2005.高斯－克吕格投影下的坐标变换算法研究［J］.计算机仿真，22(10):119－121

沈本忠.1986.椭球面横墨卡托投影坐标计算新公式［J］.西安地质学院学报，8(1):71－86

王宝军，宋国民，罗奋勇，毛玉柱，陈令羽.2010.地图投影变换接口技术［J］.测绘科学技术学报，27(6):463－466

夏兰芳，胡 鹏，黄梦龙.2007.地图投影解析变换的数值实现方法［J］.测绘科学，32(3):69－71

熊忠招.2010.浅谈UTM投影下独立坐标系统建立［J］.地理空间信息，8(2):41－43

许娅娅，黄文元.2008.山区公路测量坐标系的选择方法研究［J］.测绘通报，(7):26－28

杨启和.1981.高斯－克吕格投影和横墨卡托投影［J］.测绘通报，(6):34－37

杨启和.1986.地图投影变换理论和应用研究［J］.解放军测绘学院学报，(1):65－72

周朝宪，房志峰，于彩虹，张云国，高应波，燕丹晨，杨强.2013.UTM投影和Gauss－Krüger投影及其变换实现［J］.地质与勘探，49(5):882－889