一种相关双向波场拟合高精度探地雷达逆时偏移成像方法

本发明涉及地球物理及其数据处理，具体涉及一种用于较大规模gpr数据的加载和处理的相关双向波场拟合高精度探地雷达逆时偏移成像方法。

背景技术：

1、探地雷达是一种利用高频电磁波探测地下结构内部物质分布规律的无损检测方法，具有连续、快速、抗干扰能力强等优点，广泛应用于浅层地下结构异常体探测定位。通过探地雷达探测剖面可以粗略了解地下目标体的分布情况，然而，在复杂的地下结构剖面中异常体绕射波和反射波交错叠加，仅通过探测剖面无法明确地下结构内部的真实情况。偏移作为实测数据向地下真实结构转换的重要技术，能够清晰地反映复杂地下结构内部构造及异常体的位置等细节特征。目前，常用的偏移方法主要有kirchhoff偏移、f-k偏移、逆时偏移(rtm)等。其中，逆时偏移具有理论简单，成像精度高、无倾角限制等优点，广泛应用于地下结构异常体gpr探测成像定位。然而，传统的rtm算法在精度和效率上仍存在一些不足之处，无法满足实际工程的需求。

2、与此同时，探地雷达逆时偏移计算量大，计算过程中需要多次调用计算密集且耗时的正演模拟，常用的fdtd正演模拟算法简单明了，结构清晰，但在计算效率仍有些不足之处。

3、另一方面，传统的rtm成像算法采用零延时互相关成像条件会产生低频噪声和伪影影响成像结果的清晰度，因此去除图像中的噪声和伪影也是当下亟需解决的另一个关键难题。

技术实现思路

1、本发明要解决的技术问题是针对现有技术的种种不足，提供一种相关双向波场拟合高精度探地雷达逆时偏移成像方法。

2、为解决上述技术问题，本发明所采取的技术方案如下。

3、相关双向波场拟合高精度探地雷达逆时偏移成像方法，对于大规模gpr数据的加载和处理，引入基于gpu的多线程的并行算法路线，一方面构建算法进行正向波场模拟和逆向波场的拟合并进行互相关成像，进一步构建基于多参数的适应性图形优化算法，以获得针对不同类型的地质结构的高精度探地雷达逆时偏移图像。

4、作为本发明的一种优选技术方案，构建正向波场模拟和逆向波场的拟合算法及其互相关成像算法时，引入波场的能量守恒及其相位准确性这两个保守数据要素作为探地雷达正演方程的约束条件，由此减少误差积累和提高图像质量。

5、作为本发明的一种优选技术方案，构建正向波场模拟和逆向波场的拟合算法及其互相关成像算法，一个具体化的数据进程包括如下步骤：

6、s1.构建探地雷达正演方程；

7、s2.加载激励源信号和背景模型参数，利用gpu并行计算t＝0时刻模拟传播至t＝max时刻的正向波场，并将其存储在内存中；

8、s3.加载需要偏移的gpr剖面将其作为激励源信号，利用gpu并行计算从t＝max时刻沿时间轴反向延拓至t＝0时刻的逆向波场；

9、s4.在逆向延拓同时，从内存中取出相应时刻的正向波场，随后与逆向波场进行gpu并行互相关成像；

10、s5.重复s2-s4步骤，待所有单道数据偏移完成后进行叠加，得到原始偏移图像。

11、作为本发明的一种优选技术方案，步骤s1中，正演方程具体构建为：

12、二维各向同性介质中，maxwell方程构建为如下形式：

13、

14、

15、式中，h和e分别代表磁场和电场向量，σ，ε和μ分别表示电导率，介电常数和磁导率；

16、引入矢量磁位并令e＝-u，二维tm波的maxwell方程组表示为如下正则方程形式：

17、

18、

19、式中，az和uz分别表示场组分a和u沿坐标轴z方向的分量，▽2为拉普拉斯算子；

20、应用辛算法离散上式得到如下迭代计算公式：

21、

22、

23、式中，和分别表示场组分az和uz在ndt时刻，空间网格节点(i,j)上的离散值；

24、作为本发明的一种优选技术方案，步骤s2中，采用gpu并行运算正向延拓波场，具体的，在探地雷达正演模拟中更新场分量az和uz，其中每个网格处场分量az和uz的更新是相互独立且同步进行的，具有高度并行性，基于此引入gpu进行并行加速计算，调用大量线程同时更新所有网格点；

25、作为本发明的一种优选技术方案，步骤s3-s4中，构建如下计算策略：互相关成像条件需要同一时刻的正向延拓波场和逆向延拓波场，这两个波场一个是在时间的正向上进行而另一个是在时间的反向上进行，由此在正向延拓时存储全部时刻正向波场，在逆向延拓时从内存中取出相应时刻的正向波场，随后与逆向延拓波场进行互相关成像。

26、作为本发明的一种优选技术方案，构建基于多参数的适应性图形优化算法时，基于系列化可选用/可调整正则化函数和/或正则化参数的优化算法，实现针对不同类型的地质结构成像进行图像噪声和伪影去除，从而获得高精度探地雷达逆时偏移图像。

27、作为本发明的一种优选技术方案，构建基于多参数的适应性图形优化算法，一个具体化的数据进程包括如下步骤：

28、s6.构建tv正则化函数，并构建交替方向乘子法进行数值求解；

29、s7.对原始偏移结果进行去噪去伪影处理得到精确度更高的最终偏移图像。

30、作为本发明的一种优选技术方案，步骤s6中，构建tv正则化函数的具体数据进程为：以x为去噪后的清晰rtm剖面，y是带噪声的原始rtm剖面，表示为：y＝x+n；其中，n是噪声；

31、对rtm剖面进行去噪处理，寻找一个合适的x使x和y之间的残差最小，引入l1-tv正则化项后，构建为：

32、

33、其中，λ为正则化参数，用来保持正则部分与保真部分的平衡；||x||tv为tv正则化项；表示垂直和水平方向上的有限差分算子，||||1代表l1范数；||||2代表l2范数；

34、进一步将无约束优化问题重新表述为以下约束优化问题：

35、

36、s.t.x＝z

37、其中，z为松弛变量(slack variable)，ξ(z)＝||z||tv；

38、上式增广拉格朗日形式可表示为：

39、

40、其中，u为增广拉格朗日乘子；ρ为非负惩罚参数，决定迭代重建的收敛度，<·>表示点积；最后构建算法求解l(x,z,u)进行优化。

41、作为本发明的一种优选技术方案，求解l(x,z,u)的具体数据进程包括：

42、(1)固定zk和uk，更新xk+1：

43、

44、上式的最优解可表示为：

45、(2)固定xk+1和uk，更新zk+1：

46、

47、上式可等价写为：

48、其中，

49、(3)固定xk+1和zk+1，更新uk+1：

50、uk+1＝uk+ρ(zk+1-xk+1)

51、首先输入图像数据y，设置合适的λ，ρ，并规定一定迭代次数和限定误差，令初始值x0＝z0＝y，u0＝0，然后开始循环执行(1)，(2)，(3)，当达到最大迭代次数和限定误差后，算法停止，最后输出tv正则化结果x；数据进程完成。

52、采用上述技术方案所产生的有益效果在于：本发明创新性的运用一阶算法构建了探地雷达正演方程，对于较大规模的数据集引入了并行计算架构，给定激励源信号后利用gpu并行计算并存储t＝0时刻至t＝max时刻的正向延拓波场，同时对给定的接收剖面利用gpu并行计算t＝max时刻至t＝0时刻的逆向延拓波场，在逆向延拓同时从内存中取出相应时刻的正向波场并随后与逆向延拓波场进行互相关成像；该算法首次被引入逆时偏移中，通过提升正演模拟计算效率的从而提高rtm计算效率，在模拟过程中不会引入人为的耗散机制和虚假的激励稳定性较高，并在保证精度的同时显著提升了数据处理效率；

53、进一步，本发明还构建了全变差正则化方法对图像做进一步的去噪去伪影处理，所构建的数据进程易于与图像特征信息相结合，在gpr rtm剖面处理中构建了交替方向乘子法范式下的正则化去噪策略，形成高质量、高精度的rtm剖面，不仅能抑制图像的噪声，而且能保留图像的纹理、边缘等细节，能够有效提升探地雷达逆时偏移精度及其综合成像质量。