Python编程及Google Earth平台在地球物理勘探中的应用

王丹鹤，吕玉增，2，孙 拓

(1.桂林理工大学 地球科学学院，广西 桂林 541006；2.广西隐伏金属矿产勘查重点实验室，广西 桂林 541006)

关键字：地球物理勘探；Python；KML；Google Earth；地标

0 引言

在地球物理勘探工作中，常常会由于野外踏勘不足，造成测点、测线布置不合理、野外工作无法开展等情况，如测线经过厂矿，峭壁，水域，居民区等，迫使调整工作方案，致使预算增加，严重影响项目实施。近年来，随着地理信息技术的发展，卫星图片的分辨率，定位精度都在不断提高。因此，若能利用Google Earth的卫星影像和大数据平台，在室内对勘探区域进行云踏勘并指导地球物探勘探的全过程，必将提高地球物探勘探的工作效率，节省勘探成本。

近年来，围绕Google Earth平台在地球物理勘探中的应用，许多学者都做了相关的讨论。史来亮[1](2015)，史晓亮[2](2014)、罗文刚[3](2011)、李春芬[4](2011)、关玉东[5](2008)等人都是将Google Earth应用于地震勘探，取得了很好的效果；苗放[6](2006)介绍了Google Earth在GIS (Geographic information system)方面的应用；王路希[8](2015)、雷清[9](2019)、李良[10](2014)、叶宝莹[11](2017)分别讨论了Python在开源GIS、大地电磁测深、测井、矿产勘查[12]、环境科学[13]、地球化学[14]及地质信息提取方面的应用；国外学者运用Google Earth在地质与地球物理建模[15]以及地理信息查询[16]等。此外王艳[7](2009)在基于VC++的Google Earth KML地标文件的自动生成及应用一文中介绍了KML文件，提出了利用KML文件实现了点、线的精确定位与显示，没有详细讨论如何自动化生成。部分商业软件也能实现从经纬度坐标转换成KML文件，但是并没有开放源代码版本的出现。在现在这个阶段，开放源代码实现版本的出现尤为重要。

本文重点讨论KML(Keyhole Markup Language)文件，并利用Python编程，使用simplekml模块，开源实现自动化生成KML文件，配合Google Earth，实现对物探工作的点、线和区的设计优化，辅助数据处理解释。Python的优势就在于算法简单明了，可读性高，易于后期维护利用。

1 KML简介

KML全称为Keyhole Markup Language，是基于可扩展标记语言(Extensible Markup Language，XML)语法标准的一种标记语言，采用标记结构，含有嵌套的元素和属性。根据KML语言编写的文件则为KML文件，格式同样采用XML的文件格式，用于显示地理数据(包括点、线、面、多边形、多面体以及模型等)。后来Google将KML提交给OGC (Open Geospatial Consortium)组织，由OGC继续开发和维护。在地球物理勘探中，经常使用到KML语言中的点、线、区块等地标文件。

KML文件主要由XML标头，KML命名空间声明，包含元素的地标对象组成。包含元素的地标对象主要包括：用作地标标签的名称，附着到地标的“提示框”中显示的说明，指定地标在地球表面的位置的点(经度、纬度和高度等)。其中重要的标签有：表示点的标签，表示线的标签，为多边形标签，坐标的标签，是线型标签，是颜色标签，线宽标签。

2 Python生成KML地标文件

Python是一种广泛使用的解释型、通用型编程语言，由吉多·范罗苏姆创造。相比于C++或Java等，Python更强调代码的可读性和简洁的语法，能让开发者用更少的代码表达思想并实现。Python的另一个优点是其解释器几乎可以在所有的操作系统中运行，程序的可移植性好。

KML默认的做坐标格式是WGS 84系统下的经纬度坐标。对于北京54系统或西安80系统下的公里网格坐标而言，首先需要转换成WGS 84系统经纬度坐标，然后生成KML地标文件。接下来重点介绍如何用Python来创建KML格式的地标文件。

2.1 创建点地标KML文件

若已知若干个点在WGS 84系统下经纬度坐标，利用Python的“import”命令可方便创建KML文件。具体实现代码如下：

已知信息Python的代码实现#点号及其经纬度坐标points=[('1',#1经度,#1纬度),('2',#2经度,#2纬度),('3',#3经度,#3纬度),('4',#4经度,#4纬度),('5',#5经度,#5纬度),('6',#6经度,#6纬度),]kml=simplekml.Kml(open=1)#导入Python模块#绘制点位forindex,lon,latinpoints:pnt=kml.newpoint()pnt.name=indexpnt.coords=[(lon,lat)]kml.save("test_points.kml")#保存KML文件

2.2 创建线地标KML文件

创建线地标的时候，关键是要给出测线起止点和转折点的坐标。代码如下：

已知信息Python的代码实现#起止点经纬度坐标line_points=[(#1经度,#1纬度),(#2经度,#2纬度)]importsimplekml#导入Python模块kml=simplekml.Kml(open=1)linestring=kml.newlinestring(name='Aline')linestring.coords=line_points#设置线的颜色linestring.style.linestyle.color='ff0000ff'#设置线的宽度linestring.style.linestyle.width=10kml.save("test_line.kml")#保存KML文件

2.3 创建区地标KML文件

创建区地标的时候，需要区域边界的经纬度坐标。具体代码如下：

已知信息Python的代码实现#区域边界点经纬度坐标pol_points=[(#1经度,#1纬度), (#2经度,#2纬度), (#3经度,#3纬度), (#4经度,#4纬度), (#5经度,#5纬度), (#6经度,#6纬度), (#7经度,#7纬度), (#8经度,#8纬度), (#9经度,#9纬度), (#10经度,#10纬度), (#11经度,#11纬度), (#12经度,#12纬度)]importsimplekmlkml=simplekml.Kml()pol=kml.newpolygon(name='APolygon')#多边形边界点的坐标pol.outerboundaryis=pol_pointspol.style.linestyle.color=simplekml.Color.yellowpol.style.linestyle.width=5pol.style.polystyle.color=simplekml.Color.changealphaint(100,simplekml.Color.green)kml.save('test_polygon.kml')

生成区地标KML文件见图1。

图1 点、线、区地标的KML文件图Fig.1 KML file diagram of points，line and polygon landmarks

3 初步应用

3.1 优化勘查方案设计

当勘查区的点、线和测区的初步工作方案后，可依据上述方法生成多个KML地标文件，将生成的KML文件导入Google Earth，就可以对工区的地形、地貌特征，有一个全面的了解。如遇到矿区，水域，可以提前调整测线及测点的位置，科学合理地调整无法工作的区域。

在项目方案设计和实施阶段，利用Google Earth，基本可掌握测区的地形、地貌和构筑物等地理信息。同时，可初步对选定的物探勘探方法、方案进行优化改进。下面举例说明。

案例1：图2为某实际测区布设物探激电测深的部分测点正好落在构筑物上，通过“地标文件+Google Earth”提前预知后，就可对这些测点进行合理地规划，尽可能的避开构筑物。

图2 点地标卫星图Fig.2 Satellite map with point landmarks

案例2：图3为某测区的一条激电测深勘探线，整个测线跨过了多条公路、村庄和陡峭灰岩区，提前准确掌握这些信息，可在物探实际观测中合理地调整观测方案，高效实施，以降低野外观测成本。图4显示激电测线正好穿过一个采石场(黄线区域)，建议做弃点处理。

图3 线地标卫星图Fig.3 Satellite map with a line landmark

图4 区地标卫星图Fig.4 Satellite map with a polygonal landmark

3.2 辅助数据处理和解释

在物探数据处理阶段，往往需要测线所有测点的准确高程数据，便于异常的查证和解释。通过调用Google的应用程序接口(Application Programming Interface，API)，利用Python可准确获取测线测点的高程数据。实现流程如下：先打开测点文件，读取经纬度；然后构造超文本传输安全协议(Hyper Text Transfer Protocol Secure，https)请求；返回JavaScript对象表示法(JavaScript Object Notation，JSON)格式的文件；最后解析JSON文件，读取相应点的高程，绘制测线地形图。

下面以激电测深的一条测线为例，具体实现代码如下：

Python的代码实现importrequests#导入请求包importmatplotlib.pyplotasplt#导入绘图包points=[]lat=[] foriinrange(len(j['results'])): outputstr="{:.6f}{:.6f} {:.6f}\n".format (j["results"][i]['location']['lng'], j["results"][i]['location']['lat'], j["results"][i]['elevation'])

lon=[]#打开测点文件withopen('points.txt',encoding='utf-8-sig')asf:forlineinf: points.append(line.replace('\n',''))forpointinpoints: s=point.split(',') lat.append(s[1]) lon.append(s[0])npoints=[[x,y+'|']forx,yinzip(lat,lon)]points_loc=[]foriinnpoints:points_loc.append(i[0]+','+i[1])loc_str=''.join(points_loc)#API密钥my_key='&key=yourkey'#网址头部url='https://maps.googleapis.com/maps/api/elevation/json?locations='#链接请求r=requests.get(url+loc_str[0:len(loc_str)-1]+my_key)j=json.loads(r.text)#解析json数据#存储经纬度高程数据withopen("Elevations.dat","w")asf_w: f_w.writelines(outputstr)f_w.close()#绘图部分lle=[]ele=[]di=[]withopen("Elevations.dat","r")asf_r:forlineinf_r: lle.append(line)foriinrange(len(lle)):ele.append(float(lle[i].split()[2]))forlinrange(len(ele)):di.append(l∗50)fig=plt.figure(figsize=(16,4))axes=plt.axes()plt.title('测线地形地形起伏图',fontsize=20)axes.set_ylim([0,1600])axes.set_xlim([0,6500])axes.spines['top'].set_visible(False)plt.plot(di,ele,color='black')plt.rcParams['font.sans-serif']=['MicrosoftYaHei']plt.text(5500,1500,'地形线',fontsize=13)plt.xlabel('距离/m',fontsize=15)plt.ylabel('高程/m',fontsize=15)plt.minorticks_on()plt.savefig("topography.svg",dpi=600)plt.show()

该地形的高程数据可以用来对物探成果图进行白化处理。利用Python生成的KML文件，导入Google Earth，直观明了地查看地形地貌(图6)，辅助技术人员进行异常识别和分析。甚至可以获取测区的三维地形数据，进行三维的地形建模和正演模拟，进而对物探异常进行准确识别和综合解释。

图5 测线的地形示意图Fig.5 Topographic sketch map of the survey line

图6 带地形的区线地标图Fig.6 Topographical map with polygonal and line landmarks

4 结论

本文详细地展示了利用Python自动化生成点、线和区域的KML文件和获取测线高程数据的过程，算法简洁明了。结合Google Earth使用，可以有效地优化勘查设计方案，若测线测点遇到不可测区域，能够提前查看，做到及时调整。在勘探工作中，Python和Google Earth的联合使用可大大节省踏勘、设计等时间和生产成本，提高生产效率。另一方面，可以便利地获取高程等数据，辅助数据处理和解释，随着地形、地貌和地质信息的加入，使得整个地球物理勘探过程变得高效，有利于综合解释。

借助图形识别技术，联合程序编制、虚拟仿真、导航等技术并将其引入到地球物理勘探与解释中，提供了一个实现复杂综合信息情况下云踏勘和云设计的便捷方案。