一種基于DEM的CGCS2000城市平面坐標系統確定方法

高昭良

(1.武漢大學測繪遙感信息工程國家重點實驗室，湖北武漢 430079;2.福州市勘測院，福建福州 350003)

1 引言

從2008年7月1日起正式啟用中國大地坐標系統2000(CGCS 2000)作為國家法定的坐標系，作為我國新一代的平面基準［1］。該系統有利于GPS快速、精確地獲取高精度城市坐標和高程成果，有利于GIS與GPS應用的結合，進一步提升數字城市的綜合服務能力。

目前全國各地城市基本還是沿用原有的1954北京坐標系、1980西安坐標系統及地方坐標系統，因此，各城市都需要進行采用CGCS2000建立城市平面坐標系統。陳俊勇院士等探討國家小比例尺地圖的“西安80坐標系”更換為“大地參考系1980”(GRS80)，給出這些變化量的大小區間并討論它對地圖表示可能產生的影響［2］。李江衛、解斌等人以武漢市為例，提出了以投影改正平方和最小和投影變形大于 2.5 cm/km的點數最少為原則進行最佳中央子午線和抵償高程面選取的方法［3］。武豐雷、劉曦燦等人提出建立斜軸墨卡托投影濟南市坐標系［4］。

在討論如何確定城市平面坐標系統選擇時，相關文獻基本采用兩個方面來驗證:一是采用有限的控制點作為計算結果驗證，二是采用最小二乘法。由于城市區域范圍較大，地形起伏分布不均勻，因此，有限的控制點無法全面準確概括。并且少量的地形突變會導致最小二乘法的最優評價結果畸變。本文提出的基于DEM的CGCS2000城市平面坐標系統確定方法可以兼顧解決這兩方面問題。

2 城市平面坐標系統定義模型

一個城市坐標系統確定應根據以上基本原則，還得考慮城市具體地理范圍和形狀，最后來確定城市平面坐標系統的定義模型。根據《城市測量規范》，坐標系統應遵循以下原則:

(1)投影變形盡量小，投影長度變形值不大于2.5 cm/km(1/4 萬)［5］。

(2)坐標系與國家(省)級控制點應聯測，并建立嚴密的轉換關系［6］。

2.1 高斯投影長度變形分析

由于定義國家大地坐標系的橢球面是一個凸起的、不可展平的曲面，當采用高斯正形投影將這個曲面上的元素投影到平面上時，投影后就會發生長度變形問題。高斯投影時，先把地面觀測值歸化至參考橢球面上，再把參考橢球面上的觀測值歸化至高斯平面上。兩次歸化產生的長度變形值稱為投影變形。投影變形△D計算過程如下:

上式中各變量說明如下:

Hm:測量兩端為高出參考橢球面的平均高程，單位:m;

H0:抵償面高程，單位:m;

D:測量長度兩端點平均高程面的水平距離，單位:m;

RA:歸算邊方向參考橢球法截弧的曲率半徑，可取概值 6 371 000 m，單位:m;

ym:測距邊兩端點投影后近似橫坐標平均值，單位:m;

Rm:測距邊中點的平均曲率半徑，可取概值6 371 000 m，單位:m。

2.2 采用格網占比法確定最佳關系

對于面積較大的城市或地區，在確定CGCS2000框架下的平面坐標系時，還應根據區域地理特點以及建設發展的實際情況，合理選取最佳中央子午線和抵償高程面。不斷調整高程參考面的大地高為H，使其投影變形滿足要求，條件式可近似表示為:

若確定的中央子午線與高程抵償面在城市坐標系統覆蓋的范圍內，能夠較好地滿足式(2)或式(3)的關系，則該中央子午線與高程抵償面關系最佳。

3 基于DEM模型數據的計算模型

大部分學者計算最優值采用最小二乘法∑△D2=min［9，10］模型，由于選取的控制點有限，不能全面均勻代表該區域投影變形，因此，條件通過遍歷搜索獲得∑△D2=min的方案是否可以確定為最優，值得商榷。本文基于福州地區 1∶1萬DEM數據，提取分布均勻的高程網格，利用以上的投影計算模型，計算最佳高程抵償面。最佳高程抵償面與中央子午線的計算采取遍歷搜索尋優法:搜索范圍在中央經度［Lmin，Lmax］，抵償面高程［Hmin，Hmax］內，按適合精度的步長，以投影改正平方和最小為原則，以變形值小于 2.5 cm的格網占比(符合變形要求的格網個數/坐標定義覆蓋區域的所有格網個數:m/n，if g(i)＜2.5)then m=m+1，n為所有格網數)作為參考觀察指標，按逐步縮小搜索區間，直到理論最優值所在范圍內達到預期設定要求，最后進行幾種最優方案比選。本文將這種方法簡稱為格網占比法。

3.1 基于福州地區DEM提取500 m格網間距高程點

為了避免采用有限的控制點不夠均勻，代表性不夠，從而導致最優關系計算不夠嚴密。因此，利用福州地區 1∶1萬DEM數據，福州地區 1∶1萬DEM數據格網尺寸為 12.5 m，高程精度 1 m。坐標系統為1980西安坐標系，1985黃海高程基準。根據福州山地地形的變化情況，選擇邊長 250 m格網來提取高程，基本可以滿足計算需要。按一定的格網間距提取三維坐標點，進行最優關系計算。

3.2 最優中央子午線和抵償高程面確定算法流程

(1)基礎數據準備:按投影經度步長2分(約3 km)、抵償面高程步長 50 m設置，生成所有計算方案表1、基于CGCS2000坐標系統的格網中心點數據表2。

(2)按遍歷計算每個方案每個格網點的變形值△D，存入表3。

(3)完成一個方案的變形值計算后，將∑△D2和△D≤2.5 cm的格網占比結果計算存入表3的對應方案符合。

(4)最后遍歷搜索法對比最優方案。

具體算法流程圖如圖1所示:

圖1 計算流程圖

3.3 最優中央子午線和抵償高程面計算實驗結果

以福州地區數據為例進行實驗計算，選擇精度118°40'～119°40'范圍內，按 2 分步長遞增，高程抵償面從 0 m～600 m，按 20 m步長遞增，對福州地區范圍內的格網數據進行遍歷計算。計算結果得到最優方案不同，如下:

(1)格網占比法:按△D≤2.5 cm的網格占比最大計算結果及對比數據。

變形值符合要求的格網占比計算結果及對比數據表 表1

變形值符合要求的格網占比模型計算結果表1列出最大的5個方案，最優的方案是中央子午線為119°34'，高程抵償面為 50 m，符合每千米變形△D≤2.5 cm的網格占福州區域內格網總數的55.62%。但是每千米變形值平方和為39.78并非最小。空間分布專題圖如圖2(b)所示。

(2)最小二乘法:∑△D2最小值計算結果及對比數據。

每公里變形值平方和計算結果及對比數據表表2

∑△D2最小值模型計算結果表2列出最小的5個方案，最優的方案是中央子午線為119°30'，高程抵償面為 150 m，符合每千米變形值平方和最小為38.00。但是符合每千米變形△D≤2.5 cm的網格占福州區域內格網總數為37.90%。空間分布專題圖如圖2(b)所示。

圖2 格網占比法與最小二乘法最優方案對比圖

從圖2很明顯看出左圖的符合每千米變形△D≤2.5 cm的網格覆蓋量比右圖合理，左圖福州市中心及沿海經濟開發熱點平原區域基本符合要求，而右圖福州市中心城區大東邊區域與沿海平潭島及羅源區域大部分都不符合變形要求。顯然，按格網占比法計算模型結果更為科學合理、更符合實際應用需要。按福州地方平面坐標系統選擇的投影中央子午線為119°18'21″方 案［11］，計 算 結 果 為 符 合 每 千 米 變形△D≤2.5 cm的網格占比為46.04%，每千米變形值平方和值為64.15，并不是最優的方案。福州地區西北面基本為丘陵高山，許多高程基本在 400 m以上，經過遍歷搜索計算，驗證了該地形變化規律不適合采用高程抵償面長度歸化改正與投影變形進行抵消的計算模型。對于該高山區域的測圖，只能采用抬高投影面的方法。

4 結論

本文采用DEM數據與地理信息系統計算方法，對CGCS2000城市平面坐標系統定義的兩個計算模型進行研究，得出采用格網占比法方法得出的最優高程抵償面與中央子午線的選擇方案更具嚴密性，該方法對城市區域CGCS2000平面坐標系統定義具有較大的參考意義。

［1］陳俊勇，張鵬.國家測繪基準“十二五”重大項目的思考［J］.武漢大學學報·信息科學版，2009，34(10):1136～1138.

［2］陳俊勇.中國采用地心三維坐標系統對現有地圖的影響［J］.測繪學報，2004，32(4):283 ～288.

［3］李江衛，白潔，郭際明等.基于CGCS2000的城市平面坐標系最佳選取［J］.城市勘測，2011，3(4):120～123.

［4］武豐雷，劉曦燦，王方等.濟南市獨立坐標系的改造探討［J］.城市勘測，2010，2(2):105 ～106.

［5］CJJ T8－2011.城市測量規范［S］.

［6］吳云孫，晁定波，楊堂堂.淺談廣州市平面控制測量坐標系統的選取［J］.測繪通報，2005，23(8):42～43.

［7］陸鵬程，林冬偉.斜軸墨卡托投影模型及其應用分析［J］.鐵道勘察，2010，12(4):26～29.

［8］呂忠剛，許世寧.關于抵償高程面與移動中央子午線最佳選取問題的研究［J］.東北測繪，2003，25(2):3～8.

［9］謝征海，張澤烈.城市獨立坐標系平面控制網之擴建及改造［J］.北京測繪，2000，23(4):22～23.

［10］孔建，姚宜斌，吳寒.整體最小二乘的迭代解法［J］.武漢大學學報·信息科學版，2010，35(6):711～714.

［11］高昭良.基于相似變換的城市坐標系統轉換研究［J］.福州大學學報·自然科學版，2010，38(1):54～57.