韓軍 于亞杰 郭瑋

摘?要：充分整合歷史遙感影像數(shù)據(jù)，挖掘其服務(wù)價值對自然資源管理、經(jīng)濟社會發(fā)展具有十分重要的意義。本文以河北省自然資源廳2022—2023年度“全省歷史遙感影像開放”項目為例，對河北省近10年遙感影像成果進行收集分析，利用布爾莎七參數(shù)模型，采用C#語言編制七參數(shù)計算工具，以ArcGIS10.8作為數(shù)據(jù)處理平臺，完成全省原有參心坐標系影像數(shù)據(jù)向2000國家大地坐標系轉(zhuǎn)換，評定其轉(zhuǎn)換精度，為全省歷史遙感影像數(shù)據(jù)融合及應(yīng)用系統(tǒng)建設(shè)統(tǒng)一了數(shù)據(jù)基準。

關(guān)鍵詞：遙感影像；ArcGIS；C#；CGCS2000；坐標轉(zhuǎn)換；布爾莎七參數(shù)

在科技革新與社會效應(yīng)的驅(qū)動下，地理信息數(shù)據(jù)日益豐富，應(yīng)用前景廣泛［1］。為進一步提升測繪地理信息數(shù)據(jù)的服務(wù)價值，促進資源共享，河北省自然資源廳啟動了“全省歷史遙感影像開放”項目。由于歷史原因，歷年遙感影像數(shù)據(jù)坐標系不統(tǒng)一，無法保障整個項目質(zhì)量。本文利用C#語言開發(fā)坐標轉(zhuǎn)換軟件，通過主流的測繪地理信息數(shù)據(jù)處理平臺，實現(xiàn)大量多源遙感影像數(shù)據(jù)高效、快速、便捷地完成坐標轉(zhuǎn)換。

1?坐標轉(zhuǎn)換模型

本項目作業(yè)范圍覆蓋整個河北省，面積約為188800km2，最終成果為CGCS2000下大地坐標（BLH）。數(shù)據(jù)源主要有兩種情況：一是源坐標系均為2000國家大地坐標系，但中央子午線不同，可采用高斯反算實現(xiàn)平面坐標的向大地坐標的轉(zhuǎn)換；二是數(shù)據(jù)源為不同橢球基準，則采用布爾沙七參數(shù)模型實現(xiàn)。該模型為三維模型，在空間直角坐標系中，兩坐標系之間存在嚴密的轉(zhuǎn)換模型，不存在模型誤差和投影變形誤差，適用于橢球面3°及以上的省級及全國范圍控制點坐標轉(zhuǎn)換［24］。

布爾沙七參數(shù)模型包括3個平移參數(shù)（X0，Y0，Z0）、3個旋轉(zhuǎn)參數(shù)（εX、εY、εz，也被稱為3個歐勒角）和1個尺度參數(shù)（m）。本次采用布爾莎模型七參數(shù)法實現(xiàn)原參心坐標系到CGCS2000的轉(zhuǎn)換。如式（1）、式（2）所示［56］：

布爾莎七參數(shù)公式為：

X2Y2Z2=（1+m）X1Y1Z1+0εZ-εY-εZ0εXεY-εX0X1Y1Z1+X0Y0Z0（1）

一般展開式形式是：

X2=X0+（1+m）X1+εZY1-εYZ1Y2=Y0+（1+m）Y1-εZX1+εXZ1Z2=Z0+（1+m）Z1+εYX1-εXY1（2）

上述公式中，X0、Y0、Z0、εX、εY、εZ、m即為七參數(shù)，（X1，Y1，Z1）代表原有參心坐標系成果，（X2，Y2，Z2）代表CGCS2000坐標系下成果。轉(zhuǎn)換參數(shù)求取采用最小二乘法，至少需要3個公共點。當公共點個數(shù)較多時，可根據(jù)測量平差原理列出觀測值的誤差方程式，組成并解算法方程，求得轉(zhuǎn)換參數(shù)并評估精度。

2?已有數(shù)據(jù)資料

2.1?遙感影像數(shù)據(jù)

本項目主要數(shù)據(jù)來源為歷年全省地理國情監(jiān)測、1∶10000基礎(chǔ)測繪更新、全省航空攝影遙感資料獲取與處理等項目數(shù)據(jù)成果，具體情況如下表所示。

2.2?控制點數(shù)據(jù)

本項目選取了河北省自然資源檔案館提供的均勻覆蓋全省的國家A級、B級及省C級GPS點共計238個，作為坐標轉(zhuǎn)換參數(shù)計算依據(jù)。通過附合性檢查，精度良好。

3?技術(shù)路線與方案

3.1?坐標轉(zhuǎn)換流程

河北省橫跨38、39、40共3個分帶，因此涉及不同類型的坐標轉(zhuǎn)換。由于分區(qū)的原因，不同子午線鄰帶接邊處易出現(xiàn)轉(zhuǎn)換盲區(qū)，以分帶38（114°）向39（117°）轉(zhuǎn)換為例，首先對38（114°）的影像進行坐標轉(zhuǎn)換，同時在兩分帶接邊處，選取39（117°）的影像換帶至38（114°），與中央子午線為117°的影像一并進行坐標轉(zhuǎn)換，從而消除接邊漏洞。由于數(shù)據(jù)源存在多種分辨率，接邊處存在不滿幅的問題，需要采用數(shù)據(jù)融合方法解決。以分辨率0.8m、1.0m為例，首先將兩套影像數(shù)據(jù)按本文所述流程分別進行坐標轉(zhuǎn)換，然后將1.0m影像融合至0.8m影像，完成坐標轉(zhuǎn)換，對轉(zhuǎn)換后的影像進行坐標精度、接邊精度等方面檢查，合格后方可開展下一步工序。

針對不同橢球轉(zhuǎn)換，以1980西安坐標系向CGCS2000坐標系轉(zhuǎn)換為例（1954年北京坐標系轉(zhuǎn)換相同），主要流程如下：

（1）利用選取的均勻覆蓋全省的國家、省級高等級GPS點共計238個參與轉(zhuǎn)換參數(shù)計算與檢核，其中選定200個點為求參點，38個點為檢核點。

（2）分別對控制點1980西安坐標系、CGCS2000坐標系數(shù)據(jù)進行高斯反算，得到基于相應(yīng)橢球下的大地坐標BLH，通過數(shù)學模型將轉(zhuǎn)換至空間直角坐標XYZ，然后采用布爾沙七參數(shù)模型，基于最小二乘法求取轉(zhuǎn)換參數(shù)，并計算重合點坐標殘差。

（3）利用計算得到的七參數(shù)，將待轉(zhuǎn)換數(shù)據(jù)全部轉(zhuǎn)換至CGCS2000坐標系下。

3.2?采用C#開發(fā)坐標轉(zhuǎn)換軟件

利用BursaWolf七參數(shù)轉(zhuǎn)換模型，通過QR分解解決了矩陣求逆過程中出現(xiàn)的數(shù)值不穩(wěn)定問題，按照最小二乘法基于Visual?Studio?2015平臺實現(xiàn)坐標系之間的轉(zhuǎn)換。軟件界面如下圖所示。

主要代碼如下：

private?void?button3_Click（object?sender，EventArgs?e）

{

double［，］S=new?double［7，1］；

if（newX?！=null）

{

Matrix?Matrix=new?Matrix（）；

double［，］B=new?double［3*length，7］；for（int?i=0；i<3*length；i++）for（int?j=0；j<7；j++）B［i，j］=0；

double［，］l=new?double［3*length，1］；

for（int?i=0；i

{

B［3*i，0］=1；B［3*i+1，1］=1；B［3*i+2，2］=1；

B［3*i，3］=oldX［i，0］；B［3*i+1，3］=oldX［i，1］；B［3*i+2，3］=oldX［i，2］；

B［3*i，5］=-oldX［i，2］；B［3*i+1，6］=-oldX［i，0］；B［3*i+2，4］=-oldX［i，1］；

B［3*i，6］=oldX［i，1］；B［3*i+1，4］=oldX［i，2］；B［3*i+2，5］=oldX［i，0］；

l［3*i，0］=newX［i，0］；l［3*i+1，0］=newX［i，1］；l［3*i+2，0］=newX［i，2］；

}

para=Matrix.MultiplyMatrix（Matrix.Athwart（Matrix.MultiplyMatrix（Matrix.Transpose（B），B）），Matrix.MultiplyMatrix（Matrix.Transpose（B），l））；S=para；

xtextBox.Text=para［0，0］.ToString（）；ytextBox.Text=para［1，0］.ToString（）；ztextBox.Text=para［2，0］.ToString（）；

mtextBox.Text=Convert.ToString（（para［3，0］-1）*1000000）；

extextBox.Text=Convert.ToString（para［4，0］/para［3，0］*206265）；

eytextBox.Text=Convert.ToString（para［5，0］/para［3，0］*206265）；

eztextBox.Text=Convert.ToString（para［6，0］/para［3，0］*206265）；

}

本軟件支持空間直角坐標系與大地坐標系互換，通過源坐標、目標坐標的空間直角坐標系計算布爾沙七參數(shù)，可實現(xiàn)坐標文件的批量轉(zhuǎn)換。

3.3?基于ArcGIS平臺實現(xiàn)數(shù)據(jù)坐標轉(zhuǎn)換

坐標投影轉(zhuǎn)換坐標投影轉(zhuǎn)換實際上就是將源坐標投影信息刪除，然后重新給其定義與目標投影體系相同的坐標系統(tǒng)［7］。首先，啟動ArcGIS?Desktop平臺，添加待轉(zhuǎn)換遙感影像數(shù)據(jù)，由于影像數(shù)據(jù)量較大，因此需要選擇創(chuàng)建金字塔。其次，在Arc?Toolbox中使用數(shù)據(jù)管理工具的投影和變換工具創(chuàng)建自定義地理轉(zhuǎn)換（Create?Custom?Geographic?Transformation），設(shè)置地理變換名稱（Geographic?Transformation?Name），依次輸入源地理坐標系（Import?Geographic?Coordinate?System），輸出目標地理坐標系（Output?Geographic?Coordinate?System）和自定義地理變換方法（Custom?Geographic?Transformation），設(shè)置3.2中計算出的七參數(shù)。最后，選擇新建的地理（坐標）轉(zhuǎn)換方式，進行投影轉(zhuǎn)換。

4?精度要求與評價

4.1?精度要求

本項目最終目的是將遙感影像數(shù)據(jù)制作電子地圖，通過瓦片形式實現(xiàn)公眾服務(wù)。依據(jù)規(guī)范要求，遙感影像平面位置中誤差規(guī)定為：平地、丘陵地區(qū)不大于0.5mm（圖上），山地、高山地不大于0.75mm（圖上），以明顯地物點平面位置中誤差的2倍為其最大誤差。影像拼接數(shù)字影像圖應(yīng)與相鄰影像圖接邊，接邊誤差不大于2個像元。

4.2?精度評價

坐標轉(zhuǎn)換精度采用內(nèi)符合和外符合精度評價，依據(jù)計算轉(zhuǎn)換參數(shù)的重合點殘差中的誤差評估坐標轉(zhuǎn)換精度，殘差小于3倍點位中誤差的點位精度滿足要求。經(jīng)計算，本項目內(nèi)符合轉(zhuǎn)換各點殘差最小為：±0.002m，最大為±0.091m，內(nèi)符合中誤差為±0.032m；外符合轉(zhuǎn)換各點殘差最小為：±0.012m，最大為±0.131m，外符合中誤差為±0.058m，完全滿足本項目精度要求。

將新得到的數(shù)據(jù)與目標數(shù)據(jù)、全省行政界線、基礎(chǔ)線劃圖等參考圖層進行疊加，發(fā)現(xiàn)兩者數(shù)據(jù)吻合情況良好，沒有發(fā)現(xiàn)變形與偏移；與源數(shù)據(jù)相比，屬性信息沒有發(fā)生任何丟失，滿足了質(zhì)量要求。

結(jié)語

本文通過結(jié)合項目實例，采用ArcGIS與C#相結(jié)合，實現(xiàn)了省級歷史遙感影像數(shù)據(jù)坐標轉(zhuǎn)換的方法和流程，建立了嚴密、統(tǒng)一的坐標轉(zhuǎn)換關(guān)系模型，在保證成果質(zhì)量的情況下，提高了工作效率。經(jīng)檢驗，本次工作所采取的轉(zhuǎn)換方法與流程切實可行，對大范圍地理信息數(shù)據(jù)坐標轉(zhuǎn)換問題提供了良好的解決方案。

參考文獻：

［1］賈繼鵬，厲芳婷，侯愛羚.基于多源數(shù)據(jù)融合的基礎(chǔ)地理信息數(shù)據(jù)動態(tài)更新［J］.地理空間信息，2021，44（01）：4143+4.

［2］徐仕琪，張曉帆，周可法，等.關(guān)于利用七參數(shù)法進行WGS84和BJ54坐標轉(zhuǎn)換問題的探討［J］.測繪與空間地理信息，2007，30（5）：3338.

［3］謝艷玲，夏正清.基于ArcGIS的1980西安坐標系到2000國家坐標系的轉(zhuǎn)換研究［J］.測繪與空間地理信息，2014，37（12）：220224.

［4］成英燕，程鵬飛，秘金鐘，等.大尺度空間域下1980西安坐標系與WGS84坐標系轉(zhuǎn)換方法研究［J］.測繪通報，2007（12）：58.

［5］孔祥元，郭際明，劉宗泉.大地測量學基礎(chǔ)［M］.武漢：武漢大學出版社，2010.

［6］郭英起，唐彬，張秋江，等.基于空間直角坐標系的高精度坐標轉(zhuǎn)換方法研究［J］.大地測量與地球動力學，2012，32（3）：125128.

［7］譚玲，秦龍，毛莉莉.基于ArcGIS的DOM坐標系統(tǒng)轉(zhuǎn)換［J］.北京測繪，2017（3）：111113.

基金項目：本論文為河北省自然資源廳2023年度財政項目，立項編號：13000023P006CA410209K

作者簡介：韓軍（1976—?），男，漢族，河北石家莊人，本科，高級工程師，航空攝影室主任，主要從事測繪地理信息數(shù)據(jù)采集處理及系統(tǒng)研發(fā)工作。