基于網絡地理信息服務的選線數字地形獲取方法

聶良濤， 易思蓉， 李 陽

(1. 西南交通大學土木工程學院，四川 成都610031;2. 西南交通大學道路工程四川省重點實驗室，四川 成都610031;3. 西南交通大學高速鐵路線路工程教育部重點實驗室，四川 成都610031)

隨著BIM 技術的大力推廣，在線路勘察設計領域，基于虛擬地理環境的三維數字化選線設計方法正在逐漸代替傳統的基于等高線地形圖的二維透視設計方法［1］. 開展三維數字化選線需要建立一個逼真的虛擬地理環境，該環境的建立依賴于高精度高程數據和影像數據的支撐.在道路前期規劃階段，一般尚未開展航測工作，往往面臨數字地形資料缺乏的問題，尤其在我國一些偏遠地區項目和海外項目建設中，由于基礎資料匱乏、資料收集困難，且前期設計周期短、資金有限，不具備足夠條件和時間進行航測或外業調繪，因此，在這種情況下，數字地形成了道路前期規劃階段開展三維數字化選線的制約因素.

隨著數字地球技術的發展，網絡上部分中高分辨率衛星遙感影像和高程數據免費開放給公眾［2-8］，與商業影像和高程數據相比，通過網絡免費獲取的影像和高程數據的分辨率及精度相對低一些，但仍能夠滿足公路前期規劃階段的要求［9］. 在線路勘測設計領域，針對網絡地理信息的獲取及利用已有許多有價值的研究成果［10-13］，并應用于現場設計工作中.早期的研究主要側重于通過網絡數字地球平臺提供的API(application programming interface)接口進行二次開發以獲取數字高程及影像.然而，這種獲取方式容易受API 函數的限制，例如，Google API 取點函數存在5 000 次限制，并且通過截圖函數獲取的影像，質量也不容易控制. 也有學者針對網絡地圖服務提供的影像瓦片進行了下載研究［14-15］，但是這種針對網絡可視化的影像瓦片數據不能直接應用于道路規劃選線，還需要進行匹配與校正. 并且隨著時間的推移和技術的進步，許多網絡地圖的數學模型也有更新，因此，有必要在綜合分析當前網絡地理信息資源的基礎上，對網絡地理信息的獲取及應用于道路數字化選線系統的關鍵技術進行系統分析與探討.

1 網絡地理信息資源分析

目前全球范圍精度較高且能直接通過網頁免費下載的高程數據的典型代表有SRTM3 DEM 和ASTER GDEM (advanced spaceborne thermal emission and reflection radiometer global digital elevation model)［16］. 兩種高程數據的參數對比分析見表1.

表1 SRTM3 DEM 和ASTER GDEM 參數對比Tab.1 Comparison of SRTM3 DEM andASTER GDEM parameters

從表1 的參數對比可知，SRTM3 DEM 精度優于ASTER GDEM，但ASTER GDEM 的分辨率更高.本文參考文獻［17］的研究建議，實際中選用高程數據時，在低海拔平原地區選用精度相對較好的SRTM3 數據;在地勢變化劇烈、地形復雜地區選用分辨率高的ASTER GDEM 數據. 由于兩種數據的基準及投影方式一致，本文僅選用SRTM3 數據進行測試.

對于網絡影像資源，能夠免費獲取的種類較多，主要分為兩類:

(1)初級影像產品，例如Landsat 系列影像、EO-1 衛星影像等，其影像多光譜波段分辨率為30 m.發布這類影像的組織會提供完整的影像信息(影像的來源精度、坐標、拍攝時間以及處理技術等)，但該類影像只經過了初級矯正，應用于實際生產前還需進行專業的處理;

(2)網絡地圖服務提供的影像，例如Google Maps、World Wind、Yahoo 地圖、百度地圖等提供的衛星影像.相對于前一類，這類影像的分辨率較高，最高分辨率為分米級，并且現勢性好. 網絡地圖服務除了提供衛星影像資料之外，還能提供較全的地名、河流、既有道路及邊界等標注圖層，非常有利于選線工作.

經對比測試發現，相比其他幾種網絡地圖，Google Maps 的影像資料最為全面，因此，本文選取Google Maps 影像作為道路數字化選線系統的網絡影像數據源，并研究了相應的下載算法和快速配準算法，實現網絡地理數據本地化，以便服務于基于客戶端/服務器工作模式的道路數字化選線系統.

2 Google Maps 的影像瓦片下載

2.1 Google Maps 的數學原理

2.1.1 地圖投影

地圖投影是指建立地球球面點與地圖平面點之間一一對應關系的數學變換方法. Google Maps采用的地圖投影方式為Web 墨卡托(投影定義EPSG (Eourpean Petroleum Survey Group )3875，投影運算方法代號1024). 假設球面經緯度坐標為(λ，φ)，地圖平面坐標為(E，N)，其投影正算公式為［18］

式中:

E 為東西方向坐標;

N 為南北方向坐標;

λ 為經度，取值范圍為［-π，π］;

φ 為緯度，取值范圍為［- 0. 472 506 27π，0.472 506 27π］;

R 為球體半徑，計算過程中取WGS84 橢球的長半軸半徑6 378 137.0 m.

2.1.2 瓦片索引

Google Maps 瓦片坐標及四叉樹分割示意如圖1 所示.

圖1 Google Maps 瓦片坐標及四叉樹分割示意Fig.1 Schematic diagram of Google Maps tile coordinates and quad tree segmentation

Google Maps 采用瓦片金字塔模型對影像瓦片進行組織與管理，每個瓦片大小為256 ×256 像素.經測試，目前Google Maps 瓦片最高縮放等級n =23，但是并非每個區域都有第23 級影像，最低縮放等級n =0，此時整個地球經過投影和裁剪后得到一塊256 ×256 像素的瓦片.瓦片數據采用四叉樹進行編碼，利用瓦片坐標來標識每個瓦片，瓦片坐標由i 和j 兩個整數組成，由左至右i 依次遞增，由上至下j 依次遞增. 設左上角瓦片的坐標位置為(0，0)，右下角第n 級瓦片的坐標位置為(2n-1，2n-1)，見圖1.

為了由經緯度獲取Google Maps 的瓦片數據，利用經緯度(λ，φ)計算出瓦片坐標(i，j). 當縮放等級為第n 級時，整個投影區域被分為4n塊，i 和j 方向分別均分為2n等份，第n 級瓦片的分辨率(單位:m/像素)為

根據投影坐標系與瓦片坐標系的等比例關系，得到瓦片坐標與投影坐標的計算公式為

將式(3)與式(1)、(2)聯立，得瓦片坐標與球面經緯度坐標為

2.2 瓦片URL 地址分析

Google 提供影像瓦片下載服務的服務器有兩類，即mt 和khm. mt 服務器能提供普通地圖、衛星影像以及包含道路、界線信息的合成地圖，khm 僅提供衛星影像.通過構造URL 進行Google Maps 影像瓦片下載，首先需要了解瓦片URL 各字段的含義. Google 影像瓦片URL 地址格式舉例如下:http://mt1.google.cn/vt/lyrs=h@167＆hl =en＆src= app＆x = 216＆y = 94＆z = 8＆s = Galile(簡稱URL1)和http://khm0.google.com/kh/v =140＆src= app＆x = 216＆y = 94＆z = 8＆s = Galile(簡稱URL2).

以URL1、URL2 為例，對瓦片URL 各關鍵字段含義進行測試，測試結果見表2.

表2 Google 影像瓦片URL 字段含義Tab.2 Definition of URL fields in Google image tiles

明確Google Maps 瓦片URL 各字段含義后，就可以根據下載任務自動計算構造URL 地址，并通過libcurl 庫(應用于客戶端的開源多協議文件傳輸庫)函數進行下載. 需要注意的是，在構造瓦片URL 地址時，語言種類hl 應盡量選擇為英語或者直接缺省，因為測試發現，選擇簡體中文時影像存在地理坐標偏移.

2.3 多線程下載策略

為了平衡用戶請求，Google 的mt 服務器和khm 服務器都各有4 臺.由于高分辨率瓦片數量巨大，為了實現高速下載，宜采用多線程技術.由于不同線程鏈接的服務器可能不同，所以下載速度也存在差異，如果直接將任務均分給n 個線程，則必然使某一線程成為瓶頸. 因此，本文設計了一種“能者多勞”的多線程下載任務分配方式，首先將待下載的瓦片地址計算好，假設有這樣一個“任務池”，將每個瓦片看成1 個“任務”，將計算好的瓦片放入任務池中，當一個線程開始執行或是完成了上一個下載任務時，就會主動向“任務分派模塊”發送任務請求，然后由“任務分派模塊”從“任務池”中獲取任務，如果獲取成功，則將任務信息發送給該線程，并將該“任務”從“任務池”中刪除. 如果“任務池”中沒有任務，則向該進程發送終止運行消息，線程正常退出，不再執行任務，直至最后一個線程正常退出，向系統發送下載完成消息，否則發送異常退出消息，并記錄在下載日志文件中.

在多線程下載過程中，可以不指定構造的瓦片URL 服務器編號，這樣空閑線程會主動尋找最快鏈接上的服務器. 不過這時使用libcurl 庫函數下載瓦片，服務器將拒絕數據請求并回復403 錯誤，需要調用libcurl 庫中curl_easy_setopt()函數對CURLOPT_USERAGENT 選項進行模擬瀏覽器設置，其值可以為空，也可以設置為“User-Agent:

Mozilla/5.0 (compatible;MSIE 8.0;Windows NT 5.1;CIBA)”.使用libcurl 庫的curl_easy_getinfo()函數獲取CURLINFO_RESPONSE_CODE 的返回值，可以判斷是否下載成功，當返回值為200 時，代表下載成功.下載完成后為影像瓦片指定一個唯一索引標記Key，并將其作為該瓦片的文件名. Key的生成規則為Ln_Xi_Yj;n 為縮放等級，i 和j 為瓦片索引坐標.

3 Google Maps 瓦片與高程數據配準

下載的Google Maps 影像瓦片需要與SRTM 高程數據進行配準，并重投影后才能應用于數字化選線系統中.重投影工作相對比較簡單，可以在影像和高程數據配準后，通過GDAL 庫編寫轉換程序自動完成.

本文重點探討Google Maps 影像和SRTM3 高程數據的快速配準算法. 由于Google Maps 使用的是Web 墨卡托投影，這是一種非等長的投影方式，在高緯度南北方向的長度變形非常大，投影網格橫縱比不為1. 而SRTM3 高程數據屬于全球等間隔經緯度劃分的格網DEM，采用的是一種十進制度的投影方式，網格橫縱比為1.其投影公式為

由于Google Maps 影像和SRTM3 采用相同的參考橢球，根據兩者的投影公式可知，兩者經度在橫向上的變化是線性一致的. 將Web 墨卡托投影下的影像變換到十進制度投影下的影像，使緯度在縱向上也呈線性關系，并且保證橫縱比例為1. 因此，根據Web 墨卡托投影特性，需要對影像在縱向上進行非線性變換，變換示意圖如圖2 所示.

通過逐像元重采樣方式進行非線性變換，計算量將非常巨大.為此，借用分治法的思想，把全局影像數據的投影變換問題轉換為局部影像數據的投影變換問題，并在局部范圍內直接采用線性變換方法代替非線性變換. 算法設計思想是，尋找待變換影像需滿足的條件，使線性變換后像素誤差與非線性變換后像素誤差的差值保持在一定的范圍內.

圖2 Web 墨卡托投影到十進制投影的變換示意圖Fig.2 Schematic diagram of transformation from Web Mercator projection to decimal degrees projection

首先定義相對位置C 為影像上某點縱坐標與影像底邊縱坐標之差與整個影像縱向長度的比值，為無量綱量.

設變換前影像上某點P 的相對位置為CP，線性變換后為CPL，非線性變換后為CPN.根據投影原理有

式中:N1、N'1為變換前、后影像底邊的縱坐標;N2、N'2為變換前、后影像頂邊的縱坐標;

NP、N'P為變換前、后影像上點P 的縱坐標;φP為P 點的緯度，

φP=2arctan(exp(xP/R))-π/2;

φ、Δφ 為變換前影像底邊的緯度以及影像的緯度差.

圖3 變換前后影像上點P 相對位置示意圖Fig.3 Schematic diagram of relative locations of point P in the image before and after compression

使用線性變換與非線性變換相對位置間的差值表示兩種變換間的誤差，相對位置誤差為

式(8)對CP求偏導并令其為0，則有

該方程存在解析解，利用Mathematica 求解得:

從上式可知，CP是關于φ、Δφ 的函數，令

CP=C(φ，Δφ)，

將CP帶入式(8)得:

由式(9)可知，相對位置誤差EEC僅是緯度差Δφ 和緯度φ 的函數.

假設變換前影像高度為Hsrc，變換后影像高度為Htar，將EEC換算成以像素為單位的誤差EEP，則有

式中:Htar=Δφ/γresrad，

以弧度為單位計算的影像分辨率.

令

利用Mathematica 繪制該隱式函數的曲線，如圖4所示.

圖4 =0.5 的曲線Fig.4 Curve for=0.5

由圖4 可知，對于Google Maps 影像瓦片，當其處于低緯度區域時，中心緯度小于8°，或者底邊頂邊緯度差絕對值小于0.8°，即在曲線與坐標軸圍成的區域內，可以采用線性變換代替非線性變換.根據Web 墨卡托投影特性，經過赤道線瓦片的緯度差最大，對經過赤道的第N 級瓦片有

令360°/2n＜0.8，則有

取整n =9. 所以對于Google Maps 影像瓦片，縮放等級大于9 級后，即可采用線性變換方法代替逐像元計算的非線性變換方法.由于9 級瓦片的分辨率較低，對于實際工程意義不大，在道路數字化選線系統中，通常下載影像級別在19 級及以上，因此，完全可以采用線性變換以提高影像配準速度.

4 實驗與系統驗證

本文選取緯度范圍21°04'09″N ～21°04'49″N、經度范圍101°37'50″E ～101°38'58″E 作為測試區域，進行影像瓦片下載. 實驗環境主要配置如下:CPU:i5 3470;內存2G，DDR3;顯卡:NVIDIA GeForce 605，顯存512M;操作系統:windows 7 系統.下載影像瓦片縮放等級為第20 級，共包含瓦片數44 ×31 塊，采用多線程下載，耗時38.65 s.該區域的DEM 數據由SRTM3 DEM 重采樣生成. 將影像瓦片進行線性變換后，拼接并重投影至UTM(universal transverse mercator)投影下，生成可用于道路數字化選線系統的DEM 和DOM 數據，如圖5所示.將生成的DOM 數據與該處既有的0.2 m 分辨率的高清航飛影像資料進行局部對比，如圖6所示.

圖5 獲取的實驗區域DEM、DOM 數據Fig.5 DEM and DOM data obtained from the test region

圖6 影像局部對比Fig.6 Local contract of image

圖6 底層為既有高清影像，疊加的為Google Maps 下載并配準后的影像. 從圖6 道路接邊處可以看出，下載的影像配準效果良好，并與既有的0.2 m 分辨率影像的清晰度相當.

最后，將獲取的實驗區DEM 和DOM 數據輸入道路數字化選線系統中，進行三維地理環境創建.并在該三維地理環境下進行實驗路線選擇、實驗線三維實時建模、路線漫游，效果如圖7 所示.

圖7 網絡地理信息在道路數字化選線系統中的應用Fig.7 Application of web geographic information in the digital highway location system

5 結 論

本文在綜合分析當前開放網絡地理信息資源的基礎上，選用Google Maps 衛星影像數據和SRTM3 DEM、ASTER GDEM 作為道路數字化選線系統的數字地形網絡獲取數據源，基于Google Maps 的數學投影原理和多線程技術，研究了影像瓦片的快速下載與配準算法，得出如下結論:

(1)通過構造瓦片URL 地址，基于libcurl 開源庫，采用“能者多勞”的多線程下載策略下載Google Maps 影像瓦片，具有獲取方便、速度快捷、自動化程度高等特點.

(2)Google Maps 影像與SRTM3 DEM 配準時，當影像瓦片等級大于9 級后，可以采用線性變換代替非線性變換.

(3)將獲取的DEM 和DOM 應用于道路前期規劃數字化選線，可以減少外業工作量和工作難度，降低企業設計成本. 構建的三維虛擬地理環境逼真，滿足道路前期規劃階段開展數字化選線設計的需要.

本文以網絡地理信息服務于道路數字化選線系統為例進行研究.同理，在鐵路、電力等行業的數字化選線中，也可以應用本文的方法. 隨著遙感衛星影像質量的不斷提高，網絡地理信息資源將更多地應用于工程設計.

致謝:本文工作得到西南交通大學道路工程四川省重點實驗室開放研究基金(LHTE001201101)和西南交通大學高速鐵路線路工程教育部重點實驗室開放研究基金(2010-HRE-02)的資助.

［1］ 易思蓉，朱穎，許佑頂. 鐵路線路BIM 與數字化選線技術［M］. 北京:中國鐵道出版社，2014:1-2.

［2］ CHANDER G，MARKHAM B L，HELDER D L.Summary of current radiometric calibration coefficients for Landsat MSS，TM，ETM +，and EO-1 ALI sensors［J］. Remote Sensing of Environment，2009，113(5):893-903.

［3］ ROYD P，WULDER M A，LOVELANDT R，et al.Landsat-8:science and product vision for terrestrial global change research［J］. Remote Sensing of Environment，2014，145:154-172.

［4］ 蔣勇. 各種衛星影像數據在鐵路勘測中的應用［J］.鐵道勘察，2013(5):16-18.JIANG Yong. All kinds of satellite image data application in railway survey［J］. Railway Investigation and Surveying，2013(5):16-18.

［5］ 張朝忙，劉慶生，劉高煥，等. STRM3 與ASTER GDEM 數據處理及應用進展［J］. 地理與地理信息科學，2012(5):29-34.ZHANG Chaomang，LIU Qingsheng，LIU Gaohuan，et al. Data processing and application progress of SRTM3 and ASTER GDEM［J］. Geography and Geo-Information Science，2012(5):29-34.

［6］ BAILEY J E，SELF S，WOOLLER L K，et al.Discrimination of fluvial and eolian features on large ignimbrite sheets around La Pacana Caldera，Chile，using landsat and SRTM-derived DEM［J］. Remote Sensing of Environment，2007，108(1):24-41.

［7］ SU Y F， SLOTTOW J， MOZES A. Distributing proprietary geographic data on the World Wide Web:University of California at Los Angeles CIS database and map server［J］. Computers and Geosciences，2000，26(7):741-749.

［8］ BLOWERJ D，HAINES K，SANTOKHEE A，et al.GODIVA2:interactive visualization of environmental data on the web［J］. Philosophical Transactions of the Royal Society A: Mathematical， Physical and Engineering Sciences，2009，367(1890):1035-1039.

［9］ 王子明. Google Earth 在公路工程可行性研究中的應用［J］. 公路交通技術，2008(5):4-6.WANG Ziming. Application of google earth in feasibility study of expressway project［J］. Technology of Highway and Transport，2008(5):4-6.

［10］ 朱穎，蒲浩，劉江濤，等. 基于數字地球的鐵路三維空間選線技術研究［J］. 鐵道工程學報，2009(7):33-37.ZHU Ying，PU Hao，LIU Jiangtao，et al. Research on the digital earth-based 3D spatial technology for railway route selection［J］. Journal of Railway Engineering Society，2009(7):33-37.

［11］ 易共才，王彥軍，高宏. Google Earth 在公路工程中的應用研究［J］. 中外公路，2008(1):1-4.

［12］ 王一波，邵偉偉，羅新宇. Google Earth 數據精度分析及在鐵路選線設計中的應用［J］. 鐵道勘察，2010(5):68-71.WANG Yibo，SHAO Weiwei，LUO Xinyu. Analysis on accuracy of Google Earth data and its applications in optimized design on railway routes［J］. Railway Investigation and Surveying. 2010(5):68-71.［13］高山. 基于SRTM DEM，ASTER GDEM 地貌特征分析與鐵路選線［J］. 鐵道工程學報，2012(9):1-6.GAO Shan. Analysis of topographic feature with SRTM DEM and ASTER GDEM data and railway alignment［J］. Journal of Railway Engineering Society，2012(9):1-6.

［14］ 寇曼曼，王勤忠，譚同德. Google Map 數字柵格地圖算法及應用［J］. 計算機技術與發展，2012(4):204-206.KOU Manman， WANG Qinzhong， TAN Tongde.Research on goolge map algorithm and application［J］.Computer Technology and Development，2012 (4):204-206.

［15］ 崔金紅，王旭. Google 地圖算法研究及實現［J］. 計算機科學，2007(11):193-195.CUI Jinhong，WANG Xu. Research on goolge map algorithm and implementation［J］. Computer Science，2007(11):193-195.

［16］ NIKOLAKOPOULOS K G，KAMARATAKIS E K，CHRYSOULAKIS N. SRTM vs ASTER Elevation products， comparison for two regions in Crete，Greece［J］. International Journal of Remote Sensing，2006，27(21):4819-4838.

［17］ 杜小平，郭華東，范湘濤，等. 基于ICESat/GLAS 數據的中國典型區域SRTM 與ASTER GDEM 高程精度評價［J］. 地球科學(中國地質大學學報)，2013，38(4):887-897.DU Xiaoping，GUO Huadong，FAN Xiangtao，et al.Vertical accuracy assessment of SRTM and ASTER GDEM over typical regions of China using ICESat/GLAS［J］. Earth Science(Journal of China University of Geosciences)，2013，38(4):887-897.

［18］ International Association of Oil and Gas Producers(IOGP). Coordinate conversions and transformations including formulas， guidance note 7-2［EB/OL］.［2015-04-26］. http://www. epsg. org/Guidance Note.