本发明涉及地震预测领域，尤其涉及一种地震预测规则的优化方法和装置。
背景技术：
地震预报是指在地震发生前，对未来地震发生的震级、时间和地点进行预测预报,向社会公布，及时开展灾害预防工作,减少地震灾害损失。目前的地震预测工作可以分为:长期预报、中期预报、短期预报、临震预报和主震后余震预报。由于地球的不可入性、破坏性地震的非频发性和地震发生过程的复杂性,以目前的科学水平，无法精确地震预报，地震精准的短短期临震预报始终是世界性科学难题。
1966年邢台7.2级地震后，我国开始地震预测及相应的研究工作，开展了一些列地震预测实用化攻关，逐步形成“长、中、短、临”和震后各阶段的渐进式预测科学思路(丁国瑜等，1981；梅世蓉等，1993)，形成了以预报经验数学化的统计预测方法，如综合概率法、模式识别、模糊聚类、专家系统以及模糊综合判定等(梅世蓉等，1993)。随着统计学理论、非线性理论、震源理论、计算机科学等领域的进一步发展，大规模数值计算的普及，使得这些方法和理论呈现出相互融合、多学科交叉的特点，衍生出了许多新方法及应用。
目前的地震预报正在经历经验预报向物理预报转变的过程，而地震预测的实际工作在很大程度上还是依赖经验预报方法。指标体系建设是经验预报阶段最有效的预测方法之一，因此，每个学科，每个地区都在建立地震预报指标体系(付虹等，2019)。为深化会商制度改革，不断完善、提高震例总结的水平和实效，中国地震局监测预报司对流体、测震、电磁、形变四大学科开展了地震预测指标体系清理工作，目前四大学科都在建立各自的指标体系。
地震预测指标可以定义为，根据观测资料经过数据处理提取出来的有一定物理意义的参数，在地震发生前，预测指标会发生某种程度的异常变化。按学科分类，测震学科的预测指标包括地震空区、地震条带、前兆震群、平静、活跃、各项地震学参数(b值、q值、d值、c值)、波速比、视应力等等(焦明若等，2016)，电磁学科的预测指标包括地磁低点位移、地磁加卸载响应比、地磁逐日比、地磁垂直分量日变化等等(冯志生等，2019)，形变学科的预测指标包括定点观测和大地测量数据提取的参数，流体学科的预测指标包括水温、水位和各项地球化学测项。
一般预测指标的异常变化和未来发生的地震具有时空关联性，通过震例回溯和统计检验，根据经验确定预测指标的预测规则。选定了预测指标和预测规则，就课余确定预测预测的时空范围，就可以在实际工作中应用。
地震预报方法的有效性，不能以一次预报的成功与否来评价，需要应用统计学理论来进行检验和评分(朱令人等，2004)；显然，地震预测指标的有效性同样需要用科学的统计方法进行检验和评价。概率预测和统计检验是地震中长期地震危险性分析不可或缺的组成部分，也是设计地震预测“警报”阈值、优化预测策略和预测效能评价的有效手段(蒋长胜等，2011)。
预测指标可分为时间相关的指标和时空相关的指标，时间相关的指标仅判断时间域上地震是否发生，时空相关预测指标在三维的时空域(时间、维度、经度)上判断地震是否发生。
目前的预测规则，主要靠人工统计，具有一定主观性。
技术实现要素：
本发明的实施例提供了一种地震预测规则的优化方法和装置，能够优化地震的预测规则，提高地震预报准确度。
一种地震预测规则的优化方法，包括：
步骤1，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
步骤2，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
步骤3，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
步骤4，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
步骤5，在所述预测规则的参数空间中搜索，按照所述步骤4的方法计算预定次数的地震评价指标值；
步骤6，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
一种地震预测规则的优化装置，包括：
第一设定单元，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
第二设定单元，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
获取单元，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
第一计算单元，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
第二计算单元，在所述预测规则的参数空间中搜索，按照所述步骤4的方法计算预定次数的地震评价指标值；
选择单元，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
由上述本发明的实施例提供的技术方案可以看出，本发明实施例中，设计时空相关预测指标效能检验的整体方案和高效算法，实现了时空相关预测指标的效能评估和预测规则的优化，从而提高地震预报准确度。
本发明附加的方面和优点将在下面的描述中部分给出，这些将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明所述的一种地震预测规则的优化方法的示意图；
图2为本发明所述的一种地震预测规则的优化方法的应用场景的示意图。
图3为本发明所述的三维坐标系和关联的圆柱形区域的时空区域图示。
具体实施方式
下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
为便于对本发明实施例的理解，下面将结合附图以几个具体实施例为例做进一步的解释说明，且各个实施例并不构成对本发明实施例的限定。
如图1所示，为本发明所述的一种地震预测规则的优化方法，包括：
步骤11，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
步骤12，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
步骤13，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
步骤14，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
步骤15，在所述预测规则的参数空间中搜索，按照所述步骤14的方法计算预定次数的地震评价指标值；
步骤16，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
其中，当所述评价指标值为r值时，所述计算地震评价指标值的步骤具体为：
r值＝击中率-时空占有率。
其中，当所述评价指标值为概率增益值时，所述计算地震评价指标值的步骤具体为：概率增益＝击中率/时空占有率。
其中，所述击中率的获取步骤包括：
击中率＝预测成功地震数/目标地震总数；
其中，所述预测成功地震数的获取步骤包括：
获取所述目标地震目录中每个地震所在的时空网格点，对所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域；若判断为异常区域，表明该地震预测成功，以此统计预测成功地震数。
其中，所述时空占有率的获取步骤包括：
对所述研究区域按照所述网格化步长进行网格化后，使用蒙特卡洛方法，随机选取预定数量的时空网格点；
对选取的所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域；
将判定为异常区域的时空网格点数占抽样网格总数的比例作为时空占有率。
其中，所述所述对所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域的步骤为：
根据预测规则中的空间预测半径和预测时间窗长，确定待判断网格点有关的时空范围是时空域的一个圆柱形区域；从所述指标扫描结果筛选出所述圆柱形区域包含的网格点；判断所述网格点的预测指标值是否符合报警条件；当在筛选出的所述圆柱形区域的对应的三维坐标系中有一个格点的预测指标符合报警条件，则表明待判断的格点属于异常区域。例如：要判断某个点(设为点a)是否属于异常区域，看这个点(点a)属于哪些点(设为点b1、…、bn)的预测范围。根据预测规则中的预测时间窗长和预测半径，这些点(点b1、…、bn)处于时空范围(三维空间，时间1维度，空间2维)一个圆柱形区域，再看这些点的预测指标值(根据点b1、…、bn处点预测指标值判断)有没有达到报警条件的，只要这些点有一个达到报警条件，就说明要判断的点属于异常区域。设要判断的点为a，关联区域是指预测区域包含要判断点a的所有点的集合。
本发明还提供一种地震预测规则的优化装置，包括：
第一设定单元，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
第二设定单元，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
获取单元，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
第一计算单元，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
第二计算单元，在所述预测规则的参数空间中搜索，按照所述步骤4的方法计算预定次数的地震评价指标值；
选择单元，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
以下描述本发明的应用场景。本发明主要解决时空相关预测指标的评估问题，设计时空相关预测指标效能检验的整体方案和高效算法，实现了时空相关预测指标的效能评估和预测规则的优化。本发明用于判断网格点是否异常区域的逆向搜索算法，降低计算复杂度，提高计算效率。本发明用于统计时空占有率的蒙特卡洛抽样算法，将研究时空范围网格化，随机选取一定数量的网格点，判断每个网格点是否属于异常区域，将判定为异常区域的网格点数占抽样网格总数的比例作为时空占有率的近似值。
如图2所示，工作流程说明：
(1)设定研究区域
研究区域为三维的时空域(空间范围两个维度、时间一个维度)，设定网格化步长，包括时间的步长和空间的步长
(2)预测规则设定
预测规则包括空间半径、时间窗长、报警条件(指标的阈值)、预测震级
(3)数据准备
研究区域的预测指标扫描结果(经过网格离散化)，目标地震目录
(4)异常区域统计和时空占有率计算
对研究区域网格化后，使用蒙特卡洛方法，随机选取一定数量网格点
对抽样的网格点，使用逆向搜索算法，判定每个网格点是否属于异常区域
(5)击中率计算
对目标地震目录中每个地震所在的时空网格点，使用逆向搜索算法，判定是否属于异常区域，若判断为异常区域，表明该地震和预测成功。
击中率＝预测成功地震数/目标地震总数
(6)概率增益计算
概率增益＝击中率/时空占有率
(7)r值计算
r值＝击中率-时空占有率
(8)判断网格点是否异常区域的逆向搜索算法
根据预测规则中的预测半径和预测时间窗长，确定和要判断网格点有关的时空范围，是时空域的一个圆柱形区域，筛选出这个区域包含的网格点，判断这些网格点的预测指标值是否符合报警条件，只要在这个关联的区域中有一个格点的预测指标符合报警条件，就说明要判断的格点属于异常区域。如图3所示，为本发明所述的三维坐标系和关联的圆柱形区域的时空区域图示。
说明：将预测的时空区域网格化，根据每个格点的指标值，根据预测规则，确定其对应的预测区域，并求出超过阈值的格点对应异常区域的并集，这个计算过程计算量较大。我们使用的等价的算法(逆向搜索算法)，直接判断每个格点是否属于异常区域，具体方法是，根据预测规则，找出和这个格点关联的格点，如果这些格点中存在符合报警条件的点，就说明这个格点属于异常区域，这样就大大减少了计算量，只需要求一次最值。
(9)统计时空占有率的蒙特卡洛抽样算法
将研究时空范围网格化，随机选取一定数量的网格点，(使用逆向搜索算法)判断每个网格点是否属于异常区域，将判定为异常区域的网格点数占抽样网格总数的比例作为时空占有率的近似值
(10)预测规则的优化
在预测规则的参数空间中搜索，多次计算r值和概率增益，选取r值或者概率增益最高时对应的预测规则，作为最优的预测规则，达到预测规则优化的目的。
本发明中，具有以下有益效果：
1、本发明设计时空相关预测指标效能检验的整体方案和高效算法，实现了时空相关预测指标的效能评估和预测规则的优化。
2、本发明用于判断网格点是否异常区域的逆向搜索算法，降低计算复杂度，提高计算效率；
3、本发明用于统计时空占有率的蒙特卡洛抽样算法，将研究时空范围网格化，随机选取一定数量的网格点，判断每个网格点是否属于异常区域，将判定为异常区域的网格点数占抽样网格总数的比例作为时空占有率的近似值。此算法的目的是避免计算全部网格点，计算计算量。
4、判断网格点是否异常区域逆向搜索算法：在指标和预测规则确定后，计算异常区域是一个求大量时空范围并集的运算，计算量大且过程复杂，在发明设计了逆向搜索的算法，将求并集的运算简化为求最值的运算。
5、计算时空占有率的蒙特卡罗算法：
在网格集合中随机抽样，判断每个是否属于异常区域，计算属于异常区域的网格点占抽样总数的比例，得到时空占有率的近似值。此算法的目的是避免计算全部网格点，计算计算量。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。
1.一种地震预测规则的优化方法，其特征在于，包括：
步骤1，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
步骤2，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
步骤3，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
步骤4，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
步骤5，在所述预测规则的参数空间中搜索，按照所述步骤4的方法计算预定次数的地震评价指标值；
步骤6，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
2.根据权利要求1所述的方法，其特征在于，当所述评价指标值为r值时，
所述计算地震评价指标值的步骤具体为：
r值=击中率-时空占有率。
3.根据权利要求1所述的方法，其特征在于，当所述评价指标值为概率增益值时，所述计算地震评价指标值的步骤具体为：
概率增益=击中率/时空占有率。
4.根据权利要求2或3所述的方法，其特征在于，
所述击中率的获取步骤包括：
击中率=预测成功地震数/目标地震总数；
所述预测成功地震数的获取步骤包括：
获取所述目标地震目录中每个地震所在的时空网格点，对所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域；若判断为异常区域，表明该地震预测成功，以此统计预测成功地震数。
5.根据权利要求2或3所述的方法，其特征在于，
所述时空占有率的获取步骤包括：
对所述研究区域按照所述网格化步长进行网格化后，使用蒙特卡洛方法，随机选取预定数量的时空网格点；
对选取的所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域；
将判定为异常区域的时空网格点数占抽样网格总数的比例作为时空占有率。
6.根据权利要求4或5所述的方法，其特征在于，所述对所述时空网格点使用逆向搜索算法，判定每个时空网格点是否属于异常区域的步骤为：
根据预测规则中的空间预测半径和预测时间窗长，确定待判断网格点有关的时空范围是时空域的一个圆柱形区域；从所述指标扫描结果筛选出所述圆柱形区域包含的网格点；判断所述网格点的预测指标值是否符合报警条件；当在筛选出的所述圆柱形区域的对应的三维坐标系中有一个格点的预测指标符合报警条件，则表明待判断的格点属于异常区域。
7.一种地震预测规则的优化装置，其特征在于，包括：
第一设定单元，设定地震的研究区域，所述研究区域由：空间范围经度、纬度、预测时间窗长确定，通过三维坐标系表示；设定网格化步长，所述网格化步长包括：时间步长和空间步长；
第二设定单元，设定地震的预测规则，所述预测规则包括：空间预测半径、预测时间窗长、报警条件、预测震级；
获取单元，获取所述研究区域的预测指标扫描结果和目标地震目录；所述预测指标扫描结果为：按照预定时间步长和空间步长，对时空网络每个点扫描计算得到的数据，数据格式为由时间、空间确定的三维数组；所述预测指标是由观测数据提取的有物理意义的参数；所述目标地震目录为：研究区域内预定震级以上的全部地震数据；
第一计算单元，根据所述预测指标扫描结果和目标地震目录，计算地震评价指标值；
第二计算单元，在所述预测规则的参数空间中搜索，按照所述步骤4的方法计算预定次数的地震评价指标值；
选择单元，选取评价指标值最高时对应的预测规则的参数，作为最优的预测规则，并输出所述最优的预测规则。
技术总结