探析基于Excel實現工程測量中的坐標轉換

羅孝楠 周廣勇 謝全遠

摘要：因各個國家衛星導航使用的坐標系不一致，GPS工程測量中坐標轉換是實際應用的重要步驟，針對目前坐標轉換軟件數據錄入不便、過程不可控等問題，文中探析四參數法坐標轉換參數求解方法，并在Excel下利用少量內部函數，設計實現了四參數求解和坐標轉換電子表格。通過工程實例坐標進行驗證，求解結果準確，計算過程數據可視，編輯錄入數據方便，方便在測量工程、教學實驗中使用。

關鍵詞：衛星導航；四參數；坐標轉換；Excel

中圖分類號：TP311.1? ? 文獻標識碼：A

文章編號：1009-3044（2023）35-0118-03

開放科學（資源服務）標識碼（OSID）

0 引言

伴隨著衛星技術深入發展，在實踐中得到了廣泛應用，特別是美國全球定位系統（GPS） 的成功建立和應用，我國北斗衛星導航測量技術也取得了前所未有的進步和發展。衛星定位測量技術具有精準、快速、高效和全天候等顯著的特點，使得測繪行業經歷了一場深刻的變革。豐富的應用場景已經被應用到了現代社會經濟中的很多領域，尤其是在軍事、交通、測量等方面和領域取得了顯著的成效，凸顯了重要的應用價值。在因此全球導航定位衛星系統，在全球范圍內的眾多行業中，尤其是在測繪行業有著越來越廣泛的應用。

全球定位系統（GPS） 所采用的坐標系是WGS-84坐標系，因此直接采集的點位坐標屬于WGS-84坐標系，而目前我們國家工程測量成果普遍使用的是以CGCS2000坐標系或是地方（任意）獨立坐標系為基礎的坐標數據。因此必須將WGS-84坐標轉換到CSCS2000坐標系或地方（任意）獨立坐標系[1]。坐標轉換是從一種坐標系統映射到另一種坐標系統的過程。通過建立兩個坐標系統之間一一對應關系來實現。一般四參數應用于平面轉換和七參數應用于三維轉換。

在許多測量數據處理軟件中，如徠卡LGO、天寶TBC、南方SGO、Pinncle等數據處理軟件，都有坐標系轉換模塊，有些功能比較齊全，如在TBC軟件中包含了七參數法、格網法、多元回歸法；LGO軟件中有格網法、七參數法、三參數法、格網與參法結合法，有三維轉換也有二維轉換。在實際應用中，可以結合測區大小、范圍內已知點的數量與分布情況決定采用哪一種方法。這些軟件在學習、教學過程中還是存在一些問題，軟件比較復雜龐大、學習成本比較高，計算過程不可見，不利于對原理的理解。

Excel是微軟公司推出的一個廣泛應用于商業和科學領域的電子表格程序。它具有一些強大的功能和特點，使其成為各種應用程序中不可缺少的工具。簡潔的界面、優秀的計算功能和圖表工具，再加上出色的市場營銷，使Excel成為流行的個人計算機電子表格軟件。它的主要功能有制作表格、數字格式的轉換、制作圖表、函數功能、數據處理功能等[2]。強大的數據編輯功能，使得利用Excel編制坐標轉換程序表格方便可行。

1 四參數坐標轉換原理

四參數是用于兩個平面直角坐標系之間的互相轉換，而七參數是用于兩個三維空間直角坐標系之間的轉換。四參數可以由兩個及以上具有平面二維坐標的已知等級控制點求出，高程可以用擬合方法求出，求解較為簡單，也較容易理解；七參數轉換是空間坐標系的轉換需要在測區布設一定密度的控制網點，需要三個以上已知三維坐標的控制點，利用整個網的WGS-84坐標系下的三維約束平差結果和當地坐標系統的二維約束平差結果及各點的高程解算，參數求解比較煩瑣，理解起來相對困難。

1.1 平面坐標轉換

點位的平面坐標可以用高斯投影取得。當地面兩點的距離小于十千米時，地球曲率的影響可以忽略，此時可以采用四參數來完成兩種坐標系的轉換，精度也能達到預定的需求。四參數轉換計算如公式（1） 。

[x2y2=?x?y+mcosα-sinαsinαcosαx1y1]? ? （1）

其中，x1、y1為原始坐標，x2、y2為目標坐標，單位米，Δx、Δy為兩個坐標平移參數，即兩個平面坐標系的坐標原點之間的坐標差值。α為兩個坐標系的坐標軸間的旋轉角度，通過旋轉α，可以使兩個坐標系的X和Y軸重合在一起。M為尺度縮放因子，即兩個坐標系內的同一段直線的長度比值，實現尺度的比例轉換。通常m值非常接近于1。因為有四個參數所以稱之為四參數坐標轉換。這樣通過坐標的平移、坐標軸的旋轉和尺度的縮放，可以實現平面坐標的轉換。而這四個參數如何求解是坐標轉換的關鍵步驟。

1.2 四參數的求解

四參數求解是坐標轉換的一個逆過程，求解需要至少兩對已知點，分別知道其在目標平面的坐標和原始平面的坐標，根據間接平差方法，令[a=m*cosα-1] ， [b=m*sinα]， 按最小二乘原理推導[3]，限于篇幅直接給出公式如下：

[L=x21-x11y21-y11…x2i-x1iy2i-y1i…x2n-x1ny2n-y1n，B=1001…10…10…01…01x11-y11y11x11…x1iy1i…x1ny1n…-y1ix1i…-y1nx1n，]

[X=?x?yab，P=1001000000001001]? ?（2）

因為一個點只能建立兩個誤差方程，要解四個未知數，至少需要兩個點建立四個誤差方程，如果有多的點，就使用最小二乘法。

求得結果X=（BTPB）-1（BTPL），可以得到Δx、Δy、a、b的值，進而可以解得α=acrtan[b/（a+1）]、m=sqrt[（a+b）2+b2]，得到Δx、Δy、α、m四個參數值[4-5]。

1.3 四參數轉換的可靠精度

轉換參數的求取時，所采用的已知點的位置分布需要注意如下幾個方面：已知點的選取位置應分布合理均勻，分布于測區四周和中心。為達到既定轉換精度，盡量采用多個已知點，讓這些點位能完全并均勻覆蓋整個轉換區域。留取幾個檢查點，作為檢核。為了保證參數求解精度，需要注意以下幾點：

1） 在參數求解前，對已知坐標進行預處理，保證數據準確。

2） 控制坐標分布、密度。

3） 殘差分析。參數求解后對殘差進行綜合分析，剔除誤差較大的點，提高可靠性。

2 Excel的電子表格編輯設計

在 Excel中，函數實際上是一個預先編寫的特定計算公式。按照這個特定的計算公式對一個或多個值執行計算，并得出一個或多個運算結果，叫作函數值。使用這些函數不僅可以完成許多復雜的計算，而且還可以簡化公式的繁雜程度。使用函數可以簡化或縮短工作表中的公式，使數據處理簡單方便。

Excel每一個函數都由一系列的參數來控制其輸入輸出結果。函數的參數可以是數值、文本、單元格引用等數據類型，將函數的參數傳入函數中進行處理得到結果。其中，一些函數會有循環嵌套，即通過一個函數調用另一個函數的方式進行計算，從而實現更為復雜的計算需求。

2.1 所用函數

實現坐標轉換需要用到矩陣運算，Excel中有矩陣運算相關函數可以方便得到所需結果，常見的矩陣函數有：1） 矩陣轉置：TRANSPOSE（原矩陣）；2） 矩陣乘法：A、B矩陣相乘MMULT（矩陣A，矩陣B）；3） 求逆矩陣：MINVERSE（原矩陣）；4） 求矩陣行列式的值 MDETERM（矩陣）；5） INDIRECT函數。在使用矩陣函數時候在公式編輯欄中輸入公式：如“=MINVERSE（矩陣A）”，按下“Shift+Ctrl+Enter”組合鍵，即可得到對應的逆矩陣。

2.2 具體編輯設計

為了防止數據編輯出錯、界面簡潔，在Excel中新建兩個工作表并命名為“輸入”和“過程”。在“輸入”工作表中編輯輸入已知點和坐標轉換如圖1：

界面數據錄入分為目標坐標和源坐標，各自對應X、Y坐標。按行從上到下排列分為P1、P2、P3…Pn點。右側F列編輯殘差列，對應每個點轉換后的殘差，來監測已知點的數據質量。殘差的計算方式為：1） 由已知點求得轉換參數；2） 源坐標按第一步的轉換參數求出轉換后的坐標；3） 轉換后的坐標和已知目標坐標X、Y進行求差，再求平方和、開方即為殘差。

在“過程”工作表中編輯生成過程矩陣，其中L為2n行1列的矩陣，B為2n行4列的矩陣（n為已知點的數量）。實際應用中已知點數量為2個及以上，數量并不是固定的，因此L、B兩個矩陣宜豎向排列，以便在點數增加時進行矩陣擴展如圖2所示。

以L矩陣第1、2行為例，第一行輸入“=IF（F1>=1，INDIRECT（"輸入！D3"）-INDIRECT（"輸入！B3"），0）”，第二行輸入“=IF（F1>=1，INDIRECT（"輸入！E3"）-INDIRECT（"輸入！C3"），0）”。F1單元格為已知點數量可以用COUNT函數取得，即已知點數量大于等于1時第一行值為x21-x11否則為0。因為Excel的公式會隨著復制剪切同步變化，這會導致已知點編輯時引起后面公式變化，所以加入INDIRECT函數防止剪切單元格時公式跟隨變化。以此類推根據公式2中L矩陣編輯其余單元格。

矩陣B第一行第一列輸入“=IF（F1>=1，1，0）”即已知點數量大于等于1，第一列值為1，否則為0;第二列固定輸入0，第三列輸入“=IF（F1>=1，INDIRECT（"輸入！B3"），0）”，即已知點大于等于1時，第三列單元格值為“輸入”B3單元格的值，否則為零；第四列輸入“=IF（F1>=1，INDIRECT（"輸入！C3"）*-1，0）”，即已知點F1數量大于1時，第四列單元格值為“輸入”C3單元格值乘以-1，否則為零。據此根據需要擴展矩陣的行數，此例為5個點，當輸入點數量少于5時，空白的地方利用IF函數根據F1單元格的點數自動填充為0。即已知點數量小于等于5個時，矩陣會自動把不足5的地方填充為零，不影響最終結果的計算。可以達到已知點數量動態使用的效果。如果需要更多的已知點數量，只需要手動把矩陣范圍擴展到相應點數即可。

編輯求得B矩陣的轉置矩陣BT，因Excel的TRANSPOSE函數使用時需要固定單元格數量，所以矩陣BT采用手動引用編輯生成，以便矩陣的跟隨矩陣B動態擴展。

矩陣BT第一行從左到右等于矩陣B的第一列從上到下，矩陣BT第二行從左到右等于矩陣B的第二列從上到下，依次編輯好矩陣BT，結果如圖3所示。

編輯構造矩陣BTPB，權系數矩陣P默認為1的單位陣，在此忽略，BTPB結果為4行4列的矩陣，它的大小是不變的，和已知點數量無關。在空白處選中4行4列單元格，在編輯欄輸入“=MMULT（矩陣BT，矩陣B）”，然后同按Ctrl+Shift+回車鍵返回得到結果圖4所示：

求逆矩陣（BTPB）-1，逆矩陣大小和源矩陣大小相同，在空白處選中4行4列單元格，在編輯欄輸入“=MINVERSE（矩陣BTPB）” 然后同按Ctrl+Shift+回車鍵返回得到結果圖5所示：

編輯構造矩陣BTPL， 權系數矩陣P默認為1的單位陣，在此忽略。矩陣BTPL結果為4行1列的矩陣，它的大小也是不變的，和已知點數量無關。空白處選中4行1列單元格，在編輯欄輸入“=MMULT（矩陣BT，矩陣L）”，然后同按Ctrl+Shift+回車鍵返回得到結果圖6左邊BTPL：

求方程的最小二乘解X=（BTPB）-1*BTPL，其結果為4行1列固定大小矩陣。空白處選中4行1列單元格，在編輯欄輸入“=MMULT（矩陣（BTPB）-1， 矩陣BTPL）”，引用好各個參數矩陣的位置，然后同時按Ctrl+Shift+Enter鍵返回得到矩陣X結果圖6右邊所示，分別為Δx、Δy、a、b的值。按照公式2最終解得α=acrtan[b/（a+1）]、m=sqrt[（a+b）2+b2]，得到Δx、Δy、α、m四個參數值。

3 實例與結果分析

現有某工程已知點坐標如表1。

如圖7所示，輸入P1、P2、P3三個點得到四個轉換參數平移Δx、Δx、α、m分別為53.193m、-50.134m、6.578E-6弧度、1.0000003。殘差均小于1mm。并有殘差列可供使用者評估各個已知點的精度。轉換后的參數可以用于實際測量中。

4 結束語

本文針對工程測量中重要的坐標轉換實現問題進行研究，探討二維平面坐標轉換的實現原理、精度控制和注意事項。經過矩陣運算、實現使用Excel 內部函數進行坐標轉換。經頁面設計和調試測試達到精度要求。該方法相對于專用軟件和Excel VBA編程方法，實現簡單，無需復雜的編程過程、學習成本小；過程數據可見，過程中生成的矩陣可以分析數據，更能方便理解其轉換的意義、錄入保存數據方便，適用于多種測量、教學場景。

參考文獻：

[1] 李征航，黃勁松.GPS測量與數據處理[M].3版.武漢：武漢大學出版社，2016.

[2] 馮注龍，丘榮茂.Excel之光：高效工作的Excel完全手冊[M].北京：電子工業出版社，2019.

[3] 武漢大學測繪學院測量平差學科組.誤差理論與測量平差基礎[M].2版.武漢：武漢大學出版社，2009.

[4] 王萌，高遠，范咪娜.區域坐標轉換模型對比分析[J].礦山測量，2021，49（2）：71-76.

[5] 謝偉杰.平面四參數法坐標轉換的實例應用探討[J].工程建設與設計，2020（13）：46-47，50.

【通聯編輯：梁書】