基于ArcGIS與C#語言的省級歷史遙感影像數據坐標轉換應用研究

韓軍 于亞杰 郭瑋

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

關鍵詞：遙感影像；ArcGIS；C#；CGCS2000；坐標轉換；布爾莎七參數

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

1?坐標轉換模型

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

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

布爾莎七參數公式為：

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即為七參數，（X1，Y1，Z1）代表原有參心坐標系成果，（X2，Y2，Z2）代表CGCS2000坐標系下成果。轉換參數求取采用最小二乘法，至少需要3個公共點。當公共點個數較多時，可根據測量平差原理列出觀測值的誤差方程式，組成并解算法方程，求得轉換參數并評估精度。

2?已有數據資料

2.1?遙感影像數據

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

2.2?控制點數據

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

3?技術路線與方案

3.1?坐標轉換流程

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

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

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

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

（3）利用計算得到的七參數，將待轉換數據全部轉換至CGCS2000坐標系下。

3.2?采用C#開發坐標轉換軟件

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

主要代碼如下：

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）；

}

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

3.3?基于ArcGIS平臺實現數據坐標轉換

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

4?精度要求與評價

4.1?精度要求

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

4.2?精度評價

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

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

結語

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

參考文獻：

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

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

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

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

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

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

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

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

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