一种基于最小二乘克希霍夫偏移水工隧洞探地雷达数据成像方法与流程

1.本发明涉及水工隧洞探地雷达工程技术领域，具体涉及一种基于最小二乘克希霍夫偏移水工隧洞探地雷达数据成像方法。


背景技术：

2.近年来，国内外科学家围绕水工隧洞探地雷达数据开展了大量的处理和解释工作。其中有代表性的研究工作包括：吴丰收2009年利用探地雷达对钢筋混凝土内钢筋及缺陷进行综合探测研究分析，利用偏移技术聚焦钢筋的离散杂波。wu等人提出超曲波变换技术来压制钢筋信号。xiao等人提出利用多带通滤波技术抑制混凝土钢筋gpr 剖面中由钢筋等近地表周期性元素引起的杂波。朱自强等人于2017 年在gpr探测钢筋混凝土中引入时频分析方法s变化，可以更好的对钢筋以及缺陷进行判断且成像。zhang和liu等人使用逆时偏移结合改进的全波形反演准确地做到gpr探测中隧道衬砌和围岩中的病害识别。
3.在地质结构比较完整，速度比较简单，没有明显的横向空间变化时，探地雷达数据处理可以考虑使用stolt偏移或者fk偏移即可，在地质条件比较复杂的区域，stolt偏移或者fk偏移则不能精确的计算，但嫦娥三号数据是共偏移距数据，基于波动方程的逆时偏移等高精度算法不能应用于水工隧洞探地雷达数据，而常规的stolt偏移或者fk偏移也无法实现数据的处理，因此如何实现在水工隧洞探地雷达数据进行高精准的处理，是本领领域技术人员亟需解决的问题。


技术实现要素：

4.针对上述存在的问题，本发明提出了一种基于最小二乘克希霍夫偏移水工隧洞探地雷达数据成像方法，主要解决现有技术中存在的常规的stolt偏移或者fk偏移也无法实现探地雷达数据处理的问题。
5.为了实现上述的目的，本发明采用以下的技术方案：
6.一种基于最小二乘克希霍夫偏移的水工隧洞探地雷达数据处理方法，包括如下步骤：
7.s1)读入水工隧洞探地雷达数据，并对数据进行预处理；
8.s2)进行速度谱分析，确定水工隧洞内部主要层结构的位置；
9.s3)开展克希霍夫的水工隧洞探地雷达数据偏移处理；
10.s4)加入最小二乘迭代进行数据偏移算法；
11.s5)输出水工隧洞探地雷达数据成像结果。
12.优选的，步骤s1)中数据预处理包括以下过程：
13.(1)数据读取：根据标准存储格式读取各类商业探地雷达数据格式，并逐卷读取道集及位置信息；
14.(2)数据拼接：将读取到的数据拼接成一组完整剖面数据；
15.(3)道集调整：对得到的剖面数据根据连续同相轴进行道集调整；
16.(4)道集择取：将调整后的道集进行筛选，选择数据质量好的道集；
17.(5)时间差调整；
18.(6)无用数据删除：将调整后的数据信噪比低的无用数据部分删除；
19.(7)带通滤波：进行带通滤波，去除数据中的低频及高频噪声；
20.(8)背景去除；
21.(9)自动增益。
22.优选的，步骤s2)中进行速度谱分析，确定水工隧洞内部主要层结构的位置是通过转换数据由共偏移距道集为中心点数据开展动校正，在水平迭代处理过程中开展速度谱分析，通过时频分析技术，计算出水工隧洞内部主要雷达信号反射层的位置。
23.优选的，步骤s3)中进行克希霍夫的水工隧洞探地雷达数据偏移处理的主要流程如下：
24.线性正演模拟算子l
～
满足
25.p＝l
～mꢀꢀꢀꢀꢀꢀꢀꢀ
公式(1)
26.其中，p表示建模数据的矢量，m表示模型反射率模型矢量，克 mk＝l
～t
p希霍夫偏移使用公式(1)中的正演模型算子的转置：
27.mk＝l
～t
p
ꢀꢀꢀꢀꢀꢀꢀ
公式(2)
28.将公式(1)带入公式(2)，可得
29.mk＝l
～t
l
～mꢀꢀꢀꢀꢀꢀꢀꢀ
公式(3)
30.其中，l
～t
l
～
表示hessian矩阵，对应l
～t
l
～
i的积分算子的显式形式表达公式为：
31.k(x，x
′
)＝∫dshs(s)a
sxasx
′
×
∫drhr(r)a
xrax
′rr(τ
sx
τ
xr-τ
sx
′-τ
x
′r)
ꢀꢀꢀꢀꢀ
公式(4)
32.其中，x和x
′
分别表示与mk和m相关的位置，τ
sx
表示从源s处到x处反射器的波传播时间，τ
xr
表示从反射器到r处接收器的波传播时间，a
sx
和a
xr
表示是振幅项，r表示时间(τ
sx
τ
xr-τ
sx
′-τ
x
′r) 延迟的相关时间交叉系数，hs(s)和hr(r)分别表示源和接收器采样的函数。
33.优选的，l
～t
l
～
的积分算子的显式形式中主对角线的表达公式为：
[0034][0035]
其中，k(x0，x0)表示公式(4)中在x0点的表达形式，a
sx
和a
xr
表示是振幅项，r(0)表示零时间延迟的相关时间交叉系数，hs(s)和 hr(r)分别表示源和接收器采样的函数。
[0036]
优选的，步骤s4)中加入最小二乘迭代进行数据偏移算法，其中，克希霍夫偏移将在偏移中产生偏移畸变，并通过下列公式求解偏移畸变量m：
[0037]
p(m)＝||l
～
m-p0||2 ε2||c
～
m-c
～mapr
||2ꢀꢀꢀ
公式(6)
[0038]
其中，||l
～
m-p0||2表示数据失配函数，ε2||c
～
m-c
～mapr
||2表示正则化项，ε2表示正则化权重，c
～
表示作用于m的线性算子；
[0039]
最小化公式(6)的模型表达公式为：
[0040]
m＝(l
～t
l
～
ε2c
～tc～
)-1
(l
～t
p0 ε2c
～tc～mapr
)
ꢀꢀꢀꢀꢀ
公式(7)
[0041]
求解公式(7)：
[0042]
(l
～t
l
～
ε2c
～tc～
)m＝l
～t
p0 ε2c
～tc～mapr
ꢀꢀꢀꢀꢀ
公式(8)
[0043]
最后采用共轭梯度方案求解关于m的公式(8)，求解的m即是最后的克希霍夫偏移成像输出结果。
[0044]
由于采用上述的技术方案，本发明的有益效果是：本发明偏移计算效率高，能够适合于水工隧洞探地雷达的数据观测系统，并且对输入的雷达数据没有特殊的要求，不依赖于速度模型，处理方式方便灵活，适合于目标成像，同时采用基于最小二乘迭代方式的克希霍夫偏移可以提高目标成像分辨率和振幅保真度，在水工隧洞探地雷达数据处理和目标成像定位等方面具有很好的应用价值；
[0045]
本发明的数据预处理部分采用完整的水工隧洞探地雷达数据处理流程，包括去除重复数据，滤波，增益，去坏道集等过程。保证下一步过程克希霍夫偏移输入探地雷达数据有较高的信噪比和信号强度，同时为提高水工隧洞探地雷达数据偏移成像精度，进行了克希霍夫偏移成像技术，能够适应不同的观测系统，对输入地震数据没有特殊要求，处理方式方便灵活，非常适于做目标成像；其次，计算效率较高，对于复杂构造，横向起伏较大的地层分布，避免了常规探地雷达fk偏移方法难以获得准确的成像结果，在起伏较大的构造边缘存在绕射干扰难以消除等问题。
附图说明
[0046]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的，保护一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0047]
图1为本发明的流程图；
[0048]
图2为本发明的探地雷达数据预处理流程图；
[0049]
图3为本发明的理论水工隧洞探地雷达数据介电常数模型a)、常规fk偏移成像结果b)；
[0050]
图4为本发明的克希霍夫偏移成像结果a)、最小二乘克希霍夫偏移成像结果b)；
[0051]
图5为本发明的概念模型a)、克希霍夫偏移成像结果b)、最小二乘克希霍夫偏移成像结果c)；
[0052]
图6为本发明的常规fk偏移方法偏移成像结果a)、水工隧洞探地雷达数据克希霍夫偏移成像结果b)、水工隧洞探地雷达数据最小二乘克希霍夫偏移成像结果c)、水工隧洞探地雷达数据最小二乘克希霍夫偏移成像结果d)。
具体实施方式
[0053]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0054]
如图1-6所示，一种基于最小二乘克希霍夫偏移的水工隧洞探地雷达数据处理方法，由于发射天线与接收天线启动时间不一致，从而进行时间差调整；删除信息极弱且信噪比较低的无用数据部分；然后进行带通滤波，去除低频及高频噪声；由于水工隧洞表面反射
强烈，存在水平连续的干扰波，需要减去平均道，突出异常信息；同时为了突出下部信息，需要进行自动增益(agc)，通过转换数据由共偏移距道集为中心点道集，通过开展动校正，水平迭代处理过程开展速度谱分析，通过时频分析技术，计算出水工隧洞内部主要雷达信号反射层的位置。
[0055]
在开展克希霍夫的水工隧洞探地雷达数据偏移处理时，其克希霍夫偏移的主要流程如下：
[0056]
线性正演模拟算子l
～
满足
[0057]
p＝l
～mꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
公式(1)
[0058]
其中，p表示建模数据的矢量，m表示模型反射率模型矢量，观测数据p0由p0＝l
～0
m0描述，其中m0是真正的地球反射率模型向量，而l
～0
是实际地球模型的正演模拟算子，除非另有说明，否则假设 l
～
＝l
～0
；克希霍夫偏移使用公式(1)中的正演模型算子的转置：
[0059]
mk＝l
～t
p
ꢀꢀꢀꢀꢀꢀꢀ
公式(2)
[0060]
将公式(1)带入公式(2)，可得
[0061]
mk＝l
～t
l
～mꢀꢀꢀꢀꢀꢀ
公式(3)
[0062]
其中，矩阵l
～t
l
～
表示hessian矩阵，并将mk定义为m的l
～t
l
～
滤波版本，克希霍夫kirchhoff偏移算子将正确地重建实际的地球模型向量，通常情况下l
～t
l
～
不是单位矩阵，并具有主对角线元素是不均匀且不统一，同时主对角线之外的元素是非零的特征，对应l
～t
l
～
的和分算子的显式形式表达公式为：
[0063]
k(x，x
′
)＝∫dshs(s)a
sxasx
′
×
∫drhr(r)a
xrax
′rr(τ
sx
τ
xr-τ
sx
′-τ
x
′r)
ꢀꢀꢀꢀ
公式(4)
[0064]
其中，x和x
′
分别表示与mk和m相关的位置，τ
sx
表示从源s处到x处反射器的波传播时间，τ
xr
表示从反射器到r处接收器的波传播时间，a
sx
和a
xr
表示是振幅项，r表示时间(τ
sx
τ
xr-τ
sx
′-τ
x
′r) 延迟的相关时间交叉系数，hs(s)和hr(r)分别表示源和接收器采样的函数。
[0065]
可以从公式(4)中观察到：l
～t
l
～
的积分算子的显式形式中主对角线的表达公式为：
[0066][0067]
该表达式表明，主对角线元素(对于给定的子波源)的值取决于扩散损失和由hs(s)和hr(r)决定的波场的不连续采样。 kirchhoff偏移的scheme通常会被增加以弥补扩散损失；bleistein (1984)等的分析反演公式由此产生。通过在每个cdp箱中的数字命中计数对k(x0，x0)归一化，可以减轻不连续采样的影响。该补偿方案可以被解释为将对角矩阵(l
～t
l
～
)-1应用为预处理矩阵。
[0068]
在建模和偏移之后，l
～t
l
～
每列都是对偏移部分mk的m中固定元素的响应，对于一道，偏移响应会沿着由x描述满足τ
sx
τ
xr
＝τ
sx
′
τ
x
′r的偏移椭圆，对于完整数据，偏移响应是对各道的总响应(以各道相对应的加权系数相加)，并且对于x
′
处的每个衍射点，偏移产生在x＝x
′
附近的聚焦图像，在这种情况下，l
～t
l
～
是对角占优矩阵l
～t
≈l-1
，对于不完整的数据，偏移响应不只围绕x＝x
′
聚焦，因为缺失的数据导致偏移畸变的不完全消除，因此在l
～t
l
～
中存在明显的非对角元素。在这种情况下，l
～t
≈l-1
，需要某种矩阵求逆方法。
[0069]
在加入最小二乘迭代进行数据偏移算法时，其中，克希霍夫偏移将在偏移中产生
偏移畸变，并通过下列公式求解偏移畸变量m：
[0070]
p(m)＝||l
～
m-p0||2 ε2||c
～
m-c
～mapr
||2ꢀꢀꢀ
公式(6)
[0071]
公式(6)右边的第一项是数据失配函数，第二项是正则化项，c
～
表示作用于m的线性算子，使用先验信息为每项指定c
～
，模型的先验信息可以与m
apr
项合并，其中是正则化权重，公式(6)当c-滤波模型残差和数据残差分布是高斯分布，则最小化公式(6)的模型表达公式为：
[0072]
m＝(l
～t
l
～
ε2c
～tc～
)-1
(l
～t
p0 ε2c
～tc～mapr
)
ꢀꢀꢀꢀ
公式(7)
[0073]
求解公式(7)：
[0074]
(l
～t
l
～
ε2c
～tc～
)m＝l
～t
p0 ε2c
～tc～mapr
ꢀꢀꢀꢀ
公式(8)
[0075]
最后采用共轭梯度方案求解关于m的公式(8)，求解的m即是最后的克希霍夫偏移成像输出结果。
[0076]
其中，公式4是为了计算l
～t
l
～
是hessian矩阵，其显示形式用k 表示(即公式4)，公式4是其一般形式，在偏移方法中为了简化，即只计算公式5中的x0点的k值，公式6是计算最后偏移输出量m 的表达式，7和8是求解的过程，而其中所需的l
～t
l
～
是hessian矩阵是通过公式5求得。
[0077]
附图3是张金海等(pnas，2014)采用fk偏移方法对理论模型偏移成像结果，从中我们可以看出存在很多绕射干扰，附图4是采用本技术的克希霍夫偏移a)和最小二乘克希霍夫偏移b)方法对应附图3中相同模型的偏移成像结果，结合附图3和附图4可以得出本发明方法有更高的成像精度，对于各个目标体周围的干扰能很好的消除。
[0078]
附图5中a)是采用法文哲(grl，2015)论文中的月壤结构解释模型测试本技术中克希霍夫偏移(图b)和最小二乘克希霍夫偏移 (图c)的效果，通过月壤概念模型也验证了本发明方法能对月壤内部结构得到很好的成像结果。
[0079]
附图6是采用本技术的方法对水工隧洞探地雷达实测数据开展的偏移成像，图6a)和图6b)是原始水工隧洞探地雷达数据和去背景数据，图6c)和6d)是常规偏移和本技术采用的方法结果，可以看出本发明的结果在实测数据应用中也有相当较好的偏移成像结果。
[0080]
在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料过着特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0081]
以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。