前言
1、分析的必要性
地震，作为自然界中最为猛烈的灾害之一，其突发性和破坏力往往给人类社会带来巨大的生命财产损失。在历史的长河中，无数次的地震记录不仅见证了地球的脉动，也深刻影响了人类文明的进程。因此，对地震附近区域进行历史地震分析，其必要性不言而喻。
首先，历史地震分析是预测未来地震风险的重要依据。地震活动具有一定的周期性和复发性，通过对某一地区历史地震资料的系统梳理和深入研究，科学家能够揭示该区域地震活动的时空分布规律、震级大小、震源深度等特征，进而评估未来地震发生的可能性及其潜在影响范围。这种基于历史数据的地震预测方法，虽然无法精确到具体的时间和地点，但能够为防灾减灾工作提供宝贵的参考和指导。
其次，历史地震分析有助于提升地震应急响应能力。通过回顾历史上地震发生时的社会反应、救援行动、灾害损失等情况，我们可以总结经验教训，优化应急预案，提高应对地震灾害的效率和效果。例如，分析历史地震中人员伤亡的主要原因，可以指导我们在未来的地震灾害中更加注重建筑结构的抗震性能、公众的地震安全教育和应急演练等方面的工作。
再者，历史地震分析对于城市规划和建设具有重要意义。地震是城市规划中不可忽视的自然因素之一。通过对历史地震资料的详细分析，我们可以了解该地区的地质构造、地震活动背景以及潜在的地震风险，为城市选址、布局、建筑设计等提供科学依据。例如，在地震多发区域，应尽量避免建设高密度居住区、重要基础设施等，以降低地震灾害可能造成的损失。
2、分析的紧迫性
面对日益严峻的地震灾害形势，加强地震附近区域的历史地震分析工作显得尤为紧迫。
一方面，随着全球气候变化和地壳运动的加速，地震活动的强度和频率可能呈现上升趋势。这意味着我们需要更加深入地了解地震的成因机制、活动规律以及潜在影响，以便更好地应对未来可能发生的地震灾害。而历史地震分析作为揭示这些规律和特征的重要手段之一，其紧迫性不言而喻。
另一方面，随着城市化进程的加快和人口密度的增加，地震灾害对人类社会的影响日益凸显。一旦地震发生，其造成的损失往往难以估量。因此，我们需要通过历史地震分析来识别潜在的地震风险区域，制定科学的防灾减灾措施，减少地震灾害对人类社会的影响。这种紧迫性不仅体现在对生命财产安全的保护上，更体现在对社会稳定和经济可持续发展的维护上。
综上所述，地震附近历史地震分析的必要性和紧迫性不容忽视。我们应该高度重视这项工作，加强科研投入和人才培养，推动历史地震分析技术的不断创新和发展，为人类社会提供更加全面、准确的地震灾害预测和应对能力。
本文即在此背景下诞生，文章主要讲解如何使用Java开发语言，结合成熟的SpringBoot框架来加速WebGIS应用的开发，首先讲解在数据库中的相关设计以及实现思路，然后重点讲解如何使用Java来开发地震震中附近历史地震的应用。通过本文的相关知识介绍，相信您一定会掌握如何进行附近历史地震分析有一个深刻的认识，同时抛砖引玉，对其他的类似的WebGIS应用建设提供建设参考。
一、数据库物理模型及空间分析实现
现代的系统软件建设，大多是围绕着数据库来实现的。数据库（请注意这里的数据库可不是指狭义的关系型数据库）指的是所有可以用来存储数据的系统。对于关系型数据库和非关系统数据的知识，不属于博文的重点内容，因此不打算在此详细叙述。为了让大家对空间分析有一个具象的认识，因此这里还是把涉及的相关表进行简单介绍，其次介绍在时空数据库中如何实现相关的需求。
1、数据库物理模型
在这里我们采用与物理世界一一映射的模式来构建时空数据表结构，其实主要的表在前面的相关博客中有所介绍。为方便第一次阅读本博客的朋友也了解相关的表设计。这里再次将其物理表结构罗列出来。
其物理建表语句如下所示（需要注意的是，表结构需要在PostGIS数据库中运行，如果本地或者服务器没有PostGIS，还请先安装相对应的数据库）：
CREATE TABLE "public"."biz_ceic_earthquake" (
"pk_id" int8 NOT NULL,
"auto_flag" varchar(30) COLLATE "pg_catalog"."default",
"cata_id" varchar(30) COLLATE "pg_catalog"."default",
"cata_type" varchar(30) COLLATE "pg_catalog"."default",
"epi_depth" numeric(11,8),
"epi_lat" varchar(15) COLLATE "pg_catalog"."default",
"epi_lon" varchar(15) COLLATE "pg_catalog"."default",
"eq_cata_type" varchar(30) COLLATE "pg_catalog"."default",
"eq_type" varchar(30) COLLATE "pg_catalog"."default",
"is_del" varchar(6) COLLATE "pg_catalog"."default",
"location_c" varchar(255) COLLATE "pg_catalog"."default",
"location_s" varchar(100) COLLATE "pg_catalog"."default",
"loc_stn" varchar(20) COLLATE "pg_catalog"."default",
"m" varchar(10) COLLATE "pg_catalog"."default",
"mmb" varchar(10) COLLATE "pg_catalog"."default",
"mmb2" varchar(10) COLLATE "pg_catalog"."default",
"mml" varchar(10) COLLATE "pg_catalog"."default",
"mms" varchar(10) COLLATE "pg_catalog"."default",
"mms7" varchar(10) COLLATE "pg_catalog"."default",
"new_did" varchar(16) COLLATE "pg_catalog"."default",
"o_time" timestamp(6),
"o_time_fra" varchar(10) COLLATE "pg_catalog"."default",
"save_time" timestamp(6),
"sum_stn" varchar(10) COLLATE "pg_catalog"."default",
"sync_time" timestamp(6),
"epi_id" varchar(10) COLLATE "pg_catalog"."default",
"geom" "public"."geometry",
CONSTRAINT "pk_biz_ceic_earthquake" PRIMARY KEY ("pk_id")
);
CREATE INDEX "idx_biz_ceic_earthquake_eqidept" ON "public"."biz_ceic_earthquake" USING btree (
"epi_depth" "pg_catalog"."numeric_ops" ASC NULLS LAST
);
CREATE INDEX "idx_biz_ceic_earthquake_geom" ON "public"."biz_ceic_earthquake" USING gist (
"geom" "public"."gist_geometry_ops_2d"
);
CREATE INDEX "idx_biz_ceic_earthquake_m" ON "public"."biz_ceic_earthquake" USING btree (
"m" COLLATE "pg_catalog"."default" "pg_catalog"."text_ops" ASC NULLS LAST
);
COMMENT ON TABLE "public"."biz_ceic_earthquake" IS '中国地震台网信息表';2、空间数据库中的空间查询分析
在这里我们的需求场景是，点击一个地震发生点，我们需要在四个地图上分别展示近年3级以上的地震分布、近年5级以上的地震分布、近3年3级以上的地震分布、近3年5级以上的地震分布。除了有震级的一个查询和时间的顾虑，其实这两个查询条件都不是空间分析，而震中100*100公里的矩形范围面则是标准的空间分析查询。
为了实现上面的需求，我们将按照以下的步骤来进行实施：
第一步：首先根据地震的id查询地震的经纬度信息（实现的方法有两种，第一种是从geometry字段中获取这样就需要用到st_x和st_y这两个函数，第二种就是直接从原始表格中的经纬度字段中获取）。
这一步的关键sql如下：
select st_x(t.geom) lon,st_y(t.geom) lat from biz_ceic_earthquake t where t.pk_id = ?第二步：在上一步获取的经纬度位置后，我们来构建查询范围面。与常见的缓冲区查询（即st_buffer为例，这是一个圆形缓冲）不同，这里我们需要构建一个矩形范围面，我们需要采用PostGIS当中的空间函数： ST_MakeEnvelope。这一步生成范围面，不仅要在前端界面上展示出来，同时要在第三步的空间分析中作为空间查询参数。
这一步的关键sql如下：
select st_asgeojson(ST_MakeEnvelope( ? - (? / (111.32 * cos(radians(?)))), ? - (? / 111.32),? + (? / (111.32 * cos(radians(?)))),? + (? / 111.32), st_srid(t1.geom) )) as geomJson from biz_ceic_earthquake t1 where t1.pk_id = ?
84.18(BigDecimal), 100(BigDecimal), 40.78(BigDecimal), 40.78(BigDecimal), 100(BigDecimal), 84.18(BigDecimal), 100(BigDecimal), 40.78(BigDecimal), 40.78(BigDecimal), 100(BigDecimal), 1813908285380702210(Long)这里传入的参数包括上一步求解的经纬度，还有范围：如100公里和震发点的主键ID。为了便于前端展示，我们再使用st_asgeojson函数将geometry转换成前端可读的数据格式。
第三步：根据面进行数据过滤。这是最关键的一步，通过这一步。我们将上面生成的空间查询指定范围内的对应震级的地震信息。这一步的关键sql如下：
select t1.* from biz_ceic_earthquake t1 WHERE t1.m ::numeric > ? and t1.o_time >= ? and t1.pk_id <> ? and st_contains( ST_MakeEnvelope( ? - (? / (111.32 * cos(radians(?)))), ? - (? / 111.32),? + (? / (111.32 * cos(radians(?)))), ? + (? / 111.32),st_srid(t1.geom) ),t1.geom)
3.0(BigDecimal), 2019-07-21 22:23:54.581(Timestamp), 1813908285380702210(Long), 84.18(BigDecimal), 100(BigDecimal), 40.78(BigDecimal), 40.78(BigDecimal), 100(BigDecimal), 84.18(BigDecimal), 100(BigDecimal), 40.78(BigDecimal), 40.78(BigDecimal), 100(BigDecimal)二、Java后台程序开发
在上一节中重点对空间数据库的物理表及相关的空间查询分析sql进行了讲解。从本节开始，我们将重点介绍如何在后台结合SpringBoot框架进行后台程序的设计和开发。程序的开发架构依然采用MVC开发模式。这个开发模式的好处是简单易用好实现，复杂度低。
1、模型层设计
前面已经将数据库的物理结构模型提供出来，这里对于java的对应实体类，就不提供源代码了，大家根据数据库结构可以创建对应的实体对象。这里仅将上一节中的sql实现在Mapper层的设计进行说明。
package com.yelang.project.extend.earthquake.mapper;
import java.math.BigDecimal;
import java.util.Date;
import java.util.List;
import org.apache.ibatis.annotations.Param;
import org.apache.ibatis.annotations.Select;
import com.baomidou.mybatisplus.core.mapper.BaseMapper;
import com.yelang.project.extend.earthquake.domain.crawler.CeicEarthquake;
import com.yelang.project.extend.earthquake.domain.crawler.CeicLonLatVO;
public interface CeicEarthquakeMapper extends BaseMapper<CeicEarthquake>{
static final String FIND_ENTITY_SQL = "<script>"
+ "select st_x(t.geom) lon,st_y(t.geom) lat from biz_ceic_earthquake t where t.pk_id = #{id} "
+ "</script>";
@Select(FIND_ENTITY_SQL)
/**
* *找到地震震中的经纬度坐标，不直接用数据库中的字段，是因为数据库用的字符串来存的，不能直接计算，当然也可以使用函数进行转换，这样比较精确
* @param id 对应地震id
* @return
*/
CeicLonLatVO findEarthVoById(@Param("id")Long id);
static final String FIND_GEOJSON_SQL = "<script>"
+ "select st_asgeojson(ST_MakeEnvelope( #{lon} - (#{range} / (111.32 * cos(radians(#{lat})))),"
+ " #{lat} - (#{range} / 111.32),#{lon} + (#{range} / (111.32 * cos(radians(#{lat})))),#{lat} + (#{range} / 111.32),"
+ " st_srid(t1.geom) )) as geomJson from biz_ceic_earthquake t1 where t1.pk_id = #{id} "
+ "</script>";
@Select(FIND_GEOJSON_SQL)
/**
* *根据地震震中位置经纬度和指定的范围（单位是km）生成外包矩形框
* @param id
* @param lon 经度
* @param lat 纬度
* @param range 范围
* @return
*/
CeicEarthquake findRangGeoJsonById(@Param("id")Long id,@Param("lon") BigDecimal lon,@Param("lat") BigDecimal lat,@Param("range") BigDecimal range);
//这里注意在sql中的<>符号需要进行转义
static final String FIND_EARTH_LIST_SQL = "<script>"
+ "select t1.* from biz_ceic_earthquake t1 WHERE t1.m ::numeric > #{earthM} and t1.o_time >= #{time} "
+ " and t1.pk_id <> #{id} and st_contains( ST_MakeEnvelope( #{lon} - (#{range} / (111.32 * cos(radians(#{lat})))),"
+ " #{lat} - (#{range} / 111.32),#{lon} + (#{range} / (111.32 * cos(radians(#{lat})))),"
+ " #{lat} + (#{range} / 111.32),st_srid(t1.geom) ),t1.geom) "
+ "</script>";
@Select(FIND_EARTH_LIST_SQL)
/**
* *查询指定时间和指定范围内的地震信息列表
* @param id 排除当前的主键
* @param time 指定时间，如：2014-07-19
* @param lon 查询位置经度
* @param lat 查询位置纬度
* @param range 查询范围，单位公里
* @return
*/
List<CeicEarthquake> findEarthList(@Param("id")Long id,@Param("time") Date time,@Param("lon") BigDecimal lon,@Param("lat") BigDecimal lat,@Param("range") BigDecimal range,@Param("earthM") BigDecimal earthM);
}这里需要注意的一个地方就是，在最后一个Mapper方法中，我们有一个查询条件，即根据生成的空间范围查询时要把自己所在的点排除掉，这里采用where条件的<> 符号，如果直接在Script中写上述的脚本一定会报错的。需要进行转义。
<>2、业务层的设计与实现
由于我们想在页面打开的时候，一次性将4个地图的数据全部传递到前端，因此要求我们在业务逻辑层中将调用4次获取附近地震历史数据的的方法。数据我们使用hashMap来进行保存。具体的逻辑如下所示：
@Override
public Map<String, List<CeicEarthquake>> recentYearsMap(Long id) {
CeicLonLatVO vo = this.findEarthVoById(id);
//时间计算
Date fiveYears = DateUtils.modifyDate(DateUtils.getNowDate(), -1, 0, 0, 0, 0, 0);
Date tenYears = DateUtils.modifyDate(DateUtils.getNowDate(), -3, 0, 0, 0, 0, 0);
Map<String,List<CeicEarthquake>> result = new LinkedHashMap<String, List<CeicEarthquake>>(4);
BigDecimal range = new BigDecimal("100");
result.put("data1", this.baseMapper.findEarthList(id,fiveYears, vo.getLon(), vo.getLat(), range,new BigDecimal("3.0")));
result.put("data2", this.baseMapper.findEarthList(id,fiveYears, vo.getLon(), vo.getLat(), range,new BigDecimal("5.0")));
result.put("data3", this.baseMapper.findEarthList(id,tenYears, vo.getLon(), vo.getLat(), range,new BigDecimal("3.0")));
result.put("data4", this.baseMapper.findEarthList(id,tenYears, vo.getLon(), vo.getLat(), range,new BigDecimal("5.0")));
return result;
}在上述的方法中，主要是涉及一个时间的计算，其它的都是非常简单的。通过时间函数求去年的和近3年的时间点。
到此后台java的程序逻辑介绍完毕，控制层的逻辑比较简单不在赘述。
三、WebGIS功能设计与实现
在设计完后台程序之后，我们来进行WebGIS的功能设计与实现。将空间范围内进行展示，同时要把历史数据展示出来。
1、同时展示4幅地图
在这个场景需求中，我们需要同时展示4幅地图，用来表示不同的数据，与只展示一幅地图不同。4幅地图的初始化也不一样，篇幅有限。下面将正确的初始化方法进行介绍，防止大家在生成时走弯路。由于要展示的地图数量是4幅，因此我们采用数组的方式进行定义地图实例，同时还要展示不同地图的图例。
<div class="row">
<div class="col-sm-6">
<div id="mapid0" ></div>
</div>
<div class="col-sm-6">
<div id="mapid1" ></div>
</div>
</div>
<div class="row">
<div class="col-sm-6">
<div id="mapid2" ></div>
</div>
<div class="col-sm-6">
<div id="mapid3" ></div>
</div>
</div>var mapList = new Array();
var LegendArray = ["去年3级以上地震分布","去年5级以上地震分布","近3年3级以上地震分布 ","近3年5级以上地震分布 "];
$(function() {
initMap();
//展示震中矩形范围
showRange();
//展示近年的地震数据
showRecentYears();
//为了保证最好可视效果，重新设置视角
for(var i = 0;i< mapList.length;i++){
mapList[i].setView([lat, lon], 7);
}
});2、初始化地图
初始化地图的时候要注意，每一个地图都要新创建一个底图和标签图层的引用。否则地图的展示会有问题。
function initMap(){
for(var i = 0;i<4;i++){
$("#mapid" + i).height(wHeight / 2);//动态设置高度
var _canvasLabel = new L.CanvasLabel({
defaultLabelStyle: {
collisionFlg: true,
scale: 1,
strokeStyle: "white",
lineWidth:0,
font: "12px sans-serif"
}
});
var mapId = "mapid" + i;
var mymap = L.map(mapId,{renderer: _canvasLabel,crs:L.CRS.CustomEPSG4326,maxBounds:maxBounds,attributionControl:false});
mapList.push(mymap);
//底图
new L.tileLayer('http://localhost:8086/data/basemap_nowater/1_10_tms/{z}/{x}/{y}.jpg', {minZoom:1,
maxZoom: 16,
id: 'baseMap-nowater' + i,
tileSize: 256,
zoomOffset: -1
}).addTo(mymap);
//标签
new L.tileLayer('http://localhost:8086/data/basemap_nowater/1-10label/{z}/{x}/{y}.png', {maxZoom: 10,minZoom:1,
id: 'mapbox/label' + i ,tileSize: 256,zoomOffset: -1
}).addTo(mymap);
L.circleMarker(new L.LatLng(lat, lon), {radius: 8,
labelStyle: {
offsetX: 0, //横坐标偏移(像素)
offsetY: 15, //纵坐标偏移(像素)
text: cityInfo,
rotation: 0,
zIndex: 6,
minZoom : 2,
fillStyle: "#c50808"
},
color : "#c50808"
}).addTo(mymap);
L.control.Legend({
position: "bottomright",
collapsed: false,
symbolWidth: 12,
opacity: 1,
title:"图例",
column: 1,
legends: [ {
label: LegendArray[i],
type: "circle",
radius: 12,
color: "#c50808",
fillColor: "#c50808",
fillOpacity: 0.6,
weight: 0.5
}]
}).addTo(mymap);
mymap.setView([lat, lon], 8);
}
}3、展示100*100公里范围边界
这里我们首先在后台生成100*100公里的范围边界，并且以geojson的格式返回给前台。在Leaflet中可以直接加载GeoJSON到地图对象中。
function showRange(){
$.ajax({
type:"post",
url:prefix + "/historylist/range/" + [[${pkId + ''}]],
dataType:"json",
cache:false,
processData:false,
success:function(result){
if(result.code == web_status.SUCCESS){
for(var i = 0;i< mapList.length;i++){
L.geoJSON(JSON.parse(result.data),{style:{color:"red",weight:5,"opacity":0.6}}).addTo(mapList[i])
}
}
},
error:function(){
$.modal.alertWarning("获取信息失败");
}
});
}4、展示震中附近历年地震信息
这里调用后台的震中附件历年地震信息接口，在接口中将返回4个集合，我们只需要将后端返回的数据按照预定的方式使用Leaflet来进行加载即可。
function showRecentYears(){
$.ajax({
type:"post",
url:prefix + "/historylist/recentyears/" + [[${pkId + ''}]],
dataType:"json",
cache:false,
processData:false,
success:function(result){
if(result.code == web_status.SUCCESS){
var strokeStyleSet = "green";
recentyearsDetail(result.data.data1,mapList[0],"green");
recentyearsDetail(result.data.data2,mapList[1],"red");
recentyearsDetail(result.data.data3,mapList[2],"green");
recentyearsDetail(result.data.data4,mapList[3],"red");
}
},
error:function(){
$.modal.alertWarning("获取信息失败");
}
});
}
function recentyearsDetail(dataList,mapObj,showColor){
var strokeStyleSet = showColor;
for(var i=0;i<dataList.length;i++){
var dataInfo = dataList[i];
var marker = L.circleMarker(new L.LatLng(dataInfo.epiLat, dataInfo.epiLon), {radius: 8,labelStyle: {
offsetX: 0, //横坐标偏移(像素)
offsetY: 15, //纵坐标偏移(像素)
text: dataInfo.m + "/" + dataInfo.otime.substring(0,10),
rotation: 0,
zIndex: 1,
fillStyle:"white",
strokeStyle :strokeStyleSet
},
color : strokeStyleSet
}).addTo(mapObj);
}
}编写好上述代码，基本完成前台可视化界面的展示，下面我们结合实例来看一下实际的分析结果。
四、可视化成果展示
经过了数据库的设计与开发、Java后台程序的设计与实现以及WebGIS前端界面的设计与实现，我们基本实现了基于矩形范围面的时空空间分析。下面我们来看一下实际的效果。
1、新疆喀什叶城县
从可视化分析结果来看，叶城县应该是地震活动比较活跃的城市，虽然去年没有发生5级以上的地震，但是近三年3.0级以以上地震的次数非常多，同时可以看到仅3年内5.0级以上的破坏性地形就有好几次，最近一次是2022年的10月23日。
2、太平洋板块的花莲附近
看完喜马拉雅山脉附件的地震信息，我们再来看一下太平洋板块的，我国台湾省的花莲附件海域的地震信息，这也是一个位于明显的地震带的区域。首先来看下整体情况图：
从可视化分析结果来看，花莲县附近应该是地震活动更加活跃，不管是哪个地图都是密密麻麻的点，尤其时候右边的两个5.0级以上的地震，每一年都有发生。而且震级都不小，我们可以来看下去去年的地震信息。
通过上图可以看到，5.0级的地震相当频繁，甚至还有5.9级的大地震。是不是有点瑟瑟发抖的感觉。 而且密密麻麻的点是去年到现在一年以内的，真是伤心太平洋。
总结
以上就是本文的主要内容，文章主要讲解如何使用Java开发语言，结合成熟的SpringBoot框架来加速WebGIS应用的开发，首先讲解在数据库中的相关设计以及实现思路，然后重点讲解如何使用Java来开发地震震中附近历史地震的应用。通过本文的相关知识介绍，相信您一定会掌握如何进行附近历史地震分析有一个深刻的认识，同时抛砖引玉，对其他的类似的WebGIS应用建设提供建设参考。
