[发明专利]一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法

权利要求书

1.一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，包括以下步骤：

步骤1、加载二维模型参数；

步骤2、将探地雷达数值模拟区域分为模型区域和扩展区域；模型区域即二维模型所处区域；扩展区域是根据边界物性参数对模型区域进行扩展得到，包括过渡区域和PML区域，其中过渡区域为设置于PML区域和模型区域之间的一个有限宽度的区域；采用非结构化三角形单元对模型区域进行剖分，采用规则网格对扩展区域进行剖分，得到多个离散网格节点；

步骤3、自左向右设置M个测量点，即M个探地雷达接收天线位置，以及与M个探地雷达接收天线位置一一对应的M个探地雷达发射天线位置；初始化测量点序号m＝1；M为大于或等于1的整数；M个测量点分别设置于M个离散网格节点处；

步骤4、在第m个探地雷达发射天线位置加载激励源，计算与其对应的第m个测量点处各个时间步的电场值；

4.1、设置n为时间步序号，n的取值范围为[1，N]，N为最大时间步；设置第1个时间步及第1个时间步之前所有时间的电场与辅助场分量均为零；初始化时间步序号n＝1；

4.2、根据FDTD辅助场时间步进公式，使用第n个时间步的电场和第n-1/2个时间步的辅助场，更新第n+1/2个时间步的辅助场；其中，第n个时间步的电场，是由探地雷达数值模拟区域中所有离散网格节点第n个时间步的电场z方向分量组成的列向量；第n-1/2个时间步的辅助场，是由探地雷达数值模拟区域中所有离散网格节点第n-1/2个时间步的辅助场v方向分量组成的列向量，v＝x,y；其中，FDTD辅助场时间步进公式为：

其中，对于探地雷达数值模拟区域中任一离散网格节点，和表示其第n+1/2个时间步的辅助场v方向分量，v＝x,y，6个变量均为辅助变量，表示其第n个时间步的电场z方向分量；为差分系数，其中△t为时间步长；μ为磁导率，单位为H/m，σ_v为PML区域内v方向电导率参数，用于控制波的衰减，为大于零的正实数；α_v值的引入则是为了改善PML对低频分量的吸收特性，κ_v值的引入是为了改善PML对倏逝波的吸收特性，在PML区域内，α_v为大于零的正实数，κ_v为大于或等于1的正实数；在PML区域之外，α_v＝0，κ_v＝1；κ′_v，σ′_v，a′_v分别是κ_v，σ_v，a_v关于变量v的导数；

4.3、根据FETD主场时间步进公式，使用第n-1及第n个时间步的电场和第n+1/2个时间步的辅助场更新第n+1个时间步的电场；其中，FETD主场时间步进公式为：

其中，Eⁿ表示第n个时间步的电场，为第n+1/2个时间步的辅助场；f为激励源源向量，fⁿ表示第n个时间步激励源的值，K＝K_x+K_y，T＝T_x+T_y；t表示时间，采用Newmark-β法进行时间离散，t＝n·△t，△t为时间步长；β为Newmark-β法的参数；

其中，Ω代表整个积分区域，即整个探地雷达数值模拟区域；σ为电导率，单位为S/m，ε为介电常数，单位为F/m；N为形函数，上标T为转置符号；

4.4、令n＝n+1，返回步骤4.2，直到达到最大时间步为止；

4.5、分别从各个时间步的电场中提取出第m个测量点处各个时间步的电场值；

步骤5、判断当前测量点是否为最后一个测量点，即是否有m＝M，若是则转步骤6，否则令m＝m+1，返回步骤4；

步骤6、根据各个测量点处所有时间步的电场值得到探地雷达正演剖面图。

2.根据权利要求1所述的基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，所述步骤2中，所述模型区域采用非结构化Delaunay三角形网格进行剖分。

3.根据权利要求2所述的基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，所述步骤2中，在模型区域的剖分过程中，设置h_localv_loacl/(10f_max)，其中v_local表示局部速度，h_loal表示三角形单元的边长尺度，f_max表示子波的最大频率。