#
大家好！这几天微博上最热的话题之一就是汤加火山的喷发。网友们一方面为汤加人民和几千位华人华侨担忧祈福，另一方面也开始关心本次火山喷发对各地气候的影响。
#
不过无论看知乎还是翻公众号，我们看到的大多都是针锋相对的两种观点：要么说本次喷发会导致全球降温、“无夏之年”，要么说它对全球气候的影响微乎其微。
当然，作为接受过科学教育的理性群众，我们最信任的自然是专业学者的解释，所以看了下面这篇采访，杨老师安心了许多：
不过孙劭老师在采访中使用的论据引起了我的兴趣：“因为本次喷发的二氧化硫与菲律宾皮纳图博火山不同，所以不会产生后者那样的气候影响”。那么反过来说，如果我们能够找到一个与本次喷发相似的火山案例，是不是就可以用后者预测本次喷发的影响了呢？
换句话说，如果我们有一个记录历次火山喷发各方面特征的数据集，又有本次汤加火山喷发的对应数据，就可以计算本次喷发与各次历史喷发间的相似性，从而找到最相似的那一次。然后只要看看那一次喷发后我国乃至全球气候有什么变化，就可以合理猜测本次喷发的后果。
学过《Python实战篇 数据分析》的同学可能会觉得有点眼熟：我们上周发布的第53回内容 —— “K-近邻算法” 不就是计算相似性的么？没错，KNN的相似性计算方法正好可以用在这里，即使我们这个任务中缺少标签、更接近于“聚类”而不是“分类”问题。
想到这一点，杨老师心情顿时激动起来！要知道，当年我上高中时深受徐迟先生《祁连山下》的影响，一度坚定报考地质专业追山赶海为国探矿，直到后来老师告诉我高度近视不宜就读，才打消了这个念头。所以今天突然有机会冒充一次地质学(min)家(ke)，岂不快哉？
说干就干！第一件事情就是寻找历史上所有火山的喷发数据。经过搜索，最后选择了史密森学会的在线版火山历史数据库，下载得到9941条已经确认 (Confirmed Eruption) 的火山喷发记录：
不过很遗憾的一点是，这个数据集中不包含二氧化硫喷出量等等对气候影响巨大的指标。原因当然可以理解，毕竟远古时代有谁能去测量二氧化硫呢？
想起AI专家们常说的那句话：数据和特征决定了机器学习的上限。看来在找到更全面的数据集之前，我们只能压低期望。好在数据集还是提供了几个非常重要的特征，比如地理位置的经纬度、喷发烈度的VEI评级，还有喷发的起止日期。于是稍加整理，选出上述几个特征都不为空、喷发时间在公元元年以后的数据，最终杨老师整理出3986条记录（倒数第2条就是摧毁庞贝的那次维苏威火山喷发）：
不过随之又发现一个问题：很多喷发记录的“起始月”或“起始日”等数值都是0，而不是正常的1-12月或者1-31日。个人猜测，大概是因为古代典籍中只记录了这些喷发的年份，而现在考古和地质证据也无法精确测定月份日期，所以只能填写为0。
显然，这就是机器学习中典型的“缺失值”问题，我们可以选择删除它们或者选择一个合理的方式对其修改。考虑到数据难得、删之难过，这里我们暂时把所有为0的月份改为年中6月、把所有为0的日期都改为月中15日。整个操作使用pandas完成，处理后数据如下：
理清了数据，接下来就要研究这几个有限的特征应该怎样充分利用。基于从自己给孩子购买的科普读物中学到的全部火山知识，杨老师做了以下几个假设：
1
地理距离越近的火山喷发，对气候影响的相似度越高。原因：相邻地理区域的地貌、海流和大气环流规律接近；
2
季节越接近的火山喷发，对气候影响的相似度越高。原因：相近季节意味着气压、季风等大气传播条件接近；
3
喷发烈度（VEI评级）越接近的火山喷发，对气候影响的相似度越高。原因不用解释。
4
持续时间越接近的火山喷发，对气候影响的相似度越高。原因同（3）。
细究起来，上面的假设当然还有很多瑕疵。比如“地理距离近”并不代表二者就在同一区域，哪怕中间只隔着一条山脉，也会导致气候规律的巨大差异。不过作为业(xian)余(de)爱(wu)好(liao)，我们这一版就先不考虑那么多了。
但是有一个问题还是需要深入考虑的：喷发烈度评级 VEI 使用的是对数量表：2级的实际烈度是1级的10倍、3级烈度是2级的10倍 …… 。因此，到底采用1~8的评级数字作为比较依据、还是采用 1~100000000 的实际烈度数字，确实是一个问题。杨老师不具备相关知识，所以还是不要妄改数据，暂时就用 VEI 评级数字吧。
此外我们还要使用python的datetime模块计算出火山持续时间（注意不能使用pandas自带的datetime，因为后者无法计算古代年份），并且考虑到0级喷发记录对本次汤加火山喷发参考意义不大，所以删除这部分数据，于是得到训练数据集：
接下来，再用本次汤加火山的数据创建一个数据集。来不及详细搜索，就把火山的经纬度直接取汤加王国的大致经纬度，烈度假设6级、持续时间3天：
万事俱备、只欠东风！现在我们只要编写几个函数用来计算火山数据之间的距离，就可以开始预测。闲言少叙、先上代码：
这段代码中定义了三个函数。前两个分别用于计算两条喷发记录的直线地图距离（以经纬度为单位、不考虑地球弧度）和季节差异（即两次喷发的月份相隔多少），而最后一个则是计算所有特征的欧氏距离。由于这些特征之间数据范围各异，因此每一种距离都采用了最大值方法进行归一化。同时考虑到在我们仅有的这几个特征中，喷发烈度是最能影响天气的一个，因此给他赋予了70%的权重。
之所以要自己编写距离计算函数，是因为经纬度距离与季节差异并不是简单的减法操作。比如：东经180度和西经180度之间的差距，如果直接使用 +180 - (-180) 就会得出二者相差360度；然而实际上，这两者是同一条经线、距离为 0。同样的道理，11月和2月的间隔，并不是 11 - 2 = 9个月，而是 3 个月，所以都需要我们单独编写计算规则。
当然，上面这样的算法也有问题，因为经纬度距离是按照平面地图坐标系计算，并没有考虑到地球是一个球体，所以对于靠近南北极的高纬度地区火山，计算结果会与真实距离存在较大偏差，因此更精确的计算方法应该是使用三角函数。这个改进就留给感兴趣的同学自己搞定吧。
距离算法搞定，相似度就近在眼前了！我们使用pandas的apply函数，将汤加数据与全部历史数据逐一比较，计算出每个历史喷发记录与本次喷发的距离，放到“汤加火山相似度”列中，数字越小越相似。然后按此排序，选出最接近的10次喷发记录如下：
所以，根据著名非主流火山学家Python-Yang的测算，本次火山喷发与历史上1600年2月秘鲁“于埃纳普蒂纳火山”喷发最相似。其时我国正是万历二十八年，杨老师翻开明史，倒是没找到什么特别的自然灾害记录，仍然是往年常有的秋旱地震等等。但网上查看资料，确实有评论说1600火山喷发影响了全球气温，我国出现了桃花晚开的现象。不过考虑到本数据集中完全没有二氧化硫数据、而1600年秘鲁火山却喷发了大量的二氧化硫，因此上述模型不足为凭，仅算是我们做一个小小的数据分析练习吧。
其实即使作为数据分析练习，这个模型也还只是最粗浅的第一步，有很多问题值得感兴趣的同学深入思考。比如“经纬度距离”并不能反映真实地理距离、更不能代表区域相似性；南半球/北半球季节差异应该得到体现；地震烈度应该有更合理的表示方法；各种特征的权重设定需要客观依据 …… 而最大的问题，则是这个数据集缺少一个真正的目标值，也就是每次喷发对全球气候的真实影响指数，因此无法像有监督学习那样对我们模型进行验证。
但是不用想那么多了，谁让我们是民科呢 ……
注：本文使用的火山喷发记录数据集，与本文一起上传到官网 www.ukoedu.com 《全民一起玩 Python 实战篇 之 数据分析》第53回 “参考资料” 栏目中，有兴趣的同学可以下载练习。
end
长按前往
杨洋老师课程合集
