一種基于DEM的城市平面坐標系統(tǒng)建立方法

王友昆，謝正明，張君華，余章蓉

(1.武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079；2．昆明市測繪研究院，云南 昆明 650051； 3．云南國土資源職業(yè)學(xué)院，云南 昆明 652501；4．昆明理工大學(xué) 津橋?qū)W院，云南 昆明 650106)

城市平面坐標系統(tǒng)建立，在高斯統(tǒng)一3°帶平面直角坐標系或高斯任意帶平面直角坐標系無法滿足長度變形規(guī)定時(≤25 mm/km)，可采用城市平均高程面作為投影面(抵償高程面)。傳統(tǒng)的投影面選擇，一般基于控制點的平均高程或者控制點控制范圍的平均高程確定[1-4]。該方法計算量小，速度快，但無法準確計算和分析符合長度變形規(guī)定的區(qū)域范圍。隨著計算機性能的提升，高分辨率數(shù)字高程模型(Digital Elevation Model,DEM)已經(jīng)用于確定城市平均高程面[5]，長度變形計算[6]，以及投影面的可視化選擇[7]。這些方法大多應(yīng)用于平原等低海拔區(qū)域，主要涉及中央經(jīng)線(投影帶)的選取，沒有涉及多個投影面的選擇，也沒有充分利用GIS空間分析和統(tǒng)計分析的方法對不同投影方案進行量化指標分析和統(tǒng)計。本文針對云南高海拔區(qū)域，城市地方坐標系的建立需考慮任意帶和多個投影面的特點[8]，利用DEM成果進行長度變形計算、空間分析和可視化輸出，并以昆明市城市平面坐標系建立為例進行詳細說明。

1 長度變形原理

地面上實測的邊長，內(nèi)業(yè)計算時應(yīng)經(jīng)過高程歸化至橢球面和高斯投影至高斯平面兩個步驟，兩者對于實測邊長的影響分別為縮短和伸長[9]?！冻鞘袦y量規(guī)范》規(guī)定，城市平面坐標系建立時，坐標反算的邊長與實測邊長應(yīng)盡可能相符，即邊長歸算至參考橢球面上的高程歸化和高斯正形投影的距離改化的總和(即長度變形)應(yīng)小于25 mm/km，如式(1)所示。

(1)

式中:H為歸化高程；R為地球平均曲率半徑；ym為投影邊長端點橫坐標平均值。

按照式(1)計算分析，當城市地區(qū)高程大于160 m或其平面位置距離3°帶中央子午線的東西方向距離(橫坐標)大于45 km，其長度變形均超過規(guī)定的1/40 000[10]。

若式(1)無法滿足長度變形要求，可采用任意帶平面直角坐標系統(tǒng)。仍然無法滿足時，可采用城市平均高程面作為投影面(抵償高程面)，人為地改變歸化高程使它與高斯投影的長度改化相抵償，如式(2)所示。

(2)

式中:Hd為抵償面高程值。

為提高長度變形計算的精確性，可采用投影點所在卯酉圈半徑N替代式(2)中的地球平均曲率半徑R，N的算式為：

(3)

式中:a為橢球長半軸；e為橢球第一偏心率；B為投影點緯度。

2 DEM計算方法

2.1 計算流程

長度變形原理揭示了城市平面坐標建立時需要滿足的條件和計算方法。利用DEM成果進行城市平面坐標系建立，主要流程包括數(shù)據(jù)準備、預(yù)處理、變形計算、分析統(tǒng)計和成果輸出等步驟(見圖1)。

圖1 計算流程圖

2.2 數(shù)據(jù)準備

需要收集各類范圍線、DEM成果數(shù)據(jù)和似大地水準面精化模型。范圍線包括城市各級行政區(qū)劃面、建成區(qū)和重點規(guī)劃區(qū)范圍，用于對DEM原始數(shù)據(jù)的裁剪和分析長度變形對這些區(qū)域的影響。DEM可利用1∶5萬、1∶1萬等高精度的DEM數(shù)據(jù)，若無法獲取時可利用開源數(shù)據(jù)替代(分辨率優(yōu)于100 m，精度優(yōu)于20 m)，如航天飛機雷達地形測繪任務(wù)(Shuttle Radar Topography Mission,SRTM)的DEM成果。似大地水準面精化模型主要用于將DEM的正常高成果改化為大地高成果。

2.3 預(yù)處理

1)統(tǒng)一坐標系統(tǒng)。將以上數(shù)據(jù)統(tǒng)一為2000國家大地坐標系或WGS84(兩者橢球定義的參數(shù)差異極小，對長度變形影響不大)。

2)DEM裁剪。DEM按照格網(wǎng)存儲，為減少計算量應(yīng)利用市級行政區(qū)劃對其進行裁剪，僅保留行政區(qū)劃范圍內(nèi)的DEM成果，區(qū)域外格網(wǎng)屬性值設(shè)置為NoData。

3)大地高計算。DEM成果為正常高，而長度變形計算式(1)和式(2)中H為大地高。我國西部區(qū)域高程異常達到-40 m，對長度變形計算影響較大。因此，需要利用區(qū)域似大地水準面模型，將DEM正常高改算為大地高，如式(4)所示。

H=h+ε.

(4)

式中:H為大地高；h為正常高；ε為高程異常。

如果沒有可利用的區(qū)域似大地水準面模型，也可利用公開的全球超高階地球重力場模型EGM2008，其中國區(qū)域精度優(yōu)于1 m[11-13]，能夠滿足投影計算的要求。

4)高程統(tǒng)計。對區(qū)域范圍的DEM高程進行統(tǒng)計，便于計算投影面的大致高程。

2.4 變形計算

變形計算遵循圖 1中方案1至方案6的順序，依次嘗試計算，直至找到符合長度形變要求的城市平面坐標系統(tǒng)。如首先采用方案1進行變形計算，即采用統(tǒng)一3°帶的平面直角坐標系進行DEM變形計算，如果方案1不符合變形要求，則嘗試方案2，以此類推。是否符合變形要求，根據(jù)分析統(tǒng)計結(jié)果來判定。

1)高斯投影。按照長度變形計算的原理，需要對DEM成果進行高斯投影，中央經(jīng)線L0按照圖1的變形計算順序設(shè)置，計算DEM格網(wǎng)的橫坐標y(應(yīng)減去加常數(shù)500 km)。

2)長度變形計算。DEM為規(guī)則格網(wǎng)，格網(wǎng)自帶坐標及高程，為快速計算每一格網(wǎng)的長度變形提供了便利。橫坐標y和卯酉圈半徑N可用一維數(shù)組存儲，大地高可用二維數(shù)組存儲。利用成熟的矩陣庫(如Numpy)，可大大提高長度變形的計算效率。

3)格網(wǎng)標識。對長度變形值小于1/40 000的格網(wǎng)，標識為0，否則標識為1，并將標識后的DEM格網(wǎng)輸出為柵格，用于后續(xù)的分析和統(tǒng)計。

2.5 分析統(tǒng)計

分析統(tǒng)計主要是對DEM變形標識的格網(wǎng)，利用GIS原理對其分析和統(tǒng)計，作為城市平面坐標系統(tǒng)建立投影帶和投影面確定的依據(jù)。

1)格網(wǎng)統(tǒng)計。通過統(tǒng)計符合長度變形要求的格網(wǎng)比例來衡量區(qū)域長度變形是否符合要求。假設(shè)占比達到預(yù)定的比例(如90%以上，即大部分區(qū)域長度變形滿足1/40 000的要求)，則認為該方案符合要求。

2)柵矢轉(zhuǎn)換。為直觀地展示符合長度變形要求的范圍，以及多個投影面或關(guān)注區(qū)域(城市重點經(jīng)濟建設(shè)區(qū)域)的長度變形情況，需要將變形標識為0的柵格轉(zhuǎn)換為矢量多邊形。

3)疊置分析。如果城市平面坐標系涉及多個投影面，應(yīng)利用縣(區(qū))行政區(qū)域范圍或重點經(jīng)濟建設(shè)區(qū)域范圍，同符合長度變形要求的矢量面進行疊置分析，以此作為選擇投影面的依據(jù)。對于某投影帶，n個投影面計算的符合長度變形的多邊形記為prji(i=1,2,…,n)，其疊置分析主要步驟如下：

a)將prji中占比面積最大的作為主投影面，記為prjmax；

b)遍歷prji(i=1,2,…,n,i≠max)，分別同prjmax疊置分析，采用prji刪除prjmax，即第i個投影面排除主投影面后符合變形要求的范圍，記為Erasei；

c)將Erasei面積最大的作為次投影面prjmax-1；

d)重復(fù)步驟b)和c)，分別確定出后續(xù)的投影面prjmax-2，prjmax-3…prj1。

以上疊置分析過程中，還應(yīng)同時考慮投影面的數(shù)量、各縣(市、區(qū))建成區(qū)以及重點經(jīng)濟建設(shè)區(qū)范圍，以及單個符合長度變形區(qū)域范圍面積是否最大化等因素。

2.6 成果輸出

將DEM計算過程的所有成果，包括裁剪和投影后的DEM成果，各類矢量多邊形范圍線，不同計算方案的長度變形計算結(jié)果和統(tǒng)計分析結(jié)果，以數(shù)據(jù)庫和圖、表的方式輸出，用于成果報告的編制。

3 實例應(yīng)用

3.1 需求分析

由于歷史原因，昆明市1987昆明坐標系、2004昆明坐標系以及2013昆明坐標系等多個坐標系統(tǒng)并存，給城市建設(shè)和成果共享帶來諸多不便。按照國家及昆明市政府相關(guān)要求，昆明市啟動了建立基于CGCS2000橢球的昆明2000坐標系。昆明市轄7個區(qū)、3個縣、代管1個縣級市和3個自治縣，總面積約2.1萬km2，東西跨度150 km，南北跨度240 km，海拔范圍746～4 247 m，地勢地貌呈現(xiàn)北高南低、分布不均勻，且各縣(市、區(qū))建成區(qū)高程各不相同(見圖2)，黑色實線為昆明市主城區(qū)及各縣(市、區(qū))建成區(qū)范圍。

圖2 昆明市DEM

3.2 計算過程

1)數(shù)據(jù)準備。DEM采用STRM1數(shù)據(jù)，分辨率為30 m，高程精度為1 m；范圍線收集了昆明市及各縣(市、區(qū))行政區(qū)域范圍，以及各縣(市、區(qū))建成區(qū)范圍線；區(qū)域似大地水準面精化模型利用EGM2008成果；以上成果坐標系統(tǒng)統(tǒng)一為CGCS2000。

2)數(shù)據(jù)預(yù)處理。DEM經(jīng)過裁剪、大地高計算、高程統(tǒng)計[14]等預(yù)處理。如表1所示，昆明市高程值跨度較大，1 500～2 500 m占86.4%，高程大于2 500 m的比例達到9.7%。

表1 昆明市DEM高程值區(qū)間統(tǒng)計

3)變形計算和分析統(tǒng)計。按照圖1變形計算的順序，經(jīng)試算和分析，采用統(tǒng)一3°帶等計算方案(方案1—方案4)均無法滿足變形要求，可采用多個投影面的任意帶(單中央經(jīng)線)計算方案。為便于后期成果的接邊和管理，方案優(yōu)先選擇順序：①滿足昆明市主城區(qū)長度變形要求；②滿足各縣(市、區(qū))建成區(qū)長度變形要求；③符合要求的長度變形區(qū)域范圍能夠連片且面積最大化；④投影面數(shù)量盡量少；⑤多個投影面符合變形要求的區(qū)域面積之和達到一定的比例(如占昆明市域面積90%)。按照以上原則，利用程序?qū)Σ煌队皫?102°42′～103°09′為區(qū)間，經(jīng)差1′為間隔)、不同投影面(700～4 300 m為區(qū)間，10 m為間隔)的長度變形結(jié)果分析，選擇最優(yōu)的投影帶和投影面(見圖3)。由于涉及大量空間分析計算，為減少計算量，初算可采用經(jīng)差6′或3′，投影面間隔100 m或50 m，確定大致區(qū)間范圍后，再縮小計算間隔，減少冗余計算，提高計算效率。

圖3 投影面面積對比圖

圖3顯示不同投影帶滿足原則：①符合變形區(qū)域范圍的最大面積(虛線星號)和單個連片最大面積(實線圓點)對比圖。若僅考慮覆蓋面積最大，隨著中央經(jīng)線增大，覆蓋范圍逐漸增大，L0=102°56′時達到最大值8 615.88 km2(單個連片最大面積僅為5 598.20 km2),之后逐漸減小。若綜合考慮原則③(單個連片最大面積)，L0=102°47′時達到最大值8 025.46 km2。此時，投影面為1 890 m，作為第1投影面(記為Prj1)更為合理。之后，將其他投影面多邊形同Prj1進行疊置分析(刪減策略)，按照原則②③④進行逐一分析和比選，確定第2投影面。以此類推，直至多個投影面覆蓋面積之和滿足原則⑤。

4)成果輸出。利用ArcGIS制圖模板[15]進行成果輸出，只需修改各圖層數(shù)據(jù)源即可達到統(tǒng)一的輸出標準，便于對不同計算方案結(jié)果比選。圖 4顯示了滿足原則①②③④不同投影帶(102°42′～102°57′，3′為間隔)的計算結(jié)果。紅色實線為昆明市域范圍，黑色實線為各縣(市、區(qū))建成區(qū)范圍，黑色虛線為投影中央經(jīng)線，綠色、藍色、粉色填充區(qū)域分別為第1、第2、第3投影面覆蓋范圍。圖4(d)—圖4(f)可見，第1投影面覆蓋面積逐漸增大，但單個連片面積逐漸減小，結(jié)果同圖 3是一致的。

如果考慮原則⑤，需設(shè)置第4和第5投影面，如表 2所示，此時符合長度變形要求的總面積達到19 256.36 km2，占市域面積的91.7%。未覆蓋區(qū)域海拔大多在3 000 m以上，且位于昆明市生態(tài)紅線保護范圍內(nèi)，經(jīng)濟建設(shè)較少，可在實際工程建設(shè)時再行考慮。

圖4 投影面范圍分布圖

表2 投影面面積指標統(tǒng)計表

4 結(jié)束語

本文利用DEM數(shù)據(jù)易于讀取、計算和分析的特點，應(yīng)用于高原城市平面坐標系統(tǒng)建立時投影帶和投影面的選擇，并充分利用GIS空間分析的優(yōu)勢，對長度變形區(qū)域和行政區(qū)域、重點建設(shè)區(qū)域進行空間疊置分析，定量地分析和統(tǒng)計不同投影方案的計算結(jié)果，并實現(xiàn)了成果的可視化輸出。同時，實現(xiàn)了軟件程序化批量計算，大大提高了計算的效率。這種方法使計算者能夠直觀、便利、及時地調(diào)整計算方案，較傳統(tǒng)方法更加有效和可靠，能更優(yōu)地選擇出投影帶和投影面。同理，該方法還可應(yīng)用于大型區(qū)域工程和線狀工程的平面坐標系建立。