一种用于波形互相关双差定位算法的并行优化方法与流程

本发明属于并行优化领域，涉及一种用于波形互相关双差定位算法的并行优化方法。

背景技术：

双差定位法是一种受限制比较小的相对定位算法,其基本思想是:当两个微地震震源之间的距离相对于它们各自到观测点间的距离较小，且传播路径上的速度变化也较小时，这两个震源到同一个接收点的射线路径基本相同或者相近，这使得他们之间的走时残差与两震源点之间的距离有较强的相关性，从而降低了传播路径上速度对定位结果的影响。

波形互相关利用了“发震机理相同地震波形相似”的原理,通过对相似地震波形进行互相关分析,计算同一震相的走时差,主要有时间域和频率域两种方式。通过时间域(频率域)相关函数计算相似波形的相关系数,设定相关系数阈值,获取可用的到时差数据。

结合波形互相关的双差定位法被广泛用于大范围空间尺度上大量地震事件的重新定位中。davidp.schaff等通过计算震相的互相关系数获取震相到时差,并采用双差算法1999年发生在我国辽宁岫岩地震序列进行了重新定位,结果显示使用分布稀疏的区域台站(震中距在500～1000km)仍然可将定位精度提高到5km以内。杨智娴在国内首次将相对定位算法在我国中西部大范围内使用,经重新精确定位后显示出我国中西部地区地震震源分布更加精细的图像。朱艾斓使用双差地震定位法对川西地区1992—2002年的13367个小震进行重新定位,重新定位后,地震活动沿活动断裂成线状分布现象非常突出,呈现出与地表活动构造的密切关系。黄媛将近台记录和区域台网数据联合用于双差地震定位算法,对2003年新疆巴楚-伽师ms6.8强震后404个ml≥3.5余震序列进行双差法重新精确定位,并对其进行了理论上的可行性分析,结果获得初步成功。赵翠萍对1997—2002年新疆伽师震源区150多个地震,共计几百条波形数据进行了互相关分析,虽然由于可用台站的波形数据不多,得到的互相关走时结果不多,对定位精度的改进作用也不是很明显,但是通过这项工作,实现了批量计算波形互相关走时的可能性。

技术实现要素：

为了解决大规模应用波形互相关双差定位算法时的效率问题，本发明旨在提出一种用于波形互相关双差定位算法的并行优化方法，通过对核心步骤加速并行化，从而提高计算效率。

一种用于波形互相关双差定位算法的并行优化方法，步骤如下：

步骤1，获取总任务数和总进程数，并将总任务平均分给每个进程，并得出每个进程的任务初始位置。

步骤2，计算获取初始匹配位置as1bs2。

步骤3，计算获取结束匹配位置ae1be2。

步骤4，根据获取的as1bs2位置和ae1be2位置读取对应数据进行计算。

作为本发明的进一步改进，步骤1具体包括：

步骤101，首先选定两个台站，台站a和台站b，每个台站都有n天的数据，a台站第i天的数据可表示为ai，b台站第i天的数据可表示为bi，台站a和台站b之间的数据需要进行匹配计算，匹配方式是对每个ai，都需要和bj(j>i)进行匹配计算，匹配计算格式定义为aibj(1≤i<j≤n)，由此便可得到a和b匹配计算总的任务数为个，定义为pairs_len。

步骤102，将总进程个数定义为np，目标是求总任务平均分配给每个进程后，每个进程的任务初始匹配位置as1bs2和结束匹配位置ae1be2。

步骤103，将pairs_len整除np，得到每个进程要处理的任务数量，定义为task_numk，k表示进程索引，再将pairs_len余除np，得到暂未被分配到进程的任务数量，定义为residul_day；将residul_day个任务数分配到总进程中前residul_day个进程，所以前residul_day个进程任务数task_numk(k≤residul_day)自加1；至此完成进程任务数量分配。

步骤104，根据每个进程需要处理的任务数量task_numk，累加计算可以得到每个进程的任务初始位置task_num_startp，p表示进程索引，k＝1,2,3...p...。

作为本发明的进一步改进，步骤2具体包括：

步骤201，将ai需要做匹配计算的个数定义为a_leni，定义a_sumj为第j天之前需要做匹配计算数量。

步骤202，遍历整个a_sum数组，找到第一个大于每个进程中task_num_startp值的a_sumt，并保存t值。

步骤203，t-1即为开始匹配位置as1bs2中的s1，再由t-1值得到a_sumt-1。

步骤204，将task_num_startp-a_sumt-1即可得到每个进程在第t天的需要处理的任务数量

步骤205，再把步骤204中的task_num_startp-a_sumt-1值加上t，即可得到每个进程的初始匹配位置as1bs2中的s2，至此得到了每个进程的初始匹配位置as1bs2。

作为本发明的进一步改进，步骤3具体包括：

步骤301，定义变量nday1用于记录每个进程处理的任务aibj中imax-imin+1的值，nday1代表ai中i的个数。

步骤302，定义变量nday2用于记录每个进程处理的任务aibj中jmax-jmin+1的值，nday2代表bj中j的个数。

步骤303，定义变量day1_min，day2_min用于记录每个进程处理的任务aibj中最小的i和j；同时定义day1_start，day2_start表示步骤203、步骤205求出的as1bs2中的s1和s2。

步骤304，将day1_min+nday1-1即可得到得到ae1be2中的e1，将day2_start-day2_min+nday1即可得到ae1be2中的e2。

作为本发明的进一步改进，步骤4具体包括：

步骤401，定义结构体数组stk，k为进程索引，根据每个进程的初始匹配位置as1bs2和结束匹配位置ae1be2读入数据到结构体数组stk中。

步骤402，对结构体数组stk进行计算，公式为

t为ai、bj中数据的索引，至此得计算结果。

与常规计算相比，本发明通过大规模并行分配任务至每个进程，极大的提高了计算效率。当计算数据规模较大时，或者计算复杂度高时，仅用普通常规计算难以得到较好的效率，这时便需要设计合理高效的并行优化方案。本发明提供的技术利用已经公开的资料作为基础，构建实用的并行加速方案，大量减少了数据冗余问题，极大提高了效率。

附图说明

图1为本发明的并行优化方法的流程示意图。

图2为图1中步骤1的流程示意图。

图3为图1中步骤2的流程示意图。

图4为图1中步骤3的流程示意图。

图5为图1中步骤4的流程示意图。

图6为串行和并行模式下程序的加速比和计算时间比较图。

具体实施方式

下面通过具体的实施例并结合附图对本发明做进一步的详细描述。

实施例1

如图1所示，本发明实施例所述的一种用于波形互相关双差定位算法的并行优化方法，该方法包括：

步骤1，总任务数个和总进程数1024个，然后将总任务数平均分给每个进程，得出每个进程的任务初始位置。

步骤101，首先选定两个台站，台站a和台站b，每个台站都有731天的数据，a台站第i天的数据可表示为ai，b台站第i天的数据可表示为bi，台站a和台站b之间的数据需要进行匹配计算，匹配方式是对每个ai，都需要和bj(j>i)进行匹配计算，匹配计算格式定义为aibj(1≤i<j≤n)，由此便可得到a和b匹配计算总的任务数为个。

步骤102，进程总共1024个，目标是求总任务平均分配给每个进程后，每个进程的任务初始匹配位置as1bs2，和结束匹配位置ae1be2。

步骤103，将266815整除1024，得到每个进程要处理的任务数量为260个，定义为task_numk，k表示进程索引，再将266815余除1024，得到暂未被分配到进程的任务数量为575个，将575个任务数分配到总进程中前575个进程，所以前575个进程任务数task_numk(k≤575)自加1；至此完成进程任务数量分配。

步骤104，根据每个进程需要处理的任务数量task_numk，累加计算可以得到每个进程的任务初始位置task_num_startp，p表示进程索引。

步骤2，计算获取初始匹配位置as1bs2。

步骤201，将ai需要做匹配计算的个数定义为a_leni，设a_sumj为第j天之前需要做匹配计算数量。

步骤202，遍历整个a_sum数组，找到第一个大于每个进程中task_num_startp值的a_sumt，并保存t值。

步骤203，t-1即为开始匹配位置as1bs2中的s1，再由t-1值得到a_sumt-1。

步骤204，将task_num_startp-a_sumt-1即可得到每个进程在第t天的需要处理的任务数量。

步骤205，再把步骤204中的task_num_startp-a_sumt-1值加上t，即可得到每个进程的初始匹配位置as1bs2中的s2，至此得到了每个进程的初始匹配位置as1bs2。

步骤3，计算获取结束匹配位置ae1be2。

步骤301，定义变量nday1用于记录每个进程处理的任务aibj中imax-imin+1的值，nday1代表ai中i的个数。

步骤302，定义变量nday2用于记录每个进程处理的任务aibj中jmax-jmin+1的值，nday2代表bj中j的个数。

步骤303，定义变量day1_min，day2_min用于记录每个进程处理的任务aibj中最小的i和j；同时定义day1_start，day2_start表示步骤203、步骤205求出的as1bs2中的s1和s2。

步骤304，将day1_min+nday1-1即可得到得到ae1be2中的e1，将day2_start-day2_min+nday1即可得到ae1be2中的e2。

步骤4，根据获取的as1bs2位置和ae1be2位置读取对应数据进行计算。

步骤401，定义结构体数组stk，k为进程索引，根据每个进程的初始匹配位置as1bs2和结束匹配位置ae1be2读入数据到结构体数组stk中。

步骤402，对结构体数组stk进行计算，公式为

t为ai、bj中数据的索引，至此得到计算结果。

图6为串行和并行模式下程序的加速比和计算时间比较，左轴为加速比，右轴为计算时间，单位小时。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。