鈾礦地質(zhì)成果數(shù)據(jù)2000國家大地坐標(biāo)系轉(zhuǎn)換方法研究

徐志強(qiáng)，王貴，王龍輝

（核工業(yè)二〇八大隊(duì)，內(nèi)蒙古 包頭 014010）

隨著我國社會(huì)的發(fā)展，以往參心坐標(biāo)系逐漸不能滿足科學(xué)研究、經(jīng)濟(jì)建設(shè)、國防軍事及社會(huì)需要。地心坐標(biāo)系的使用有利于現(xiàn)代空間技術(shù)維護(hù)和坐標(biāo)系的快速更新，高精度控制點(diǎn)三維坐標(biāo)的測(cè)量，提高工作效率。一些測(cè)繪前沿高新技術(shù)，如網(wǎng)絡(luò)RTK 技術(shù)、北斗GPS 系統(tǒng)等的建立都是使用2000 國家大地坐標(biāo)系統(tǒng)（CGCS2000）。各行業(yè)對(duì)2000 國家大地坐標(biāo)系的運(yùn)用是必然趨勢(shì)，因此參心坐標(biāo)系的行業(yè)資料和空間成果數(shù)據(jù)轉(zhuǎn)換為2000地心坐標(biāo)系就極為重要。成熟空間技術(shù)的廣泛應(yīng)用，需要高度精確的、變化發(fā)展的、實(shí)時(shí)的、統(tǒng)一的大地坐標(biāo)系統(tǒng)，作為鈾礦地質(zhì)工作者各項(xiàng)成果表達(dá)的基礎(chǔ)性保障［1］，地心坐標(biāo)系在鈾礦地質(zhì)行業(yè)的應(yīng)用已經(jīng)勢(shì)在必行。經(jīng)國務(wù)院批準(zhǔn)我國自2008 年7 月1 日啟用2000 國家大地坐標(biāo)系［2］，2018 年國家地質(zhì)行業(yè)要求全面實(shí)行2000 坐標(biāo)系。鈾礦地質(zhì)行業(yè)當(dāng)前的坐標(biāo)轉(zhuǎn)換主要是用簡單的測(cè)量軟件進(jìn)行單點(diǎn)轉(zhuǎn)換，文件輸入輸出格式繁瑣，圖件轉(zhuǎn)換主要是手工轉(zhuǎn)換為主，工作效率很低。

經(jīng)過數(shù)十年的積累，鈾礦地質(zhì)成果數(shù)量龐大。僅鄂爾多斯盆地北部的數(shù)據(jù)，電子資料也有上百個(gè)GB。如果用簡單的測(cè)量軟件對(duì)這些數(shù)據(jù)進(jìn)行單點(diǎn)轉(zhuǎn)換，再手工轉(zhuǎn)換圖件，其工作量和花費(fèi)的工作成本將是龐大的數(shù)字。因此需要我們找到一種快速、實(shí)用、易學(xué)、準(zhǔn)確的方法，并推廣應(yīng)用，為鈾礦地質(zhì)成果數(shù)據(jù)坐標(biāo)系的轉(zhuǎn)換提供技術(shù)支持。

1 鈾礦地質(zhì)成果資料分析

通過對(duì)鄂爾多斯盆地北部成果資料進(jìn)行梳理，截止2018 年底，鄂爾多斯盆地北部目前資料主要包括：區(qū)域地質(zhì)資料、航測(cè)遙感資料、地面物化探資料、鈾成礦環(huán)境綜合研究與編圖資料、鈾資源調(diào)查評(píng)價(jià)資料、鈾礦勘查資料、專題研究資料。鈾礦地質(zhì)成果數(shù)據(jù)的電子資料有約140 GB。這些資料90%以上均為MapGIS 制作，實(shí)行2000 國家大地坐標(biāo)系后，80 坐標(biāo)系和54 坐標(biāo)系的資料均需要進(jìn)行坐標(biāo)轉(zhuǎn)換，由于資料涉及不同單位，不同坐標(biāo)系，不同比例尺，種類繁多，且數(shù)量巨大，轉(zhuǎn)換難度較大，如何高效、準(zhǔn)確的進(jìn)行坐標(biāo)轉(zhuǎn)換，成為目前迫在眉睫的任務(wù)。

從空間數(shù)據(jù)的角度分析鈾礦成果資料，可分為兩大類。一類為坐標(biāo)點(diǎn)數(shù)據(jù)，以數(shù)字形式存在于成果資料的報(bào)告和附表中，一般為大地經(jīng)緯度或者是高斯坐標(biāo)投影下的X、Y、H 坐標(biāo)。坐標(biāo)點(diǎn)數(shù)據(jù)的整理方法是把每一個(gè)鈾礦地質(zhì)項(xiàng)目資料中的坐標(biāo)點(diǎn)數(shù)據(jù)整理成一個(gè)EXCEL 數(shù)據(jù)文件，以方便坐標(biāo)點(diǎn)數(shù)據(jù)的坐標(biāo)轉(zhuǎn)換。第二類是圖件，鈾礦地質(zhì)成果有很多都是以矢量圖件形式表示的，如地質(zhì)圖、砂體分布圖、剖面圖等，并且其大部分是MapGIS 數(shù)據(jù)格式。對(duì)于圖件的整理方法是把圖件按地區(qū)和不同比例尺分類整理，并且按項(xiàng)目放入文件夾內(nèi)，方便圖件的批量坐標(biāo)轉(zhuǎn)換工作。

2 坐標(biāo)轉(zhuǎn)換模型介紹及確定

2.1 常用轉(zhuǎn)換模型介紹

參心坐標(biāo)向地心坐標(biāo)的轉(zhuǎn)換，實(shí)際上是不同測(cè)量坐標(biāo)基準(zhǔn)下的坐標(biāo)轉(zhuǎn)換［3］，其實(shí)質(zhì)上是兩個(gè)空間直角坐標(biāo)系之間的數(shù)學(xué)轉(zhuǎn)換關(guān)系。如果建立了兩個(gè)空間坐標(biāo)系之間的數(shù)學(xué)模型，并且研究出了數(shù)學(xué)算法，就可以實(shí)現(xiàn)兩個(gè)坐標(biāo)系之間的坐標(biāo)轉(zhuǎn)換。測(cè)繪行業(yè)常用的數(shù)學(xué)模型有三維布爾莎模型、二維七參數(shù)轉(zhuǎn)換模型、平面四參數(shù)轉(zhuǎn)換模型［4］。

2.1.1三維七參數(shù)坐標(biāo)轉(zhuǎn)換模型（布爾莎模型）

式中：ΔB—同一點(diǎn)為在兩個(gè)坐標(biāo)系下的緯度差/rad；ΔL—同一點(diǎn)為在兩個(gè)坐標(biāo)系下的經(jīng)度差/rad；ΔH—同一點(diǎn)為在兩個(gè)坐標(biāo)系下的大地高差/m［5］；ρ″—180×3600/（rad·s-1）；Δa—橢球長半軸差/m；Δf—扁率差；ΔX、ΔY、ΔZ—平移參數(shù)/m；εx、εy、εz—旋轉(zhuǎn)參數(shù)/rad；m—尺度參數(shù)；e—第一篇心率；N—卯西圈曲率半徑；M—子午圈曲率半徑。

2.1.2二維七參數(shù)轉(zhuǎn)換模型

式中：ΔB—公共點(diǎn)在兩個(gè)大地坐標(biāo)系的緯度差/rad；ΔL—公共點(diǎn)在兩個(gè)大地坐標(biāo)系的經(jīng)度差/rad；Δa—橢球長半軸差/m；Δf—橢球扁率差；ΔX、ΔY、ΔZ—平移參數(shù)/m；εx、εy、εz—旋轉(zhuǎn)參數(shù)/rad；m—尺度參數(shù)［5］；ρ″—180×3600/（rad·s-1）；e—第一篇心率；N—卯西圈曲率半徑；M—子午圈曲率半徑。

2.1.3平面四參數(shù)轉(zhuǎn)換模型

平面四參數(shù)模型是二維坐標(biāo)轉(zhuǎn)換，其通過最小二乘數(shù)學(xué)算法實(shí)現(xiàn)。平面直角坐標(biāo)轉(zhuǎn)換模型［5］：

式中：x0、y0—平面坐標(biāo)平移參數(shù)/m；α—旋轉(zhuǎn)參數(shù)；m—尺度參數(shù)；x2、y2—目標(biāo)大地坐標(biāo)系中的平面直角坐標(biāo)/m；x1、y1—原始坐標(biāo)系中的平面直角坐標(biāo)/m。

2.2 轉(zhuǎn)換模型的確定

以上3 個(gè)模型各有不同，三維模型精度相對(duì)比較高。高精度測(cè)繪對(duì)高程精度要求很高，一般采用三維七參數(shù)模型。平面四參數(shù)法考慮了兩個(gè)平面移參數(shù)ΔX、ΔY 和一個(gè)旋轉(zhuǎn)參數(shù)α 及一個(gè)比例參數(shù)K，能夠在一定的平面范圍內(nèi)高精度的完成坐標(biāo)轉(zhuǎn)換，在公共點(diǎn)個(gè)數(shù)為3～5 時(shí)，模型轉(zhuǎn)換就可以達(dá)到較高的精度，并且在整個(gè)平面參數(shù)轉(zhuǎn)換過程中高程不參與。對(duì)高程要求不高的行業(yè)，如地質(zhì)行業(yè)，平面四參數(shù)法是一個(gè)方便快捷的方法。

3 坐標(biāo)轉(zhuǎn)換技術(shù)流程和轉(zhuǎn)換模型

3.1 技術(shù)流程

大部分鈾礦成果圖件對(duì)平面精度要求較高，基本上在厘米級(jí)，而部分對(duì)高程要求相對(duì)要低，分米級(jí)即可以滿足生產(chǎn)需要，結(jié)合對(duì)坐標(biāo)轉(zhuǎn)換不同模型的特點(diǎn)分析，確定采用基于最小二乘法的平面四參數(shù)模型，實(shí)現(xiàn)鈾礦勘查成果空間數(shù)據(jù)的2000 國家坐標(biāo)系的轉(zhuǎn)換工作。

第一步：整理分析單點(diǎn)空間資料，如鉆探工程點(diǎn)坐標(biāo)、取樣位置坐標(biāo)、礦體拐點(diǎn)坐標(biāo)、礦權(quán)拐點(diǎn)等坐標(biāo)資料，成果圖件按不同的比例尺分類；

第二步：再到相應(yīng)的政府部門收集控制點(diǎn)資料，控制點(diǎn)要在整個(gè)區(qū)域均勻分布且有足夠數(shù)量，確保有原坐標(biāo)系與2000 坐標(biāo)系的重合點(diǎn)，同時(shí)收集一定數(shù)量的用于檢查點(diǎn)的控制點(diǎn)；

第三步：基于EXCEL 平臺(tái)下，用VBA編寫的平面四參數(shù)轉(zhuǎn)換程序計(jì)算出平面四參數(shù)，再用VBA 編寫的程序把原坐標(biāo)系的單點(diǎn)坐標(biāo)通過平面四參數(shù)轉(zhuǎn)換為2000 國家大地坐標(biāo)；

第四步：用程序?qū)⒉煌壤叩脑鴺?biāo)系下的矢量圖件，結(jié)合平面四參數(shù)及MapGIS二次開發(fā)包轉(zhuǎn)換為2000 國家大地坐標(biāo)系下的矢量圖件（圖1）。

圖1 技術(shù)流程示意圖Fig.1 Schematic diagram of technical process

3.2 基于EXCEL 平臺(tái)下的平面四參數(shù)模型的軟件實(shí)現(xiàn)

在EXCEL 平臺(tái)下利用VBA 編寫了一個(gè)四參數(shù)計(jì)算的軟件，可以把EXCEL 和四參數(shù)模型結(jié)合起來，使得軟件的操作更加方便快捷。VBA 程序用到數(shù)組運(yùn)算。主要程序代碼及界面如下：

3.3 MapGIS 圖件四參數(shù)的整體轉(zhuǎn)換的程序?qū)崿F(xiàn)

為提高圖件坐標(biāo)轉(zhuǎn)換的效率，在基于MapGIS SDK 軟件開發(fā)包的基礎(chǔ)上進(jìn)行二次開發(fā)，用EXCEL VBA 開發(fā)語言編寫函數(shù)，直接把圖件進(jìn)行坐標(biāo)轉(zhuǎn)換。矢量圖形坐標(biāo)轉(zhuǎn)換調(diào)用MapGIS SDK 類庫MapGIS BasCom1 里的旋轉(zhuǎn)角度函數(shù)、比例尺度函數(shù)以及x、y 平移量函數(shù)，完成圖件矢量數(shù)據(jù)在平面四參數(shù)下的2000 國家大地坐標(biāo)轉(zhuǎn)換。程序主要代碼如下:

4 坐標(biāo)轉(zhuǎn)換方法的實(shí)例應(yīng)用及精度分析

4.1 實(shí)例應(yīng)用

本文研究的方法已經(jīng)用于實(shí)際生產(chǎn)中，現(xiàn)以某鈾礦整裝勘查項(xiàng)目為例，用本文提到的方法對(duì)鈾礦項(xiàng)目的成果資料進(jìn)行2000 國家大地坐標(biāo)進(jìn)行轉(zhuǎn)換。

4.1.1平面四參數(shù)計(jì)算控制點(diǎn)、檢查控制點(diǎn)及圖件內(nèi)檢查點(diǎn)的選取

1）平面四參數(shù)計(jì)算控制點(diǎn)的選取

工作區(qū)面積不要超過30 km2（個(gè)別小比例尺圖件可以更大一些），國家三角點(diǎn)均勻分布在工作區(qū)內(nèi)，最好是4 個(gè)方向都有，并且控制面積比工作區(qū)面積大一些。觀察個(gè)別誤差比較大的三角點(diǎn)，不能參與平面四參數(shù)的運(yùn)算。

2）檢查控制點(diǎn)的選取

為評(píng)定平面四參數(shù)的精度，需要選取一定數(shù)量的控制點(diǎn)檢查平面四參數(shù)轉(zhuǎn)換的準(zhǔn)確性，用檢查點(diǎn)的原坐標(biāo)通過四參數(shù)計(jì)算出2000 國家坐標(biāo)與檢查點(diǎn)已知的2000 國家坐標(biāo)進(jìn)行對(duì)比分析。檢查控制點(diǎn)在工作區(qū)范圍內(nèi)要均勻分布，盡量在工作區(qū)的中部和4 角都有檢查點(diǎn)。

3）圖件內(nèi)檢查點(diǎn)的選取

為確定圖件坐標(biāo)轉(zhuǎn)換精度是否滿足規(guī)范要求，需要在每一個(gè)圖件范圍內(nèi)選取5 個(gè)以上的檢查點(diǎn)來評(píng)定圖件的坐標(biāo)轉(zhuǎn)換精度。圖件控制點(diǎn)在圖件范圍內(nèi)要均勻分布，盡量在工作區(qū)的中部和4 角都有檢查點(diǎn)。

4.1.2平面四參數(shù)計(jì)算及項(xiàng)目單點(diǎn)坐標(biāo)的計(jì)算

根據(jù)選取好的平面四參數(shù)計(jì)算控制點(diǎn)，把控制點(diǎn)的舊坐標(biāo)和2000 國家大地坐標(biāo)輸入到EXCEL VAB 程序相應(yīng)的位置上，就可以算出工作區(qū)內(nèi)的平面四參數(shù)。再把項(xiàng)目成果涉及到的所有單點(diǎn)資料和檢查控制點(diǎn)舊坐標(biāo)及圖件內(nèi)檢查點(diǎn)舊坐標(biāo)輸入程序相應(yīng)的位置上，程序可以計(jì)算出所有的單點(diǎn)坐標(biāo)對(duì)應(yīng)的2000國家大地坐標(biāo)。圖2 為實(shí)例計(jì)算程序界面示意圖。

圖2 平面四參數(shù)計(jì)算程序示意圖Fig.2 Schematic diagram of four planar parameter calculation program

4.1.3圖件的實(shí)例坐標(biāo)轉(zhuǎn)換

選取鄂爾多斯市巴音青格利鈾礦地質(zhì)普查項(xiàng)目成果數(shù)據(jù)實(shí)例分析，項(xiàng)目成果圖件共有MapGIS 的圖件66 張，按比例尺不同分為4類，分別是1:25 萬、1:10 萬、1:5 萬、1:1 萬，本文選取一張1:1 萬的圖件為例進(jìn)行坐標(biāo)轉(zhuǎn)換。

首先，對(duì)圖面現(xiàn)有的圖件進(jìn)行坐標(biāo)檢查，在MapGIS 里面設(shè)好圖形比例尺，使圖上坐標(biāo)和實(shí)際坐標(biāo)一致，并且讀取圖上的工程點(diǎn)坐標(biāo)與工程點(diǎn)的實(shí)際坐標(biāo)對(duì)比是不是在限差范圍內(nèi)；其次，對(duì)圖框進(jìn)行檢查，確定現(xiàn)有圖框是否與舊坐標(biāo)系統(tǒng)一致，是否和圖面的坐標(biāo)系一致；然后，用EXCEL VBA 程序打開選擇需要轉(zhuǎn)換圖件的點(diǎn)、線、面文件，對(duì)其進(jìn)行圖形的2000 國家大地坐標(biāo)轉(zhuǎn)換。圖3 為圖件轉(zhuǎn)換的程序界面。最后，刪除舊坐標(biāo)系圖框，按圖的范圍用MapGIS 軟件生成2000 坐標(biāo)圖框并且與轉(zhuǎn)換后的2000 圖形套合。圖4為2000 坐標(biāo)圖框生成界面。

圖3 批量MapGIS 圖件轉(zhuǎn)換程序示意圖Fig.3 Schematic interface of batch MapGIS map conversion program

4.2 實(shí)例坐標(biāo)轉(zhuǎn)換精度分析

4.2.1實(shí)例平面四參數(shù)模型精度分析（單點(diǎn)轉(zhuǎn)換精度分析）

選取5 個(gè)檢查控制點(diǎn)，用此方法轉(zhuǎn)換后得到的2000 坐標(biāo)與收集的2000 坐標(biāo)進(jìn)行精度對(duì)比分析，較差為5 mm 左右，滿足規(guī)范要求（表1）。

4.2.2實(shí)例圖形轉(zhuǎn)換精度分析

讀取轉(zhuǎn)換后的圖形的檢查點(diǎn)2000 坐標(biāo)與單點(diǎn)計(jì)算的檢查點(diǎn)2000 坐標(biāo)進(jìn)行誤差較對(duì)，并轉(zhuǎn)換為圖面精度，與圖面限差比較，能夠滿足 《國土資源數(shù)據(jù)2000 坐標(biāo)系轉(zhuǎn)換技術(shù)要求》對(duì)圖形空間數(shù)據(jù)轉(zhuǎn)換精度要求，精度應(yīng)小于圖上的0.1 mm（表2）。

圖4 生成2000 坐標(biāo)系圖框示意圖Fig.4 Schematic frame interface of generating CGCS2000 coordinater

表1 檢查控制點(diǎn)程序轉(zhuǎn)換2000 坐標(biāo)與收集2000 國家坐標(biāo)對(duì)照表Table 1 Check and control points of conversing CGCS2000 and collecting 2000 national coordinate

表2 圖形檢查點(diǎn)四參數(shù)轉(zhuǎn)換結(jié)果與圖形轉(zhuǎn)換結(jié)果對(duì)照表Table 2 Comparison of four-parameter conversion results and graph conversion results

5 結(jié)論

當(dāng)前鈾礦地質(zhì)行業(yè)坐標(biāo)轉(zhuǎn)換主要是用簡單的測(cè)量軟件進(jìn)行單點(diǎn)轉(zhuǎn)換，文件輸入輸出格式繁瑣，圖件轉(zhuǎn)換主要是手工轉(zhuǎn)換為主，工作效率很低。本文在EXCEL 軟件平臺(tái)下，使用VBA 編寫空間數(shù)據(jù)平面四參數(shù)轉(zhuǎn)換程序，采用EXCEL 結(jié)合VBA 在MapGIS 二次開發(fā)平臺(tái)上編寫程序，實(shí)現(xiàn)對(duì)MapGIS 格式圖件的2000 國家坐標(biāo)系的批量轉(zhuǎn)換：

1）實(shí)例證明平面四參數(shù)坐標(biāo)轉(zhuǎn)換模型能夠滿足鈾礦地質(zhì)項(xiàng)目2000 國家坐標(biāo)轉(zhuǎn)換的要求。

2）本文在EXCEL 平臺(tái)下開發(fā)的程序，通過實(shí)例驗(yàn)證在生產(chǎn)中精度可靠且容易操作，可以快速準(zhǔn)確地完成國家鈾礦地質(zhì)成果空間數(shù)據(jù)的2000 國家大地坐標(biāo)系轉(zhuǎn)換。該方法在大量鈾礦地質(zhì)成果數(shù)據(jù)的坐標(biāo)轉(zhuǎn)換中將發(fā)揮積極的作用，值得進(jìn)一步推廣。