基于DEM實(shí)現(xiàn)城市高斯投影變形可視化分析

尚金光，張小波，陳軍勝，馮威

(1.成都市勘察測繪研究院，四川 成都 610081； 2.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 610031)

1 引 言

平面控制測量中，將地面長度投影到參考橢球面以及將參考橢球面長度投影到高斯平面均會引起長度變形[1]。關(guān)于投影變形分析的常規(guī)思路，多基于存檔控制點(diǎn)資料進(jìn)行，如李春華、李江衛(wèi)等利用城市C級GPS點(diǎn)分別對成都市和武漢市的控制網(wǎng)進(jìn)行長度投影分析[2，3]。后來，王忠禮博士認(rèn)為獨(dú)立坐標(biāo)系構(gòu)建中的平均高程取值含糊，提出采用地圖產(chǎn)品獲取測區(qū)平均高程[4]，對抵償高計(jì)算進(jìn)行了優(yōu)化，初步將面狀分析引入到投影變形中。隨著分辨率高、覆蓋均勻、地貌特征詳細(xì)的DEM在城市精準(zhǔn)規(guī)劃工作中的應(yīng)用[5]，更加直觀的面域投影分析將逐步更替?zhèn)鹘y(tǒng)的點(diǎn)狀分析。結(jié)合日益擴(kuò)張的各大城市建設(shè)幅面，提出基于DEM獲取整個城市滿覆蓋的面狀“連續(xù)”點(diǎn)位來實(shí)現(xiàn)高斯投影變形的精細(xì)分析評估。文章首先介紹關(guān)于DEM解析以及獨(dú)立坐標(biāo)系投影變形的分析，通過自編軟件實(shí)現(xiàn)自動計(jì)算與可視化。以成都市為例，結(jié)合行政界線，分析了投影變形分布情況、最佳投影參數(shù)驗(yàn)算以及簡陽納入成都的投影分帶等問題。

2 DEM解析及投影變形計(jì)算理論

2.1 DEM讀取與格網(wǎng)點(diǎn)解析

如圖1所示，DEM格網(wǎng)是由一連串的四邊形相互連接而成，每個四邊形由4個頂點(diǎn)組成[6]。

文本形式的格網(wǎng)數(shù)據(jù)(*.asc)由兩部分組成：前6行表示文件頭的基本信息，從第7行起為數(shù)據(jù)點(diǎn)高程值。如表1所示為成都市DEM文件頭信息。

圖1 DEM格網(wǎng)組織示意圖

DEM文件頭基本信息 表1

解析格網(wǎng)點(diǎn)的大地坐標(biāo)(Li，Bj)可通過其對應(yīng)的行列號i，j和左下角坐標(biāo)計(jì)算求得，公式為：

(1)

高程值與格網(wǎng)點(diǎn)一一對應(yīng)，從而形成滿覆蓋的“控制點(diǎn)”。

2.2 投影變形及有效控制區(qū)分析

在大地測量計(jì)算中，兩項(xiàng)投影引起長度變形規(guī)律為[7～11]：①實(shí)測邊長歸算到參考橢球面，會引起長度縮短；②參考橢球面上的邊長歸算到高斯投影面，會引起長度增長。

規(guī)范中的綜合投影長度變形是指以上兩項(xiàng)歸算變形的綜合影響，即要求不大于 2.5 cm/km。不在國家標(biāo)準(zhǔn)分帶中央子午線附近的城市，會受投影變形的影響，而選擇城市中心作為中央子午線，平均高程面作為抵償面，通過橢球膨脹等方式[12，13]建立城市獨(dú)立坐標(biāo)系。因而，以抵償面建立的獨(dú)立坐標(biāo)系統(tǒng)投影相對變形(每公里長度變形值)公式為：

(2)

式中，Hm為大地高與抵償面高程之差，ym為平均橫坐標(biāo)。RA為測線所在方向參考橢球面法截弧曲率半徑，Rm為測線兩端平均緯度處參考橢球面平均曲率半徑。實(shí)際計(jì)算中，一般認(rèn)為RA≈Rm[9]。

成都市域平均緯度30°34′，則Rm≈6367.777 km。代入式(2)，可推導(dǎo)出成都市平面坐標(biāo)系統(tǒng)閾值邊界，相對抵償面高H與中央子午線帶寬y的關(guān)系式：

(3)

以y為橫軸，H為縱軸，極限值關(guān)系曲線如圖2所示：

圖2 投影歸算有效區(qū)圖示

由圖2可見，投影歸算長度變形是存在正負(fù)變形邊界的，在滿足城市測量要求的情況下，理論最大抵償范圍約 128 km，即“投影變形有效控制區(qū)”。而成都市域東西寬 180 km，已然超出有效控制區(qū)理論范圍，因而需要進(jìn)行更加詳細(xì)的分析。

3 數(shù)據(jù)處理基本流程及軟件實(shí)現(xiàn)

如圖3所示，數(shù)據(jù)處理的核心就是將數(shù)字高程格網(wǎng)(DEM_GRID)換算成投影變形格網(wǎng)(DEF_GRID)。

圖3 DEM投影變形流程示意

數(shù)據(jù)處理及可視化的基本流程概括為以下7個步驟：

①準(zhǔn)備DEM和驗(yàn)算數(shù)據(jù)；

②利用ArcGIS軟件將DEM轉(zhuǎn)成格網(wǎng)文件；

③解析DEM GRID以及后續(xù)的各項(xiàng)計(jì)算；

④高斯投影計(jì)算每個格網(wǎng)點(diǎn)的橫坐標(biāo)；

⑤計(jì)算每個格網(wǎng)點(diǎn)的綜合投影長度變形量；

⑥定制繪圖用的格網(wǎng)文件、譜色文件等；

⑦自動調(diào)用Surfer.dll進(jìn)行各類可視化。

每個格網(wǎng)點(diǎn)的變形量計(jì)算結(jié)束后，可通過“返賦值”的方式將DEM格網(wǎng)的高程值替換為變形量。然后，引用Surfer.dll建立對象實(shí)例，進(jìn)行賦值和操作。如表2為本文調(diào)用函數(shù)：

Surfer可視化函數(shù) 表2

這其中，譜色文件(*.lvl)的制作最為重要。建議通過設(shè)定投影變形上下限 100 mm～0 mm、譜色間隔 5 mm，動態(tài)計(jì)算漸變填充色(百分制分子)。

為了實(shí)現(xiàn)以上功能，作者采用C#語言開發(fā)了“基于DEM的高斯投影變形計(jì)算分析軟件(DEM_GSPA V1.0)”，實(shí)現(xiàn)了DEM讀取、計(jì)算、統(tǒng)計(jì)、成圖一體化操作，界面如圖4所示：

圖4 DEM_GSPA軟件界面示意圖

4 成都市投影長度變形精細(xì)分析

4.1 成都市概況及數(shù)據(jù)準(zhǔn)備

如圖5所示，成都市西鄰龍門山山脈、東有龍泉山山脈，東西跨度約 180 km，海拔最低 334 m，最高 5 353 m。僅看跨度和高差，已然超出投影變形閾值的理論最大可控范圍。2005年3月，國家測繪局，批準(zhǔn)成都市平面坐標(biāo)系統(tǒng)是成都市唯一合法的地方獨(dú)立坐標(biāo)系統(tǒng)，劃分為三個投影帶(中帶、東帶、西帶，如圖5所示)。近年，隨著成都市建設(shè)范圍逐步外擴(kuò)以及簡陽納入成都管轄，投影分帶問題再次列為年度督辦任務(wù)。以往簡單的控制點(diǎn)位分布分析過于蒼白，有必要從全市域角度進(jìn)行一次直觀精細(xì)分析，以明確各行政區(qū)投影變形情況，為坐標(biāo)換帶/分區(qū)轉(zhuǎn)換等工作打下基礎(chǔ)。

圖5 成都市域數(shù)字高程模型

本文使用的數(shù)據(jù)均來自成都市地理信息中心，如表3所示。

數(shù)據(jù)列表 表3

考慮到計(jì)算效率問題，實(shí)際計(jì)算前將 10 m分辨率DEM抽稀成 50 m。作者從SRTM下載了90m分辨率的“srtm_57_06.asc”，精度滿足計(jì)算要求[3，14]，DEM獲取困難的區(qū)域可考慮。

4.2 投影變形有效控制區(qū)分析

采用DEM_GSPA軟件，導(dǎo)入DEM數(shù)據(jù)，設(shè)置已知的投影參數(shù)[2]，疊加行政區(qū)劃界線，一鍵生成中帶投影變形分布圖。如圖6所示，中帶投影變形分布與成都市的地貌特征較為吻合，呈現(xiàn)中間小、東西大的規(guī)律，即成都平原變形小、山區(qū)變形大。宏觀上，高程起伏成為成都市投影變形變化快慢的決定性因素。

將2.5 cm/km的限差等值線作為分界線，超限區(qū)域全部變成灰色，剩下的區(qū)域則為1.2節(jié)所提到的投影變形有效控制區(qū)。有效控制區(qū)的覆蓋范圍對投影參

圖6 中帶投影變形分布圖(單位/mm)

數(shù)的選取、優(yōu)化具有非常重要的指導(dǎo)意義。如圖7所示，中帶的有效控制區(qū)在103.5°～104.5°之間，跨度約 100 km。其中，中心城區(qū)、新都、龍泉驛、青白江、天府新區(qū)、雙流區(qū)、溫江、郫縣、新津、崇州、都江堰和彭州的非山區(qū)等在中帶控制范圍內(nèi)，而西邊的大邑、邛崍、蒲江，東部的金堂、簡陽均大范圍超限，這與參考文獻(xiàn)[2]通過控制點(diǎn)分析結(jié)果是一致的。

圖7 中帶投影變形有效控制區(qū)(單位/mm)

4.3 簡陽納入成都市的投影分析

2016年5月16日，經(jīng)國務(wù)院批準(zhǔn)，同意將資陽市代管的縣級簡陽市改由成都市代管(川府函〔2016〕90號)。簡陽市有天府國際機(jī)場、空港新城、龍簡新城、簡陽中心城等重點(diǎn)建設(shè)區(qū)域，時間緊迫，需要盡快將簡陽市統(tǒng)一至成都市平面坐標(biāo)系統(tǒng)。根據(jù)規(guī)劃建設(shè)需求，簡陽納入成都的投影參數(shù)選擇的原則為：

簡陽參數(shù)選擇原則 表4

由此只需要考慮中帶、東帶或者組合情況。圖7顯示簡陽超過一半的區(qū)域不在中帶有效控制區(qū)范圍內(nèi)，東部邊緣變形量達(dá)到 7.17 cm/km，因而需要轉(zhuǎn)到東帶分析。如圖8所示，在東帶投影參數(shù)下，91.2%的區(qū)域不超限，東部邊緣略超限(2.94 cm/km)，但此區(qū)域建設(shè)量極少。

圖8 東帶投影變形有效控制區(qū)(單位/mm)

因此，東帶初步成為比較好的選擇方案。中東帶組合方案的核心就是實(shí)地尋找與投影分界線吻合的地理參照。如圖9所示，疊加簡陽規(guī)劃要素，圖上主要地理參照(沱江、成渝客專、二繞、三繞等)均與投影帶分界線不相吻合，此方案不具有實(shí)際意義。因而，最終確定簡陽區(qū)域采用成都市平面坐標(biāo)系統(tǒng)東帶投影參數(shù)。

圖9 簡陽中東帶投影分界線選擇示意圖

5 數(shù)據(jù)驗(yàn)算分析

5.1 投影變形正確性驗(yàn)算

如表5所示均勻選取具有代表性的10個C級點(diǎn)與參考文獻(xiàn)[2]投影計(jì)算結(jié)果進(jìn)行比對，計(jì)算結(jié)果完全一致。

投影變形計(jì)算正確性檢驗(yàn) 表5

5.2 最佳投影參數(shù)驗(yàn)算

以DEM格網(wǎng)點(diǎn)來求取城市最佳投影參數(shù)。按照1′的經(jīng)差、0.1 m的高差為步長，進(jìn)行遍歷搜索[3]。原則依次為：①超限點(diǎn)數(shù)最少；②最大變形量最小；③變形量平方和最小。搜索判斷條件如下[3，15，16]：

(4)

為達(dá)到對比效果，參照2005年成都市控制網(wǎng)改擴(kuò)建項(xiàng)目中的點(diǎn)位分布，利用Surfer拾取相同點(diǎn)位、采用DEM_GSPA搜索最優(yōu)投影參數(shù)進(jìn)行比對。

投影參數(shù)驗(yàn)算(保密原因，*表示) 表6

經(jīng)過1 754次搜索計(jì)算，如表6所示，最佳投影參數(shù)與既有參數(shù)非常接近，中央子午線僅差1′，投影面大地高僅差 2 m，變形量考核指標(biāo)也非常接近。數(shù)據(jù)證實(shí)成都市測繪基準(zhǔn)既有投影參數(shù)(中帶)十分接近程序搜索最優(yōu)值，說明參數(shù)已經(jīng)最優(yōu)，同時驗(yàn)證了軟件計(jì)算的正確性。

6 結(jié) 語

基于DEM精細(xì)數(shù)據(jù)的投影變形分布分辨率高、覆蓋面廣，為投影變形的精細(xì)分析及城市坐標(biāo)系的建立或改造提供了翔實(shí)數(shù)據(jù)基礎(chǔ)。對變形格網(wǎng)進(jìn)行可視化后，可根據(jù)限差等值線建立有效控制區(qū)，結(jié)合行政區(qū)劃，疊加國土、規(guī)劃、交通等部門專題要素，進(jìn)行城市投影參數(shù)分析與調(diào)整。基于DEM進(jìn)行投影變形計(jì)算，相對于傳統(tǒng)計(jì)算方法，結(jié)果更精細(xì)直觀。其他城市均可參照此方法，建立城市精細(xì)結(jié)構(gòu)投影變形格網(wǎng)，從而檢驗(yàn)、改進(jìn)投影參數(shù)，更好地為城市規(guī)劃建設(shè)服務(wù)。