基于网络地理信息服务的选线数字地形获取方法

聂良涛， 易思蓉， 李 阳

(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.