CN103969627B - 基于fdtd的探地雷达大规模三维正演模拟方法 - Google Patents

Abstract

本发明公开了一种基于FDTD的探地雷达大规模三维正演模拟方法，它通过网络进行各节点数据的交换，实现各节点通讯和同步；选择任一节点中的一个测点基于FDTD正演并行计算，首先计算电场值，电场计算完成后加入激励源计算磁场，通过通讯、电场计算、磁场计算不停的迭代，计算到设定的时间总步长后停止，保存输出的数据，该数据就是这个测点的雷达反射信号；对所有场值进行初始化更新，使得整个计算域恢复到计算开始前的状态，然后移动发射天线和接收天线的位置到下一个测点重新开始正演计算；若所有测点计算完成，把每个节点的数据汇集以进行成像和图像后处理。本发明能够适用于各种尺度复杂地电介质和模型的计算，具有通用性强的特点，特别适合于大尺度、高精度的正演模拟。

Description

基于FDTD的探地雷达大规模三维正演模拟方法

技术领域

本发明涉及用于作为探地雷达理论研究手段的探地雷达正演模拟方法。

背景技术

探地雷达(GroundPenetratingRadar，GPR)方法是一种用于确定地下介质分布的光谱(1MHz-1GHz)电磁波技术。探地雷达是一种较新的地球物理方法，在近10年的时间内，逐渐成熟起来。探地雷达的发展既伴随着各种各样军用和民用的需求，又得到高新技术发展的推动

目前，有关探地雷达探测与解释技术大都建立在二维的基础上，二维雷达探测基于地质体在垂直探测剖面的方向是均匀不变的假设，即地质体是无限延伸的，但在工程实际中，常见的地质体如岩溶洞穴是不规则的三维体，因此，二维雷达探测通常能探测到地质体是否存在及并对其定位，但难于对探测对象提供更为准确的信息，如探测对象的形状、空间信息。因此，为提高探地雷达探测的准确度，开展三维探测与解释技术就显得特别重要。

探地雷达三维正演模拟是通过数值计算的方法来模拟雷达波在含有被测物体的地电介质中传播的全过程。它是分析探测问题、研究电磁波在介质中传播规律的有效手段，是三维探测与解释技术的基础。例如地下洞穴的探测、钢筋混泥土中缺陷的探测可归结为三维形体的探测问题，而近地土壤结果、岩石分层、地下水位面等具有三维曲面的特征。研究雷达波在这些问题的介质中传播，有利于提高探测的效果和解释的准确性。探地雷达波具有高频特征，波长较短，介质吸收强烈，加之受地面干扰大，使得探测剖面较为困难。因此，数值模拟复杂形体存在时的雷达波场特征，对认识实际的雷达记录，识别目标有着重要的意义。

探地雷达的正演模拟技术一直备受国内外很多学者的关注，尤其是近二十年来，随着计算机技术的飞速发展，各国学者纷纷把正演模拟作为探地雷达理论研究的主要手段。目前，正演模拟算法主要有射线追踪法、有限元法(FEM)以及时域有限差分法(FiniteDifferenceTimeDomain，FDTD)。自从1994年用于截断计算空间的完全匹配层(PerfectlyMarchedLayer,PML)技术被发现后，FDTD逐渐发展成为日趋成熟的计算电磁学方法，并由于其独特的优势，目前成为最主流的探地雷达正演算法。和其他算法相比，它具有以下几个优势:(1)时域计算，能设置超宽带激励源；(2)吸收边界条件优良；(3)算法复杂度较低；(4)能处理复杂的媒质。

尽管FDTD方法有许多优势，但是在进行大尺寸探地雷达的正演模拟时会面临诸多问题。首先，FDTD算法必须存储整个空间的电介质和电磁场信息，为了避免网格色散误差以及对缝隙等薄结构进行建模，空间步长必须足够小，通常取最小波长的1/20。当计算空间是三维且尺度远大于波长时(即电大尺度),需要的存储空间完全超出目前单机的极限，使得计算无法进行下去；其次，FDTD方法要获得稳定的数值解，时间步长必须满足，即Courant稳定条件，这导致三维情况下探地雷达单个测点所需的计算时间很长；第三，探地雷达进行三维探测时为了提高精度，往往需要众多的测线，正演模拟需要依次计算每一个测点完整的雷达波过程，当测点分布密集，数量庞大时，需要的计算时间是极其漫长的，甚至根本无法计算。因此，目前只能进行较小规模的模拟，对大尺度、大规模的正演模拟还无能为力。随着探地雷达技术的不断进步，探测精度、解释能力要求的不断提高，迫切需要进行大尺度的三维正演模拟。这对于探地雷达理论研究，探测方式的研究具有重要意义。

发明内容

本发明目的是为了克服现有技术的不足而提供一种适用于大规模三维正演模拟的基于FDTD探底雷达正演方法。

为达到上述目的，本发明采用的技术方案是：一种基于FDTD的探地雷达大规模三维正演模拟方法，它包括如下步骤，

(1)采用基于消息传递的并行运算模式的MPL并行库，建立并行计算节点拓扑，确定节点通讯组的大小，将目标空间区域离散化，形成多个子区域，每个子区域对应一个节点，且在子区域上根据精度要求设置多个测点；

(2)初始化MPL并行库，并根据正演需求设置计算域大小、网格精度、时间总步长、激励源参数、地电模型、地下介质模型、吸收边界参数等，根据这些参数计算需要的时间步长、空间步长等中间数据，获得当前节点的ID，并在各节点分配内存空间并初始化；

(3)通过网络进行各节点数据的交换，实现各节点通讯和同步；选择任一个测点基于FDTD正演对所有节点并行计算，首先计算电场值，电场计算完成后加入激励源计算磁场，通过通讯、电场计算、磁场计算不停的迭代，计算到设定的时间总步长后停止，保存输出的数据，该数据就是这个测点的雷达反射信号；

(4)步骤(3)完成后，对包括电场和磁场以及天线激励源等所有场值进行初始化更新，使得整个计算域恢复到计算开始前的状态，然后移动发射天线和接收天线的位置到下一个测点重新开始正演计算；

(5)若所有测点计算完成，把每个节点的数据汇集以进行成像和图像后处理；若测点未计算完成，重复进入步骤(3)(4)。

优化地，步骤(4)中，所述的FDTD正演并行计算中，地下物体物性参数的分配、任务分配和最终数据的收集是采用主从式方法，将所述的节点设为若干子节点和一个主节点，所述的主节点将任务分解，分配到各个子节点，各子节点自行计算，所有测点计算完成后，把每个子节点的数据汇集到主节点，在主节点进行成像和图像后处理。

优化地，所述的FDTD通过分治策略来进行并行计算，即通过区域分割，每个计算节点分配一个区域来实现大尺寸的计算，每个节点并行地完成各自的计算任务，包括计算过程中的各节点的数据交换(施行通信同步)，最后将结果汇总起来。

优化地，所述的区域分割采用三维分割的方式。

优化地，步骤(4)中，所述的数据交换采取交换磁场的方式，相邻节点的计算区域均包含一个重合网格，当一个节点进行FDTD电场迭代计算时，把相邻节点该位置上的磁场值传递给该节点作为其边界磁场。

优化地，在执行步骤(4)前，还需要在计算空间外围设置用于吸收传播到边界的电磁波的吸收边界。

优化地，所述的吸收边界设置方式为：按照指数分布的方式从内到外设置具有电导率、磁导率的吸收层，以使得电磁波逐层衰减。

优化地，步骤(4)中，所述的激励源的设置用超宽带Recker子波。

由于上述技术方案运用，本发明与现有技术相比具有下列优点：本发明把并行机计算节点进行分组，每一组由若干个节点组成，依靠着这些节点进行单独的大尺度三维FDTD并行计算。每组负责计算若干个测点，把这些组的计算数据组合起来，最后可以得到整个三维探地雷达剖面。这样不仅能充分保证节点之间负载的平衡，还能进行各组之间无耦合的并行计算，最大限度降低因为计算节点故障而需要重新计算造成的时间和昂贵机时费的损失。该方法能够适用于各种尺度复杂地电介质和模型的计算，具有通用性强的特点，特别适合于大尺度、高精度的正演模拟，有望提高探地雷达探测的可靠性、准确度及雷达剖面的分辨率，更有利于探地雷达资料的解释。

附图说明

图1为本发明基于FDTD探地雷达三维正演模拟计算流程图；

图2为FDTD的电磁场空间离散策略示意图；

图3为本发明主从式并行方法原理图；

图4为FDTD三种区间分割方式图；

图5为基于FDTD各节点并行方法原理图；

图6为发明相邻节点间通讯交换数据的示意图；

图7为本发明采用的激励源波形图；

图8为图7波形的频谱图；

图9为本发明并行程序用不同节点数计算的加速比测试图。

具体实施方式

下面将结合附图对本发明优选实施方案进行详细说明：

图1给出了本发明整个三维正演并行计算流程图。首先将目标空间离散化，形成多个节点，每个节点上根据精度要求设置多个测点；初始化MPI并行环境，建立节点拓扑，确定通讯组的大小，获得当前节点的ID。对正演计算各种物理量进行设置、分配内存空间、初始化。然后根据正演需求设置计算域大小、网格精度、时间总步长、天线参数、地电模型、地下介质模型等，根据这些参数计算需要的时间步长、空间步长等中间数据。前期并行框架、正演计算域框架搭建完毕后，开始第一个测点的正演计算。在迭代的开始，通过网络进行各节点数据的交换(边缘磁场的交换)，完成通讯和同步。接着通过FDTD算法计算电场值。电场计算完成后加入激励源，然后通过FDTD算法计算磁场。通过通讯、电场计算、磁场计算不停的迭代，计算到时间总步长以后停止。保存输出的数据，该数据就是这个测点的雷达反射信号。这个测点完成后，对所有场值进行初始化更新，使得整个计算域恢复到计算开始前的状态，然后移动发射天线和接收天线的位置到下一个测点重新开始正演计算。这样对所有的测点进行扫描，最终获得每个测点的雷达反射数据。所有测点计算完成后，把每个节点的数据汇集，进行成像和图像后处理。

再对上述流程进一步说明之前，我们简要介绍下探地雷达GRP正演计算基本方程：

GPR所研究的问题一般只涉及线性各向同性且与时间无关的媒质，含电和磁衰减的GPR电磁波传播基本规律可用经典的Maxwell方程来描述：

$\begin{array}{c}\▿\×E=-\mathrm{\μ}\frac{\∂H}{\∂t}-J\\ \▿\×H=-\mathrm{\ϵ}\frac{\∂E}{\∂t}-J_m\\ J=\mathrm{\σE}\\ J_m={\mathrm{\σ}}_{m}H\end{array}---\left(1\right)$

其中：E为电场强度，V/m；H为磁场强度，A/m；ε为介电常数，F/m；μ为磁导率，H/m；J是电流密度，A/m2；Jm为磁流密度，V/m2；σ为电导率，S/m；σm为等效磁阻率，Ω/m。在导出差分方程之前，需要写出与方程(1)等价的电磁场各分量分别满足的方程，在直角坐标系中，可得：

$\begin{array}{c}\frac{\∂H_s}{\∂y}-\frac{\∂H_y}{\∂c}=\mathrm{\ϵ}\frac{\∂E_x}{\∂t}+\mathrm{\σ}{E}_x\\ \frac{\∂H_x}{\∂c}-\frac{\∂H_s}{\∂x}=\mathrm{\ϵ}\frac{\∂E_y}{\∂t}+\mathrm{\σ}{E}_y\\ \frac{\∂H_y}{\∂x}-\frac{\∂H_x}{\∂y}=\mathrm{\ϵ}\frac{\∂E_s}{\∂t}+\mathrm{\σ}{E}_s\\ \frac{\∂E_s}{\∂y}-\frac{\∂E_y}{\∂c}=-\mathrm{\μ}\frac{\∂H_x}{\∂t}-{\mathrm{\δ}}_m{H}_x\\ \frac{\∂E_x}{\∂c}-\frac{\∂E_s}{\∂x}=-\mathrm{\μ}\frac{\∂H_y}{\∂t}-{\mathrm{\σ}}_m{H}_y\\ \frac{\∂E_y}{\∂x}-\frac{\∂E_x}{\∂y}=-\mathrm{\μ}\frac{\∂H_s}{\∂t}-{\mathrm{\σ}}_m{H}_s\end{array}----\left(2\right)$

为了建立差分方程，首先要在变量空间将连续变量离散化。通常用某种形式的网格划分变量空间，且只取网格点上的未知量作为计算对象。1966年K.S.Yee首次提出了电磁场数值计算的一种新方法，即FDTD，电磁场E、H分量在空间和时间上采取交替抽样的离散方式，每一个E(或H)场分量周围有四个H(或E)分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场，使计算过程完全模拟电磁波逐步传播以及与目标物体相互作用的过程。图2给出了FDTD的电磁场空间离散策略。

时间和空间的离散取以下符号表示：

f(x,y,z,t)＝f(i△x,j△y,k△z,n△t)＝fⁿ(i,j,k)(3)

其中△x，△y，△z，△t分别是x、y、z方向的网格步长以及时间步长。则以方程组(2)中Ez计算方程为例，为了提高计算精度电流密度的表达式用时间中心差分，其余用Yee氏网格差分格式，可得：

${\mathrm{Ez}}_{i,j,k}^{n+1}={\mathrm{Ez}}_{i,j,k}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ}}(\frac{{\mathrm{Hy}}_{i,j+1/2,k}^{n+1/2}-{\mathrm{Hy}}_{i,j-1/2,k}^{n+1/2}}{\mathrm{\Δx}}-\frac{{\mathrm{Hx}}_{i,j,k+1/2}^{n+1/2}-{\mathrm{Hx}}_{i,j,k-1/2}^{n+1/2}}{\mathrm{\Δy}})-\frac{\mathrm{\σ\Δt}}{2}({\mathrm{Ez}}_{i,j,k}^n+{\mathrm{Ez}}_{i,j,k}^{n+1})$

整理得：

${\mathrm{Ez}}_{i,j,k}^{n+1}=\frac{(2\mathrm{\ϵ}-\mathrm{\σ\Δt})}{(2\mathrm{\ϵ}+\mathrm{\σ\Δt})}{\mathrm{Ez}}_{i,j,k}^n+\frac{2\mathrm{\Δt}}{(2\mathrm{\ϵ}+\mathrm{\σ\Δt})}(\frac{{\mathrm{Hy}}_{i,j+1/2,k}^{n+1/2}-{\mathrm{Hy}}_{i,j-1/2,k}^{n+1/2}}{\mathrm{\Δx}}-\frac{{\mathrm{Hx}}_{i,j,k+1/2}^{n+1/2}-{\mathrm{Hx}}_{i,j,k-1/2}^{n+1/2}}{\mathrm{\Δy}})---\left(4\right)$

如果三个方向网格步长相等，即△x＝△y＝△z＝△h，上式可简写成：

${\mathrm{Ez}}_{i,j,k}^{n+1}=A\·{\mathrm{Ez}}_{i,j,k}^n+B\·({\mathrm{Hy}}_{i,j+1/2,k}^{n+1/2}-{\mathrm{Hy}}_{i,j-1/2,k}^{n+1/2}-{\mathrm{Hx}}_{i,j,k+1/2}^{n+1/2}+{\mathrm{Hx}}_{i,j,k-1/2}^{n+1/2})---\left(5\right)$

其中它们是空间分布的介电常数和电导率的函数。

类似地：Hz可离散化成：

${\mathrm{Hz}}_{i+1/2,j+1/2,k+1/2}^{n+1/2}=C\·{\mathrm{Hz}}_{i+1/2,j+1/2,k+1/2}^{n-1/2}+D\·({\mathrm{Ey}}_{i,j,k+1/2}^n-{\mathrm{Ey}}_{i,j,k+1/2}^n-{\mathrm{Ex}}_{i+1,j,k+1/2}^n+{\mathrm{Ex}}_{i,j,k+1/2}^n)---\left(6\right)$

其中它们是空间分布的磁导率和等效磁阻率的函数。方程组(2)中其余方程的离散化与(5)和(6)式无原则性差别，这里不再赘述。

FDTD方法是以差分方程组的解代替偏微分方程组的解，只有离散后差分方程组的解是收敛和稳定的，这种代替才有意义。为此，需满足如下Courant稳定性条件:

$\mathrm{\Δt}\≤\frac{1}{v_{\max }}{[{\left(\frac{1}{\mathrm{\Δx}}\right)}^2+{\left(\frac{1}{\mathrm{\Δy}}\right)}^2+{\left(\frac{1}{\mathrm{\Δz}}\right)}^2]}^{-1/2}---\left(7\right)$

其中是计算域的最大波速。

上述引出了FDTD基本计算方法，本发明要研究的是实现探地雷达三维大规模正演模拟，对于大尺度，计算空间相对于电磁波波长来说较大，单个测点的计算规模成了首要问题，需要研究其并行化策略，进行高效计算。下面将重点介绍本发明GPR正演计算的并行化策略。

目前对并行化策略有很多种分类方法，大致有：主从式、单程序流多数据流、数据流水线、分治策略和混合方法。GPR正演并行计算是一个复杂的过程，通过一种并行化策略往往不能解决所有的问题，本发明采用混合的方法，即一个程序包含了不同的并行化策略。

首先地下物体物性参数的分配、任务分配和最终数据的收集是采用主从式方法，如图3所示。主从式方法又称作任务播种，其基本思想是主节点将任务分解，分配到各个子节点，各节点自行计算，主节点然后负责收集各个子节点的求解结果最后汇总得到问题的最总结果。

FDTD是通过分治策略来进行并行计算，即通过区域分割，每个计算节点分配一个区域来实现大尺寸的计算，每个节点各自管理存储系统，各节点之间通讯来交换数据。图4给出了一维、二维、三维三种不同的任务分割方式。选择何种分割方式与计算区域的形状有关，同时会导致不同的节点间通讯量。例如图4相同大小的正方形计算域，一维分割中间有七个相同大小的数据交换面，二维分割和三维分割分别只有四个和三个数据交换面。一维分割往往会导致较大的通讯量，降低计算效率。我们采用三维的任务分割方式，最大限度地发挥并行机的计算效率。

各个节点间使用相同的程序，但是操作在各自不同的数据上，这称之为单程序多数据流SPMD方法，这种方法首先需要将应用程序的数据预先分配给各个计算节点，然后各节点并行地完成各自的计算任务，包括计算过程中的各节点的数据交换(施行通信同步)，最后将结果汇总起来。整个流程如图5所示。

各节点在计算过程中以及计算完成后，需要进行数据交换，三维并行FDTD的数据交换方式有多种，可以同时交换边界重合网格的电场和磁场，也可以只交换一种场分量。本发明采取交换磁场的方式，图6给出了节点间通讯交换数据的示意图。相邻节点的计算区域均包含一个重合网格，当一个节点进行FDTD电场迭代计算时，需要用到边界磁场，则把相邻节点该位置上的磁场值传递给该节点。同理，也需要把该节点上的边界磁场传递给相邻节点使得相邻节点边界电场的迭代计算能持续下去。采取本发明的并行策略能很大程度减小程序的复杂性，并提高计算效率。

对于节点间通讯，我们采用应用最广泛的MPI(MessagePassingInterface)并行库。MPI是基于消息传递的并行运算模式，是美国Argonne国家实验室于1992-1994年开发的,经过多个组织不断修改完善，现在已经成为并行计算领域事实上的标准之一。它有如下特点：(1)标准性，MPI是唯一的并行数据传输标准并支持所有的计算机系统平台；(2)可移植性，当运行应用程序在不同的平台上时，不需要修改源程序；(3)性能优化，可根据硬件设备进行优化。MPI包括超过200个函数，功能强大，对商业和科学研究均免费使用。

GPR正演计算的计算空间是有限的，为了模拟电磁波在地下和空气中无反射自由传播，需要在计算空间外围设置吸收边界，用于吸收传播到边界的电磁波，同时保证在边界上反射回计算空间的电磁波足够小。设置有吸收边界的电磁波传播就如同开放边界一样，不会被计算域固有的外边界反射。一般常用的吸收边界有二阶Mur和完全匹配层(PerfectlyMarchedLayer,PML)，其中二阶Mur反射较大，尤其是大角度入射时。PML原理是通过设置一定厚度的具有电导率、磁导率的吸收层来衰减电磁波，只要电导率和磁导率满足一定关系，PML吸收边界从理论上是完全无反射的，不受到电磁波入射角度限制。但在匹配层的实现过程中，由于空间的离散化等原因，理论上的完全无反射无法实现。本发明是设置从内到外按指数分布的电导率、磁导率，让电磁波逐层衰减。

在GPR正演计算过程中，需要用到激励源。激励源的设置对于GPR数值模拟至关重要。本发明选用超宽带Recker子波，它的表达式如下：

$A\left(t\right)=A_0[1-2{\left(\mathrm{\π}{f}_{m}t\right)}^2]e^{-{\left(\mathrm{\π}{f}_{m}t\right)}^2}$

其中f_m是超宽带激励源中心频率，A₀是脉冲输入最大值。它的波形如图7所示，它具有正负脉冲值。当中心频率f_m＝1.6GHz时，Recker子波的频率谱如图8所示，其半高宽的高频约为f_H＝2.62GHz，低频约为f_L＝0.78GHz，带宽为1.84GHz。

激励源时间离散化后按时间步加入，为确保离散化后频谱基本不变，以及考虑到FDTD本身的稳定性，时间步长必须足够小。因此整个GPR正演需要的时间较长，往往要数天或数周才能完成所有的测点计算。

验证：

上述对本发明基于FDTD探底雷达三维正演模拟方法做了介绍，下面我们通过并行程序的加速比来测试本发明方法的实施效果：

并行计算的加速比指的是同一个算法在单处理器上运算的时间除以在多个处理器上运行的时间。FDTD算法具有优良的并行化特征，对于三维问题，每个节点进行三维立体区域的计算，而和其他节点通讯时只交换边界二维平面的数据，因此可以得到较高的加速比和效率。我们采用Linux集群，由于条件限制，仅使用64台刀片服务器作为计算节点，另外配置2个系统管理节点和2个存储节点(2TB)，每台服务器配备2颗3.0GHzXeon处理器，峰值计算能力每秒超过2万亿次。图9给出了并行程序用不同节点数计算的加速比。从图可以看到，随着节点数的增加，加速比基本上是线性变化的，这说明该方法具有较好的可扩展性。

并行效率是指加速比除以节点数，对于分布式存储系统，节点间通讯量、同步、等待是影响效率的主要因素。从图9中可以看出，尽管随着节点的增加，效率有所下降，但效率依然很高，可以保持94％以上。对于节点计算能力相同的情况下，并行效率主要由并行机的硬件数据交换能力决定。随着计算机技术的发展，网络带宽和交换能力不断提高，并行效率可以进一步提高。

本发明基于FDTD对探地雷达进行大规模三维正演模拟方法，能够适用于各种尺度复杂地电介质和模型的计算，具有通用性强的特点，该方法特别适合于大尺度、高精度的正演模拟，有望提高探地雷达探测的可靠性、准确度及雷达剖面的分辨率，更有利于探地雷达资料的解释。将此方向用于GPR三维探测数据处理中，可以提高探测的精度和速度，促进GPR三维探测技术的实用化。

上述实施例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围，凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。

Claims (8)

1.一种基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：它包括如下步骤，

（1）采用基于消息传递的并行运算模式的MPL并行库，建立并行计算节点拓扑，确定节点通讯组的大小，将目标空间区域离散化，形成多个子区域，每个子区域对应一个节点，且在子区域上根据精度要求设置多个测点；

（2）初始化MPL并行库，并根据正演需求设置计算域大小、网格精度、时间总步长、激励源参数、地电模型、地下介质模型、吸收边界参数，根据这些参数计算需要的时间步长、空间步长这两项中间数据，获得当前节点的ID，并在各节点分配内存空间并初始化；

（3）通过网络进行各节点数据的交换，实现各节点通讯和同步；选择任一个测点基于FDTD正演对所有节点并行计算，首先计算电场值，电场计算完成后加入激励源计算磁场，通过通讯、电场计算、磁场计算不停的迭代，计算到设定的时间总步长后停止，保存输出的数据，该数据就是这个测点的雷达反射信号；

（4）步骤（3）完成后，对包括电场和磁场以及天线激励源的所有场值进行初始化更新，使得整个计算域恢复到计算开始前的状态，然后移动发射天线和接收天线的位置到下一个测点重新开始正演计算；

（5）若所有测点计算完成，把每个节点的数据汇集以进行成像和图像后处理；若测点未计算完成，重复进入步骤（3）（4）。

2.根据权利要求1所述的基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：步骤（4）中，所述的正演计算中，地下物体物性参数的分配、任务分配和最终数据的收集是采用主从式方法，将所述的节点设为若干子节点和一个主节点，所述的主节点将任务分解，分配到各个子节点，各子节点自行计算，所有测点计算完成后，把每个子节点的数据汇集到主节点，在主节点进行成像和图像后处理。

3.根据权利要求2所述的基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：所述的FDTD通过分治策略来进行并行计算，即通过区域分割，每个计算节点分配一个或大于一个子区域来实现大尺寸的计算，每个节点并行地完成各自的计算任务，包括计算过程中的各节点的数据交换，最后将结果汇总起来。

4.根据权利要求3所述的基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：所述的区域分割采用三维分割的方式。

5.根据权利要求1所述的一种基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：步骤（3）中，所述的数据交换采取交换磁场的方式，相邻节点的计算区域均包含一个重合网格，当一个节点进行FDTD电场迭代计算时，把相邻节点的边界磁场传递给该节点作为其边界磁场。

6.根据权利要求1所述的一种基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：在执行步骤（4）前，还需要在计算空间外围设置用于吸收传播到边界的电磁波的吸收边界。

7.根据权利要求6所述的一种基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：所述的吸收边界设置方式为：按照指数分布的方式从内到外设置具有电导率、磁导率的吸收层，以使得电磁波逐层衰减。

8.根据权利要求1所述的一种基于FDTD的探地雷达大规模三维正演模拟方法，其特征在于：步骤（4）中，所述的激励源的设置用超宽带Recker子波。