最優化理論在橢球參數逆向求定中的應用

劉敏，王磊

(廣州市城市規劃勘測設計研究院，廣東 廣州 510060)

1 引 言

2008年我國正式啟用的2000國家大地坐標系(CGCS2000)是經國務院批準使用的新一代國家大地坐標系。具有三維、高精度、動態等特點，能更好地滿足各領域業務工作需要，更好地為經濟建設、社會公眾服務。統一采用2000坐標系，有利于各系統、部門之間資源和成果的共建共享，有利于促進各級基礎地理信息數據成果的有效銜接，有利于推動地方“多規合一”、不動產統一登記等重大工作順利開展。建立基于CGCS2000橢球體的地方坐標系，保持空間關系連續性，是大多數城市的優選策略。

許多城市獨立坐標系統建立較早，橢球定義并不明確，或者是由于城市的擴展，經過多次擴展后的控制網投影自由度并不滿足1/40000要求，也存在原有投影參數仍處于保密狀態，為解決以上3類問題，本文提出一種基于最優化理論來逆向求定投影自由度，通過固定CGCS2000橢球體的扁率建立基于CGCS2000橢球體的地方坐標系。

圖1 大地坐標與平面坐標轉換

如圖1所示，5個投影參數分別為，橢球參數(長軸a，扁率f，)、中央子午線l0、橫坐標偏移量△Y、縱坐標偏移量△X。

一般情況下，是已知5個投影參數來求定一組大地坐標對應的平面坐標(反之亦然)。比如國家坐標系統，前兩個用固定的橢球體區分，中央子午線采用標準經度(如3度帶采用3*帶號)。以上5個投影參數只要有1個與國家坐標系統不同，就可以稱之為獨立坐標系統。

本文的技術路線是：在已知一組同名點兩種坐標情況下，逆向求定5個投影參數。

從投影平面坐標系互換到地理坐標系之間的五個投影參數的組合情況，理論上說有無數個合理實現兩種坐標的轉換(在規定的合理精度內)。

在眾多的方案或組合中什么樣的方案是最優的，最接近事實的，就需要引入最優化理論。

2 逆向求定投影參數的技術路線

地理坐標系統GCS中用經緯度來確定橢球面上的點位。投影坐標系始終基于地理坐標系，即：“投影坐標系=地理坐標系+投影算法函數”[1]。通用的GIS一般提供投影參數設置，GCS在建立地理信息平臺或者數據信息共享上，比較有優勢，但是在地方城市建設管理上，必須考慮與平面坐標系統銜接問題。

為從形式上讓平面坐標與某個橢球體建立聯系，從正形角度出發應該優先考慮固定橢球扁率，其次固定長軸，如圖2所示。

3 最優化理論及MATLAB實現

最優化理論解決問題包含兩個步驟：

(1)建立數學模型：用數學語言來描述最優化問題。模型中的數學關系式反映了最優化問題所要達到的目標和各種約束條件。

圖2 投影自由度逆向求定的技術路線

(2)數學求解：數學模型建好以后，選擇合理的最優化方法進行求解。

(1)首先對圖像進行掃描，當掃到當前像素為1時選為種子，并給它賦予一個標簽，然后按照四鄰域準則，把與之相鄰的所有值為1的像素點都壓到棧內保存。

最優化理論常用的算法有兩類，一類是搜索算法，另一類是迭代算法。

迭代算法是從參數的某一初始猜測值出發，然后產生一系列的參數點，如果這個參數序列收斂到使目標函數極小的參數點，那么對充分大的參數就可用來求解優化解。其中最小二乘法是最常用的方法，就是對由若干個函數的平方和構成的目標函數求極小值的問題。本文采用的就是通過迭代尋找min{sum[ Fun(x)2+Fun(y)2}的最小值。在同名點對中找出在符合最小二乘條件下的一組投影參數，并實現編程化處理。

通過編程工程和函數比選，本文采用MATLAB優化工具箱提供的Lsqnonlin函數來解決非線性最小二乘法問題。其算法就是采用Trust-region-reflective(信賴域反射算法)和Levenberg-marquardt(加權信賴域算法)。

信賴域法是從給定的初始解出發，通過逐步迭代，不斷改進，直到獲得最優解為止。具有收斂性較好、可靠性高，有效性強等特點。Lsqnonlin函數處理的是非線性最小均方差問題。函數從初值開始，最后到找滿足函數Fun(x)均方誤差和最小的x值返回。選定合適的初值可以保障非線性最小二乘法收斂速度加快。

4 高精度Gauss-Krueger公式

4.1 傳統的Gauss-Krueger公式

Gauss-Krueger投影是一種等角橫軸切橢圓柱投影，在投影面上，中央子午線和赤道的投影都是直線，其他子午線和平行線是復雜曲線。傳統投影公式是利用黎曼柯西方程給出的實數形經差的冪級數函數展開。該公式適用范圍比較窄(經差在4°～6°)在反算過程中的不準確性表現明顯，經MATLAB程序測試，該公式正反算精度不高，在非線性最小二乘過程中收斂值偏大，不適用于本文要求。但這個公式在國內外教科書被廣泛采用。

4.2 Karney-Krueger公式

Karney-Krueger按精度也分兩類公式：

第一類稱為擴展的Krueger(1912年)公式。該公式是利用計算機代數系統Maxima擴展了Krueger(1912年)的原始系列(第三扁率n的4階)至8階甚至30階。在第三扁率n高階至6時，距離中央子午線 4 200 km，正反算的精度達到優于 5 nm[6]。該公式精度高，投影范圍廣(經差可以擴展至75°)，同時適用于高緯度地區。

擴展的Krueger公式，原理是利用等量緯度與歸化緯度關系，把等量緯度和經差作為變量拓展至復數域。再借助雙曲正弦函數及其與三角函數的關系，把投影的復變函數表達式拆為實虛分離的實數形式。

推導過程大致為可以用圖3表達：

圖3 Karney-Krueger第一類公式原理

第二類是采用雅克比橢球函數進行橢球積分，利用計算機代數系統Maxima可以達到任意精度。

該公式精度可以達到任意精度(在計算機內存允許的范圍內)，投影范圍更廣(適應全橢球體)。

本文采用的是擴展至6階的Karney-Krueger第一類公式進行高精度正反算。采用高精度的高斯投影正反算公式是本文逆向求定投影參數的關鍵所在。

5 實例驗證

5.1 迭代初值的選擇

根據公開文獻報道，廣州市坐標系統采用克拉索夫斯基橢球(1954年北京坐標系)為基礎，并采用不同于國家標準經線(114°)作為中央子午線，并對橫坐標偏移量△Y、縱坐標偏移量△Y同時進行了偏移，根據本文對投影自由度的定義，該市坐標系是典型的非國家坐標系。

橢球體長軸a，扁率f可以選用1975IAG橢球或CGCS2000橢球作為初始值帶入。中央子午線以城市中心經度113.50°為初值；橫坐標偏移量△Y、縱坐標偏移量△X選用城市中心同名點橫縱差值代入，如圖4所示。

圖4 起算點與檢測點分布示意圖

5.2 投影自由度迭代技術

根據本文技術路線，圖3中8個起算點和6個檢測點同時具有CGCS2000坐標和廣州市坐標系統[1]。

第一步，根據技術路線，采用MATLAB編寫程序，分三個函數執行。具體代碼可聯系本文作者或參見相關文獻[4]。其中納米級高斯正反算公式采用Karney-Krueger 第三橢球扁率階數6計算。得出以下5個投影自由度：

橢球體長軸a，橢球體扁率f、中央子午線l0、橫坐標偏移量△Y、縱坐標偏移量△Y；(因多方面原因，以下部分數據有所省略)

6378126.*** ；1/298.255456291；113.17**；411*3.***；-25296*5.***

第二步，為形式上建立獨立坐標與CGCS2000之間的關系，可以固定扁率為CGCS2000扁率，其他4個自由度分別為：

6378126.5**；1/298.257222101；113.17**；411*3.*** ；-25296*5.***；可以從第一步和第二步得出的橢球體長軸a變化、坐標偏移量△Y、縱坐標偏移量△Y均在分米級上，說明整個非線性最小二乘的優化解是收斂的。

第三步，為固化與CGCS2000橢球體的長軸關系可以約定在原有橢球長軸6378137修改至整數長度，再進行一次非線性最小二乘。

第四步，將以上各步驟得出的5個投影自由度代入納米級高斯反算公式，與已知點進行比較，坐標偏差均小于 2 mm，證明在全市域以上幾套投影自由均是滿足要求的，如要在形式上與CGCS2000建立聯系，可以分別采用第二步以后的各組投影自由度。

6 結 論

隨著計算機技術的發展，在城鄉規劃管理中直接(后臺數據庫或間接使用)采用CGCS2000大地坐標(BLH)是歷史發展的必然趨勢。

許多獨立坐標無確切的數學模型，一般四參數或七參數與CGCS2000建立聯系，本文通過逆向求定投影參數，建立地方坐標與CGCS2000簡單的高斯投影關系，有著積極的現實意義，也為大地坐標作為地理信息系統后臺數據庫奠定了數學基礎。

[1] 劉敏. 城市CORS站轉換到CGCS2000的技術研究[J]. 城市勘測，2013(2)：113～115.

[2] Karney，C.F.F. (2011)，'Transverse Mercator projection with an accuracy of a few nanometres'，Journal of Geodesy.

[3] 施一民. 現代大地控制測量[M]. 北京：測繪出版社，2003，220～225.

[4] 劉萬華，葉水全. 構建獨立坐標系與CGCS2000 坐標系轉換關系的研究[J]. 測繪與空間地理信息，2016，39(1) ：216～218.

[5] 邊少鋒，劉強，李忠美. 不分帶的高斯投影實數公式[J]. 測繪通報，2016，1(1) ：6～9.

[6] Lee，L.P. Conformal projections based on Jacobian elliptic functions，Cartographica[J]. 1976，13，128～129.

[7] Yang，Q.，Snyder，J.P. and Tobler，W.，Map Projection Transformation [M]. 2000，100-108，Taylor & Francis，London.

[8] Deakin，R.E. and Hunter，M.N. Geometric Geodesy Part B，School of Mathematical and Geospatial Sciences[J]. 2010，28～30.

[9] 劉強，邊少鋒，李忠美. 球面高斯投影及其變形的閉合公式[J]. 海軍工程大學學報，2015，27(1)：46～48.