劉敏，王磊

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

1 引 言

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

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

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

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

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

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

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

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

2 逆向求定投影參數(shù)的技術路線

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

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

3 最優(yōu)化理論及MATLAB實現(xiàn)

最優(yōu)化理論解決問題包含兩個步驟：

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

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

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

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

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

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

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

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

4 高精度Gauss-Krueger公式

4.1 傳統(tǒng)的Gauss-Krueger公式

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

4.2 Karney-Krueger公式

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

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

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

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

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

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

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

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

5 實例驗證

5.1 迭代初值的選擇

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

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

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

5.2 投影自由度迭代技術

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

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

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

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均在分米級上，說明整個非線性最小二乘的優(yōu)化解是收斂的。

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

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

6 結 論

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

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

[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] 施一民. 現(xiàn)代大地控制測量[M]. 北京：測繪出版社，2003，220～225.

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

[5] 邊少鋒，劉強，李忠美. 不分帶的高斯投影實數(shù)公式[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.