CN110095773A - 探地雷达多尺度全波形双参数反演方法 - Google Patents

探地雷达多尺度全波形双参数反演方法 Download PDF

Info

Publication number
CN110095773A
CN110095773A CN201910477198.5A CN201910477198A CN110095773A CN 110095773 A CN110095773 A CN 110095773A CN 201910477198 A CN201910477198 A CN 201910477198A CN 110095773 A CN110095773 A CN 110095773A
Authority
CN
China
Prior art keywords
model
inverting
parameter
objective function
inversion
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910477198.5A
Other languages
English (en)
Other versions
CN110095773B (zh
Inventor
冯德山
王珣
李广场
刘强
曹岑
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Zhejiang East China Engineering Safety Technology Co Ltd
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang East China Engineering Safety Technology Co Ltd, Central South University filed Critical Zhejiang East China Engineering Safety Technology Co Ltd
Priority to CN201910477198.5A priority Critical patent/CN110095773B/zh
Publication of CN110095773A publication Critical patent/CN110095773A/zh
Application granted granted Critical
Publication of CN110095773B publication Critical patent/CN110095773B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种探地雷达多尺度全波形双参数反演方法，包括将观测数据的GPR时间域共剖点记录作为输入，加载子波，建立初始模型并设定迭代精度；对建立的初始模型进行正演；确定多频率反演顺序；构造Wiener滤波器；采用Wiener滤波器滤波后得到起始频带和原始频带；数据在低频带进行反演再在高频带进行反演；构建反演目标函数；对目标函数进行求解；引入伴随场；计算更新步长并更新当前模型；以得到的模型作为新的初始模型重复上述步骤直至反演结束得到最终的反演结果。本发明能够提高雷达反演剖面的分辨率，从而提高GPR资料解释准确度，而且可靠性高，精确度高。

Description

探地雷达多尺度全波形双参数反演方法
技术领域
本发明具体涉及一种探地雷达多尺度全波形双参数反演方法。
背景技术
探地雷达(Ground penetrating radar，GPR)是一种对地下异常体的结构、特性或物体内部不可见目标体进行定位的电磁无损探测技术。实际的GPR探测中仅仅依靠获取的雷达剖面，即可大致推断探测目标体的位置、大小与尺寸，但多局限于粗略估计，尚不能给出异常体的介电常数、电导率等参数，从而精准界定目标体的属性。
为了获得地下目标体的更精确信息，实现高精度成像，层析成像、速度分析及全波形反演等速度建模方法被提出。Williamson&Worthington，Zhou&Liu应用射线类方法进行层析成像，由于只利用了信号信息的一小部分，仅能对大于信号主波的异常体成像，分辨率近似为第一菲涅尔带的直径，成像精度有限；Cai等利用二维Maxwell方程进行电磁波速度反演，其结果优于射线追踪成像，但计算量增加若干数量级。而时域全波形GPR反演是新兴起的重要成像手段，它起源于近代地震数据处理领域，相比于基于射线法的旅行时层析成像法，全波形反演法直接利用振幅、走时和相位等雷达波场信息，将正演的雷达记录与实际观测的雷达波形进行匹配，通过使两者的数据之差最小建立目标函数，从而寻找最佳模型参数，将分辨率提高到波长的四分之一。显然，全波形反演能够获取介质的介电常数、电导率等物性参数，有助于岩性的精确描述和异常体预测，提高GPR解释精度。
目前全波形反演已在钻孔雷达成像中得到广泛关注与研究，分为时间域全波形钻孔雷达数据反演及频率域全波形钻孔雷达数据反演。与钻孔GPR相比，地面GPR反演受测量方式制约，反演的不适定性更显著，计算量大，反演效率过低、易受初始模型约束，成像效果仍不甚理想。为了进一步提高地面GPR反演精度，众多学者开展了深入研究，取得了一系列成果：Moghaddam等开展小目标体介电常数波形反演成像研究；Cui等利用高阶Born近似研究了地下目标体反演技术；Rucker等用模拟退火法进行GPR资料的非线性反演，但算法的效率过低；Bouajaji等对地表GPR测量进行了地质解释，以实现对介质介电常数二维切片的定量成像；Lavoué等采用频率域的拟牛顿法对多偏移距GPR数据进行了二维介电常数和电导率的全波形反演。Feng等将地面GPR全波形反演与地震全波形反演相结合，对地下构造进行了定量成像。此外，全波形反演也被用于反演混凝土中氯离子导电梯度与湿度评估、土壤细粒土的改良特性研究等典型工程实例中，但受制于反演速度、源的提取、噪声等影响，距离实践推广仍存在一定的距离。
发明内容
本发明的目的在于提供一种能够提高雷达反演剖面的分辨率，从而提高GPR资料解释准确度的探地雷达多尺度全波形双参数反演方法。
本发明提供的这种探地雷达多尺度全波形双参数反演方法，包括如下步骤：
S1.将观测数据的GPR时间域共剖点记录作为输入，加载子波，建立初始模型并设定迭代精度；
S2.对建立的初始模型进行正演；
S3.确定多频率反演顺序；
S4.构造Wiener滤波器；
S5.采用步骤S4构造的Wiener滤波器进行滤波后，得到起始频带和原始频带；数据先在低频带进行反演，然后再在高频带进行反演；
S6.构建反演目标函数；
S7.对目标函数进行求解；
S8.引入伴随场；
S9.计算更新步长，并更新当前模型；
S10.以步骤S9得到的模型作为新的初始模型，重复步骤S2～S9直至反演结束，得到最终的反演结果。
步骤S2所述的对建立的初始模型进行正演，具体为采用基于交错网格的时域有限差分法进行GPR正演，同时时间离散精度和空间离散精度均设为2阶，截断边界处采用非分裂的卷积完全匹配层。
步骤S4所述的构造Wiener滤波器，具体为采用如下算式作为Wiener滤波器：
式中fwiener(ω)为频率域Wiener低通滤波器，Wtarget(ω)为目标激励源子波频谱，*为共轭符号，Woriginal(ω)为原始激励源子波频谱，δ为防止分母为零的常量小数。
步骤S6所述的构建反演目标函数，具体为采用如下步骤构建反演目标函数：
(1)构建数据目标函数Φd(m)：
式中i为源的编号，j为接收器的编号，M为源的个数，N为每个源的接收器个数，Ei(m,rj,t)为第i个源对猜测模型正演计算的模拟数据，m为模型介质参数向量且m＝(μ(r),ε(r),σ(r))T，μ为磁导率，ε为介电常数，σ为电导率，为第i个源激发在rj处接收到的观测数据，rj为第j个接收器空间坐标向量；
(2)引入相对电导率σr＝σ/σ0，取参考介质的电导率σ0＝1/η0；引入无量纲比例因子β，将模型参数m设定为相对介电常数和相对电导率的线性组合(εr,σr/β)；
(3)加入先验模型约束的正则化项；在引入一个全变差正则化后，形成新的目标函数Φ(m)：
Φ(m)＝Φd(m)+λΦm(m)
式中Φd(m)为数据目标函数，λ为正则化因子，Φm(m)为模型参数目标函数且Φm(m)＝R(εr)+R(σr/β)，Ω为模型空间，为微分算子，R为全变差正则化算子且
步骤S7所述的对目标函数进行求解，具体为采用局部求导类的有限内存BFGS法求解，每次迭代需计算目标函数的Fréchet导数：
式中δm＝(δε(r),δσ(r),δμ(r))T，
步骤S8所述的引入伴随场，具体为引入伴随场定义算子L*为算子L的伴随算子，根据伴随作用：
<L*w,δu>＝<w,Lδu>
式中<,>表示内积；
目标函数对参数μ，ε和σ的梯度分别为gμ，gε和gσ：
重写尺度变化之后的梯度向量为：
式中，ε0与η0为常量，β为参数调整因子；
全变差正则化算子的梯度可表示为：
目标函数Φ(m)的梯度写为：
步骤S9所述的计算更新步长，具体为采用非精确线搜索方法Wolfe准则计算更新步长。
本发明提供的这种探地雷达多尺度全波形双参数反演方法，提出了一种多尺度全变差正则化GPR全波形反演算法，并将该算法成功用于雷达数据的介电常数与电导率双参数同步反演中；通过引入参数调节因子，实现了相同的迭代步内同时更新介电常数和电导率值，结合并行方式在普通微机上实现了高精度的GPR反演；将多尺度串行反演策略应用到探地雷达全波形反演中，它将反演问题分解为不同尺度，先将雷达数据进行低频带反演，再将低频带的反演结果作为高频带反演的初始模型，通过采用2～3个低频带到高频带的逐频反演，从而逐步搜索到全局极值点，避免陷入局部极值，与单尺度的反演结果对比表明，多尺度串行反演具有更高的反演精度；在双参数全波形雷达反演中引入了参数调节因子β及正则化参数λ，通过选取合适的参数调节因子，能够有效解决介电常数与电导率在数值上相差较大的问题；而正则化参数的引入，它可以使数据拟合与模型约束之间达到更好的平衡，使反演更加快速稳定收敛。
附图说明
图1为本发明方法的方法流程示意图。
图2为本发明方法的模型1的参数分布示意图。
图3为本发明方法的模型1的100MHz中心频率的雷克子波的时域共剖点记录示意图。
图4为本发明方法的模型1的无噪数据反演结果示意图。
图5为本发明方法的相对介电常数一维深度剖面图。
图6为本发明方法的电导率的一维深度剖面图。
图7为本发明方法的模型1使用缩放因子β＝1.0和正则化权重λ＝5.0的无噪数据反演结果示意图。
图8为本发明方法的SNR＝25dB的白噪声下模型1的时域共剖点记录示意图。
图9为本发明方法的模型1使用缩放因子β＝1.0和减小的正则化权重的噪声数据(SNR＝25dB)反演结果示意图。
图10为本发明方法的“中”字型模型的参数分布示意图。
图11为本发明方法的模型2的200MHz中心频率的雷克子波的时域共剖点记录示意图。
图12为本发明方法的Wiener滤波操作示意图。
图13为本发明方法的三种不同反演策略的结果示意图。
图14为本发明方法的单道反演切线示意图。
图15为本发明方法的不同反演策略下目标函数及模型重构误差对比示意图。
图16为本发明方法的不同参数调节因子下目标函数及模型重构误差对比示意图。
图17为本发明方法的不同参数调节因子下反演结果示意图。
图18为本发明方法的不同正则化参数下反演结果示意图。
图19为本发明方法的不同参数调节因子及正则化参数下目标函数及模型重构误差对比示意图。
具体实施方式
如图1所示为本发明方法的方法流程示意图：本发明提供的这种探地雷达多尺度全波形双参数反演方法，首先，进行时域共剖点记录、子波及初始模型的加载；其次，采用Wiener低通滤波器将GPR观测数据与子波滤到从低频带到高频带的不同频率段；然后，引入参数调节因子及正则化因子，建立目标函数，并采用伴随状态法求取模型参数梯度公式；利用非精确线搜索方法计算更新步长，更新当前模型，直到达到反演终止条件，得到全波形反演的最终结果；具体包括如下步骤：
S1.将观测数据的GPR时间域共剖点记录作为输入，加载子波，建立初始模型并设定迭代精度；
S2.对建立的初始模型进行正演；具体为采用基于交错网格的时域有限差分法(FDTD)进行GPR正演，同时时间离散精度和空间离散精度均设为2阶，截断边界处采用非分裂的卷积完全匹配层；
探地雷达FDTD正演的公式体系如下：
由电磁场与电磁波的理论可知，二维GPR波满足的Maxwell方程可表示为：
Lu＝j
式中L为正演算子，u为波场向量，j为场源：
u＝(Hx Hz Ey)T
j＝(0 0 Jy)T
式中，上角标T表示转置，Hx、Hz为磁场强度分量(A/m)，Ey为电场强度分量(V/m)，Jy为电流密度分量(A/m2)。系数矩阵A，B，C，D为：
式中ε为介电常数(F/m)，μ为磁导率(H/m)，σ为电导率(S/m)；由于GPR反演中需要多次调用GPR正演，因此，正演的效率及并行加速的便利性非常关键，为此，本发明选用基于交错网格的时域有限差分法(FDTD)进行GPR正演，时间空间离散精度均为2阶，截断边界处采用非分裂的卷积完全匹配层；
S3.确定多频率反演顺序；
雷达数据中的低频分量主要包含地下较大的构造体信息，而高频分量包含较小构造体的细节信息。因此，本发明将高/低频分量结合起来进行多尺度反演，采用2-3个低频带到高频带的逐频带反演，根据不同尺度上的反演目标函数的特征去求解反演问题，从而逐步搜索到全局极值点，避免陷入局部极值；
S4.构造Wiener滤波器；
本发明采用Wiener低通滤波器，对观测数据与激励源子波滤波得到低频带信息：
式中fwiener(ω)为频率域Wiener低通滤波器，Wtarget为目标激励源子波频谱，Woriginal(ω)为原始激励源子波频谱，上标*表示共轭，δ为一个防止分母为零的常量小数，如果该值选取过大，会导致滤波子波形态与目标子波不同，这里设δ＝10-4。可以将数据变换到频率域，低通滤波到目标频段，再反变换回时间域。需要注意的是，震源函数和观测数据都要进行低通滤波；
S5.采用步骤S4构造的Wiener滤波器进行滤波后，得到起始频带和原始频带；数据先在低频带进行反演，然后再在高频带进行反演；
S6.构建反演目标函数；具体为采用如下步骤构建反演目标函数：
(1)构建数据目标函数Φd(m)
探地雷达全波形反演实质是利用已知的实测数据来重构地下介质中的模型参数：磁导率μ，介电常数ε与电导率σ的空间分布；根据正演模拟数据与实测数据之间的拟合最优，数
据目标函数可以定义为：
式中i为源的编号，j为接收器的编号，M为源的个数，N为每个源的接收器个数，Ei(m,rj,t)为第i个源对猜测模型正演计算的模拟数据，m为模型介质参数向量且m＝(μ(r),ε(r),σ(r))T，μ为磁导率，ε为介电常数，σ为电导率，为第i个源激发在rj处接收到的观测数据，rj为第j个接收器空间坐标向量；
(2)引入相对电导率，将模型参数m设定为相对介电常数和相对电导率的线性组合
实际的GPR全波形双参数反演过程中，通常为非磁介质，为了更准确地对模型参数定性，需要对介电常数、电导率同时反演。由于介电常数、电导率在数量级上相差很大，给反演计算带来了诸多不便。因此，如何设计一个能处理不同单位和敏感度的多参数反演策略，是GPR全波形双参数反演的关键。
考虑到相对介电常数可以根据真空介电常数来定义εr＝ε/ε0，因此，类似地可以引入相对电导率σr＝σ/σ0，取参考介质的电导率σ0＝1/η0，可以保证相对介电常数εr和相对电导率σr处于同一个量级，然后采用Lavoué等的做法，在反演过程中引入无量纲比例因子β，将模型参数m设定为相对介电常数和相对电导率的线性组合(εr,σr/β)，重写尺度变化之后的模型向量和梯度向量的明确表达式为：
上式中ε0与η0为常量，β为参数调整因子。这样，在优化过程中通过控制σr对εr的权重，避免由相对电导率与相对介电常数定义不准确引起反演过程的不稳定性；
(3)加入先验模型约束的正则化项
GPR全波形反演是一个高度非线性、不适定问题，对初始模型具有很高的依赖性，易陷入局部极小值及周波跳跃现象；为了全波形反演具有更好的稳定性、更高精度，本发明引入了全变差模型约束及多尺度反演策略；
反演不适定性最常用的解法为Tikhonov正则化方法，该方法通过加入先验模型约束的正则化项，使反问题更加稳定。但由于光滑性的缘故，Tikhonov正则化方法易导致目标区域与背景区域边界模糊。而全变差模型约束方法能有效改善Tikhonov正则化反演边界值的过度光滑，与背景区域区别更加明显，重构的目标体边缘轮廓更加清晰。引入了一个全变差正则化后，形成新的目标函数为：
Φ(m)＝Φd(m)+λΦm(m)
式中，m＝(εr,σr/β)T为尺度变换之后模型参数，Φd(m)为数据目标函数，λ为正则化因子，Φm(m)为模型参数目标函数，可表示为：
Φm(m)＝R(εr)+R(σr/β)
式中R为全变差正则化算子：
式中，Ω为模型空间，代表微分算子。全变差正则化算子的梯度可表示为：
目标函数Φ(m)的梯度可写为：
S7.对目标函数进行求解；具体为采用局部求导类的有限内存BFGS法求解；
全波形反演就是寻求目标函数的Φ(m)极小值的模型介质参数向量m；由于全波形反演计算量太大，为了减少计算量，本发明采用局部求导类的有限内存BFGS法(L-BFGS)求解。每次迭代需计算目标函数的Fréchet导数：
式中δm＝(δε(r),δσ(r),δμ(r))T，
对于每次迭代k：
m(k+1)＝m(k)+a(k)p(k)
式中，a(k)代表搜索步长，由一种不确定线性搜索法Wolfe准则决定；p(k)代表更新方向，可由双循环递归计算；L-BFGS法的每次迭代都需要计算目标函数的Fréchet导数：
δui是第i个源的场向量ui的Fréchet导数。场向量ui的Fréchet导数为：
δui(m,r,t)＝u′miδm
因为在t＝0时刻电磁波尚未传播，初始条件δui(m,r,0)＝0；
S8.引入伴随场；具体为引入伴随场
推导δu(m,r,t)与u(m,r,t)的关系，δu是δm引起的u的变化量，则有：
联立上两式可得：
式中δC和δD分别为C和D变分：
忽略高阶项上式可化为：
相应的δui与ui满足如下关系：
为了使目标函数的梯度可以显示表达，引入伴随场定义算子L*为算子L的伴随算子，详细的L*算子表达式参见附录A，根据伴随作用：
<L*w,δu>＝<w,Lδu>
式中<,>表示内积，根据内积的定义上式可以化为：
定义伴随场wi(m,r,t)方程满足如下微分方程和终止条件：
wi(m,r,t)＝0
式中iy是一个y方向的单位向量，δ(r-rj)为Dirac函数。将上式联立可得：
对上式左端的广义函数求体积积分得：
则有：
将上式右端项展开，并改写为内积形式有：
Φm′δm＝<gμ,δμ>+<gε,δε>+<gσ,δσ>
式中，gμ，gε和gσ代表目标函数对参数μ，ε和σ的梯度，分别为：
重写尺度变化之后的梯度向量为：
式中，ε0与η0为常量，β为参数调整因子；
全变差正则化算子的梯度可表示为：
目标函数Φ(m)的梯度写为：
S9.计算更新步长，并更新当前模型；具体为采用非精确线搜索方法Wolfe准则计算更新步长；
S10.以步骤S9得到的模型作为新的初始模型，重复步骤S2～S9直至反演结束，得到最终的反演结果。
以下结合一个实施例对本发明方法进行进一步说明：
(1)Overthrust模型构建及无噪声数据反演
如图2所示为本发明方法的模型1——Overthrust模型的参数分布示意图：“x”和“·”符号代表发射机和接收机的位置；图2(a)为模型1的相对介电常数分布，图2(c)为模型1的电导率分布；图2(b)为反演初始模型的介电常数分布，图2(d)为反演初始模型的电导率分布；
该模型显然更符合实际地质情况，模型的长度与深度分别为5m和10m，其中图2(a)的相对介电常数模型的变化范围设为3-30区间，图2(c)的电导率模型的变化范围在0-20mS/m区间。地表之上设置50cm的空气层，上层复杂分层区域代表地表相对干燥的砂层已回填覆盖区域；深度约3.5米的界面代表了地下水位。地下具有尖锐、变化范围较大的结构，模型中存在多个弯曲的褶皱及强衰减层，电导率约为10mS/m，这可能掩盖下面的结构。
以网格步长为0.05m对该模型进行离散，不含PML区域共计111×201个网格。限定反演中空气层中的参数，同时避免源和接收器位置处的奇异性，则待求的地下介电常数和电导率的参数为20301个。采用共偏移距的GPR数据采样方式，设置21个源，水平间隔0.5m，101个接收器，水平间隔0.1m，源和接收器均位于地面之上两个网格点-0.1m深处。反演的介电常数与电导率的初始模型设置如图2(b)与图2(d)所示，它们是图2(a)与图2(c)中的真实模型高斯平滑后的结果，仅仅描绘了介质参数变化的大致趋势。
如图3所示为本发明方法的模型1的100MHz中心频率的雷克子波的时域共剖点记录示意图：图3(a)为真实模型，图3(b)为初始模型，图3(c)为真实模型正演剖面图3(a)与初始模型正演剖面图3(b)两者的残差；
分析图3(a)可知，横坐标3.0-7.0m，双程走时20-40ns处的多道反射(绕射)波，推断由模型浅部多个褶皱的界面反射引起，而30-50ns处的同相轴不连续现象对应于模型中的断层。再观察图3(b)可知，该共激发点记录中图像相对干净，仅在0-25ns处观测到直达波，说明初始模型中由于缺乏对比度，无法观测到反射波。图3(c)中的多道双曲同相轴，主要由褶皱界面和构造引起的反射(绕射)雷达记录。
将多尺度串行双参数全波形反演应用于Overthrust模型，反演时采用60MHz与100MHz两个尺度，先考虑不加正则化项，当参数调节因子β分别取0.5和1.0。如图4所示为本发明方法的模型1的无噪数据反演结果示意图：图4(a)为使用比例因子β＝1.0反演无噪数据获得的介电常数模型(进行458次总迭代)；图4(b)为使用比例因子β＝0.5反演无噪数据获得的介电常数模型(进行833次总迭代)；图4(c)为使用比例因子β＝1.0反演无噪数据获得的电导率模型(进行458次总迭代)；图4(d)为使用比例因子β＝0.5反演无噪数据获得的电导率模型(进行833次总迭代)；
其中β＝0.5反演的归一化的最终数据拟合差为9.4394e-06，β＝1.0反演的归一化最终数据拟合差为3.3415e-05。由于参数因子β将会改变反演过程中电导率与介电常数的反演权重，因此，取不同的参数调节因子，得到的反演效果也将产生细微的差别。如图4(d)所示，当β＝0.5时，反演所得的电导率模型更平滑，幅值变化较小，仅能反映出浅部的分层情况。其图4(b)中的介电常数的反演剖面中，浅部褶皱界面得到了较好地恢复，浅部各层的厚度及深部的分界面细节也得到了较好反映，但是介电常数的反演数值较真实模型要小；如图4(c)所示，当β＝1.0时，反演所得的电导率模型各层界面更加分明，浅部低电导率的反演振荡更大，反演结果不稳定性更强，在一定程度上高估了电导率的变化情况。其图4(a)中的介电常数的反演剖面显示，浅部褶皱界面构造明显，细微处较真实值小，深部分界面模糊，介电常数反演数值较真实模型小。
如图5所示为本发明方法的相对介电常数一维深度剖面图：图5(a)为沿距离2.5m处的相对介电常数一维深度剖面图；图5(b)为沿距离5.0m处的相对介电常数一维深度剖面图；图5(c)为沿距离7.5m处的相对介电常数一维深度剖面图；
如图6所示为本发明方法的电导率的一维深度剖面图：图6(a)为沿距离2.5m处的电导率一维深度剖面图；图6(b)为沿距离5.0m处的电导率一维深度剖面图；图6(c)为沿距离7.5m处的电导率一维深度剖面图；
分析图5的相对介电常数反演曲线可知，两个β的结果在2.5m之上与真实模型基本一致，2.5-5m仅能体现界面的起伏情况，反演结果与真实模型存在一定偏差；总体而言β＝0.5的反演结果与真实结果更为接近。分析图6的电导率反演曲线可知，两个β的结果在总体上的变化趋势与真实结果基本一致，两者的值与真实结果均出现了错位及偏离现象；总体而言β＝0.5的反演结果比真实结果小，起伏平缓，β＝1.0的反演结果与真实结果变化程度基本相当，但是错位现象更加明显。
由于介电常数主要体现的是波动效应，对界面反映更加明显；电导率主要体现扩散效应，其体积效应更加明显。因此，无法从反演数据中恢复含有精确界面信息的电导率结果，可以以考虑以介电常数为主要参数，得到一个分辨率更高的地下构造重建结果。
显然取β＝0.5时的介电常数反演结果更为准确，而取β＝1.0时电导率成像构造更准确，是否能在[0.5，1]区间内取一个折中的参数调节因子，提高反演的精确度。可根据最终的数据拟合差作为判定标准，开展大量的β选取实验来求取更为精确的反演结果。但实际数据中多含有噪声，太小的数据拟合差可能导致对噪声的过拟合。
为了求取最优的反演效果，一个好的解决方案为：取参数调节因子β＝1.0并加载按等比例递减的正则化因子λ＝λ0qk，归一化的最终数据拟合差为1.4789e-05，得到高分率介电常数反演模型和电导率反演模型。如图7所示为本发明方法的模型1使用缩放因子β＝1.0和正则化权重λ＝5.0的无噪数据反演结果示意图：图7(a)为使用缩放因子β＝1.0和正则化权重λ＝5.0，通过无噪声数据反演获得的介电常数模型；图7(b)为沿距离2.5m，5m和7.5m处相对介电常数的一维深度剖面图，其中，实线表示真实模型，虚线表示反转模型，点画线表示初始模型(共770次迭代)；图7(c)为使用缩放因子β＝1.0和正则化权重λ＝5.0，通过无噪声数据反演获得的电导率模型；图7(d)为沿距离2.5m，5m和7.5m处电导率的一维深度剖面图，其中，实线表示真实模型，虚线表示反转模型，点画线表示初始模型(共770次迭代)；
对比分析图7(a)与图4(b)、图4(a)的介电常数反演结果，图7(a)反演剖面在深处的重建结果比图4(b)中β＝0.5反演结果要较差，但比图4(a)图中β＝1.0反演结果较好，这是由于全变差正则化的引入，提高了浅部界面的分辨率。而图7(b)的电导率的反演结果在深处的重建结果比图4(d)图中β＝0.5的反演结果要好，却比图4(c)中β＝1.0反演结果更为平缓。显然，这是一种折中方案，可使反演更加稳定，减少了界面内部的振荡，使界面更加清晰。
(2)噪声数据反演
为了验证该算法对噪声的适应性，将白噪声添加到正演数据中，加噪后数据的信噪比(SNR)为25dB。
如图8所示为本发明方法的SNR＝25dB的白噪声下模型1的时域共剖点记录示意图：图8(a)为真实模型，图8(b)为初始模型，图8(c)为残差；
如图9所示为本发明方法的模型1使用缩放因子β＝1.0和减小的正则化权重的噪声数据(SNR＝25dB)反演结果示意图：图9(a)为使用缩放因子β＝1.0和减小的正则化权重，通过噪声数据反演获得的介电常数模型；图9(b)为沿距离2.5m，5m和7.5m处相对介电常数的一维深度剖面图，其中，黑色曲线表示真实模型，红色曲线表示反转模型，蓝色曲线表示初始模型(共202次迭代)；图9(c)为使用缩放因子β＝1.0和减小的正则化权重，通过噪声数据反演获得的电导率模型；图9(d)为沿距离2.5m，5m和7.5m处电导率的一维深度剖面图，其中，黑色曲线表示真实模型，红色曲线表示反转模型，蓝色曲线表示初始模型(共202次迭代)；
由图8可知，噪声的加入使雷达剖面数据更为杂乱，加大了反演难度。在图9(a)介电常数图像中，该反演曲线没有对噪声信息进行拟合，得到了较为准确的浅部构造信息，尽管浅部层中含有一些振荡，但浅部的分层仍然能准确地重建，分辨率随着深度显着降低，深部的界面无法识别，并且高介电常数层不明显出现。图9(b)为电导率反演剖面图，图中可见，地下结构的主要趋势在电阻率部面图中也得到较好的反映，说明该反演方法具有一定的抗噪性，但反演的电导率深度曲线与真实模型存在偏差，剖面图像也稍显模糊。显然，含噪数据的重建结果达不到无噪数据的反演精度。
本发明方法的优势如下：
GPR时间域全波形反演计算量巨大，内存要求高，在普通微机上计算难度大。为进一步优化反演算法，本发明应用多尺度串行反演及全变差模型约束策略，将反演问题分解为2-3个从低频至高频的逐频带多尺度反演，从而避免反演陷入局部极值，有利于搜索全局极值点；并通参数调节因子β及正则化参数λ的选取实验，说明合适的参数调节因子及正则化参数，能有效保证双参数反演的收敛速度与稳定性，提高雷达剖面的反演精度。
①多尺度串行反演策略优势
为了分析多尺度串行策略及不同参数对反演结果影响，设置图10所示的模型2。图10为该“中”字型模型的参数分布示意图：包括两个交叉十字的异常体，“x”和“·”符号代表发射机和接收机的位置；图10(a)为相对介电常数分布，图10(b)为电导率分布；
模型背景为均匀介质，其相对介电常数和电导率分别为εr＝4和σ＝3mS/m，在该均匀介质中包含了两个形状相同的“中”字型异常体。左边异常体1的外沿轮廓的相对介电常数为8，电导率为10ms/m，“中”字型内部空腔区域相对介电常数为1，电导率为0ms/m，而右边异常体2的介电常数与电性参数与异常体1内外取相反的参数。模拟区域的边界设有40个“x”型点源和200个“·”表示的接收器，激励源为200MHz的Ricker子波，时窗长度为75ns，网格离散步长为0.05m，时间步长为0.1ns。反演中采用均匀介质作为初始模型。
为了计算目标函数中的残差及多尺度串行反演中的高低频雷达数据，对图10中预设模型开展正演。正演时将Ricer子波激励源置于(2.5m，0.0m)，接收天线布设200道，其中1-50道、51-100道、101-150道、151-200道的接收天线分别为上边界z＝0m、右边界x＝5m、下边界z＝5m、左边界x＝0m处，均匀等间距自左至右接收雷达信号。
如图11所示为本发明方法的模型2的200MHz中心频率的雷克子波的时域共剖点记录示意图：图11(a)为真实模型，图11(b)为初始均匀模型，图11(c)为残差；
如图12所示为本发明方法的Wiener滤波操作示意图：图12(a)中，点划线是200MHz的原始Ricker子波，实线是80MHz的目标Ricker子波，点划线是Wiener滤波所得的波形，点划线与实线较好地重合；图12(b)为Wiener滤波后得到的800MHz时域共炮点记录；
首先分析多尺度串行策略的影响，以模型2为例进行单尺度和多尺度反演计算。反演过程中，为了使分析更纯粹，令β＝1.0，λ＝0(未加模型约束)，仅仅探讨单尺度与多尺度对反演效果的影响。单尺度和多尺度反演都在个人计算机上GPU并行运算，反演终止条件有两个，一个是相对目标函数Φk/Φ0的比值小于1e-5，二是步长为零。
如图13所示为本发明方法的三种不同反演策略的结果示意图：图13(a)为应用原始数据进行单尺度反演得到的相对介电常数结果；图13(b)为应用滤波数据低频带进行单尺度反演得到的相对介电常数结果；图13(c)为进行多尺度反演得到的相对介电常数结果；图13(d)为应用原始数据进行单尺度反演得到的电导率结果；图13(e)为应用滤波数据低频带进行单尺度反演得到的电导率结果；图13(f)为进行多尺度反演得到的电导率结果；
多尺度反演时，使用多个频带的源和数据，分别采用Wiener低通滤波后的80MHz起始频率和200MHz的原始频带。数据首先在低频带反演20次；然后再进行高频带反演，直至达到反演停止标准。在图13(a)的单尺度介电常数反演剖面中，右下方异常体2形态重构较准确，由于该异常体介电常数较低，反演效果较好，而左上方异常体1畸变非常严重，高介电常数甚至被反演成低介电常数。而图13(d)的电导率反演剖面的两个异常体的形态非常模糊，异常体的形态存在较大扭曲，大小、空间位置都存在较大偏差。而多尺度的低频介电常数反演剖面图13(b)与图13(e)电导率反演剖面都能刻画出“中”字型异常体的大致轮廓，但异常体形态出现了一定程度的畸变，相较图13(b)的介电常数反演剖面，图13(e)中的电导率反演剖面畸变更为严重，“中”字形出现了扭曲、变形，右下方的“中”甚至难以分辨出该异常体的原来轮廓，重建效果并不理想。图13(c)与图13(f)的多尺度串行反演的介电常数与电导率剖面，都能很好地重构出两个“中”字型异常体，异常体的形态、大小、空间位置都能准确反映，反演效果最好。对比单尺度反演、多尺度仅低频反演和多尺度综合反演可以发现：双参数反演时，介电常数的反演结果较电导率的反演效果更好，而多尺度串行反演由于充分结合了高、低频的信息，对异常体的大小、形态、空间位置重构更准确，分界面更加明显，异常体的轮廓更清晰可辨，反演出的背景及异常体的电性参数与真实模型基本一致。
为了更精确地描述几种反演方式在反演精度方面的异同，在模型图10中选取两个截面，具体可见图10(a)与图10(b)中两条虚线AA’与BB’，截面位置分别位于深度z＝1.6m、z＝3.25m处。
如图14所示为本发明方法的单道反演切线示意图：图14(a)为z＝1.6m的介电常数截面，图14(b)为z＝3.25m的介电常数截面，图14(c)为z＝1.6m的电导率截面，图14(d)为z＝3.25m的电导率截面差；
图14(a)中可见，虚线表示的单尺度介电常数全波形反演曲线与实线表示的真实模型幅值存在较大的偏差，曲线更尖细，表明介电常数重构不精确，同时第1个与第3个尖峰出现位置与真实模型存在偏差，说明对异常体形态归位亦不准确，而点划线表示的Wiener低通滤波后的仅80MHz的多尺度反演曲线，该曲线波峰位置与实线表示的真实模型能较好地吻合，说明异常体位置反演比较准确，但幅值存在一定偏差，说明介电常数重构结果仍有待提高。而点划线表示多尺度综合反演曲线与实线表示的真实模型曲线能较好地拟合，无论是幅值大小、曲线拐点位置都能较好地对应，说明多尺度综合反演具有较高的精度。而图14(c)为与图14(a)相对应的电导率反演结果，总体来看，电导率的反演精度较介电常数的反演精度略差。图14(b)为z＝3.25m处不同反演方式的曲线对比，图中可见，虚线表示的单尺度及点划线表示的多尺度串行反演与实线表示的真实模型拟合较好，而点划线的幅值与实线偏离较大。图14(d)为对应的电导率反演结果，虚线和点划线与实线存在较大的偏差，说明电导率的反演精度不高，而点划线和实线位置及幅值都能大致拟合，说明多尺度电导率串行反演的精度能够得到有效保证。由此可见，两个不同深度截取的单道反演曲线对比表明，相对于单尺度反演，多尺度串行策略的反演结果与真实模型更为接近，特别是介电常数与真实模型基本吻合，电导率反演结果趋势一致，界面稍显模糊。
如图15所示为本发明方法的不同反演策略下目标函数及模型重构误差对比示意图：图15(a)为归一化的目标函数收敛曲线图，实线表示多尺度串行反演，虚线表示单尺度反演；图15(b)为模型重构误差曲线，点划线表示单尺度相对介电常数重构误差曲线，点线表示单尺度电导率重构误差曲线，实线表示多尺度相对介电常数重构误差曲线，虚线表示多尺度电导率重构误差曲线；
由图15可知，多尺度串行的目标函数比单尺度目标函数低一个数量级，进一步说明了多尺度串行反演结果具有更高的精度。结合4条误差曲线可知，相比单尺度反演曲线，多尺度串行反演结果与真实模型更为接近，可以收敛到全局最小，具有比单尺度波形反演成像更快的收敛速度，拟合效果更好。而且，相比电导率重构误差，介电常数重构误差更少，结果与真实模型更为接近。
②全变差模型约束策略优势
a.参数调节因子选取
双参数反演过程中，由于介电常数与电导率在数值上相差较大，同步反演时需要引入参数调节因子β，合理的β值可以保证反演快速稳定的收敛到全局极小值，因此，最优的β值选取非常重要。
如图16所示为本发明方法的不同参数调节因子下目标函数及模型重构误差对比示意图：图16(a)为不同参数调节因子下目标函数收敛曲线图；图16(b)为模型重构误差曲线，实线表示相对介电常数重构误差曲线，虚线表示电导率重构误差曲线；
图中可见，当β＝1.0时，最终的数据目标函数以及介电常数和电导率的模型重构误差达到最小。而图16(a)中的最终数据目标函数随β因子变化趋势与图16(b)中的重构误差趋势基本相同。
为了进一步探讨参数调节因子β对反演效果的影响，在不加载模型约束情况下(β＝0)采用相同多尺度策略及终止条件，仍以图10中模型为例，分别取β＝0.5，1.0，2.0，进行全波形双参数反演。
如图17所示为本发明方法的不同参数调节因子下反演结果示意图：图17(a)为β＝0.5的相对介电常数反演结果；图17(b)为β＝1.0的相对介电常数反演结果；图17(c)为β＝2.0的相对介电常数反演结果；图17(d)为β＝0.5的电导率反演结果；图17(e)为β＝1.0的电导率反演结果；图17(f)为β＝2.0的电导率反演结果；
分析图17(a)、图17(b)、图17(c)可知，不同的β取值对介电常数的反演结果仅有细微的影响。观察电导率反演剖面图17(d)、图17(e)、图17(f)，随着β取值的增大，电导率参数在反演中所占权重会逐步增大，当β取值过大时，电导率反演结果会发生振荡，而当β＝1.0时，图17(e)的反演效果最佳，该反演结果与图16中得出的结论一致，说明前文公式中相对电导率选取的σ0对σr＝σ/σ0进行归一化是较为合适的选择。
b.正则化参数选取
全变差正则化方法求解雷达全波形反演，可以使反演更稳定。由于正则化参数控制着模型范数Φm和数据拟合残差范数Φd在目标函数中所占的比例，它的选取对解的性态起着关键的作用。如果正则化参数λ取的太小，则意味着模型范数在代价函数中所占的比例少，对反演结果平滑程度的影响也小，此时噪声得不到很好的抑制，解的不稳定性仍然存在；如果λ选取的太大，噪声将会得到很好的抑制，稳定性得到了保证，却使反演结果过于平滑，细节不易突出，加大了原问题与新问题之间的偏离程度，所求的根本不是原问题的解了，因此，一个好的正则化参数能合理的平衡两者之间的关系。
常用的λ值取值方案有两种：(a)固定正则化参数，譬如取λ＝1.0，(b)设定一个渐变的正则化参数，如：设λ＝λ0qk，λ0为初始正则化因子，q为等比因子，k为迭代次数，本发明取λ0＝1.0，q＝0.97，使正则化因子等比例逐步递减。考虑到GPR反演开始于给定的初始模型，此时，目标函数中模型范数Φm较小，甚至可能为0，而模型数据的拟合残差范数Φd一般都很大，随着反演迭代的进行，Φm会逐渐增加而Φd逐渐下降，在这个过程中正则化因子λ需要平衡目标函数中的两部分Φd和Φm，使两者不至于相差过大。因此，选取第2种渐变的正则化参数更为合理。
为了说明正则化方法的效果及正则化参数选取对反演结果的影响，仍以“中”字型模型为例，设置3组双参数反演实验，(a)参数调节因子β＝1.0，正则化因子λ＝1.0；(b)参数调节因子β＝1.0，正则化因子按等比例递减λ＝λ0qk；(c)参数调节因子β＝2.0，正则化因子按等比例递减λ＝λ0qk。
如图18所示为本发明方法的不同正则化参数下反演结果示意图：图18(a)为β＝1.0，λ＝1.0的相对介电常数反演结果；图18(b)为β＝1.0，λ等比递减的相对介电常数反演结果；图18(c)为β＝2.0，λ等比递减的相对介电常数反演结果；图18(d)为β＝1.0，λ＝1.0的电导率反演结果；图18(e)为β＝1.0，λ等比递减的电导率反演结果；图18(f)为β＝2.0，λ等比递减的电导率反演结果；
综合对比6幅反演剖面图可知，参数设置为条件(b)时的图18(b)与图18(e)反演图像最清晰，异常体重构效果最好；参数设置为条件(c)时的图18(c)与图18(f)反演效果其次；参数设置为条件(a)时的图18(a)与图18(d)反演效果最差，尤其是图18(d)的电导率反演剖面图像失真较厉害，“中”字型出现了严重畸变，异常体轮廓也不清晰。再将图18的反演结果与图16和图17的反演结果进行对比，可知：若不加载正则化项，由于目标函数仅强调数据的拟合，导致目标函数的非线性较强，当尺度因子选取不合理时，拟合数据的模型可能包含较多高波数虚假构造，如图16(f)所示，如果数据中含有噪声可能对数据拟合过度而拟合了其中的误差部分；若正则化因子取得过大，则过分强调模型的光滑度而忽略数据的拟合，使得反演提早终止迭代，反演所得的重构模型较简单、光滑，但对数据的拟合度很低，如图18(d)所示。而加载一个等比例递减的正则化因子λ比固定正则化因子反演效果更好，它可以使数据拟合与模型约束之间达到更好的平衡，使反演更加快速稳定收敛。同时，正则化因子的引入，可以降低反演中尺度因子β的敏感度，扩大最佳尺度因子β的选择区间，在一定程度上压制背景的非物理振荡，提高反演剖面的分辨率。
如图19所示为本发明方法的不同参数调节因子及正则化参数下目标函数及模型重构误差对比示意图：图19(a)为不同参数调节因子及正则化参数下的目标函数收敛曲线图，虚线为β＝1.0、λ＝1.0下的反演结果，实线为β＝1.0、λ等比递减下的反演结果，点划线为β＝2.0、λ等比递减下的反演结果；图19(b)为不同参数调节因子及正则化参数下的模型重构误差曲线，带“x”号标记虚线为β＝1.0、λ＝1.0下的相对介电常数重构误差曲线，带“x”号标记实线为β＝1.0、λ等比递减下的相对介电常数重构误差曲线，带“x”号标记点划线为β＝2.0、λ等比递减下的相对介电常数重构误差曲线，虚线为β＝1.0、λ＝1.0下的电导率重构误差曲线，实线为β＝1.0、λ等比递减下的电导率重构误差曲线，点划线为β＝2.0、λ等比递减下的电导率重构误差曲线；
由图19可知，实线表示的参数调节因子β＝1.0，正则化因子按等比例递减λ＝λ0qk条件下的反演目标函数最小，收敛最好。图19(b)中，带“x”号标记线图表示的介电常数重构误差曲线较不带标记线图表示的电导率重构误差曲线整体要小。

Claims (7)

1.一种探地雷达多尺度全波形双参数反演方法，包括如下步骤：
S1.将观测数据的GPR时间域共剖点记录作为输入，加载子波，建立初始模型并设定迭代精度；
S2.对建立的初始模型进行正演；
S3.确定多频率反演顺序；
S4.构造Wiener滤波器；
S5.采用步骤S4构造的Wiener滤波器进行滤波后，得到起始频带和原始频带；数据先在低频带进行反演，然后再在高频带进行反演；
S6.构建反演目标函数；
S7.对目标函数进行求解；
S8.引入伴随场；
S9.计算更新步长，并更新当前模型；
S10.以步骤S9得到的模型作为新的初始模型，重复步骤S2～S9直至反演结束，得到最终的反演结果。
2.根据权利要求1所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S2所述的对建立的初始模型进行正演，具体为采用基于交错网格的时域有限差分法进行GPR正演，同时时间离散精度和空间离散精度均设为2阶，截断边界处采用非分裂的卷积完全匹配层。
3.根据权利要求1所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S4所述的构造Wiener滤波器，具体为采用如下算式作为Wiener滤波器：
式中fwiener(ω)为频率域Wiener低通滤波器，Wtarget(ω)为目标激励源子波频谱，*为共轭符号，Woriginal(ω)为原始激励源子波频谱，δ为防止分母为零的常量小数。
4.根据权利要求1～3之一所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S6所述的构建反演目标函数，具体为采用如下步骤构建反演目标函数：
(1)构建数据目标函数Φd(m)：
式中i为源的编号，j为接收器的编号，M为源的个数，N为每个源的接收器个数，Ei(m,rj,t)为第i个源对猜测模型正演计算的模拟数据，m为模型介质参数向量且m＝(μ(r),ε(r),σ(r))T，μ为磁导率，ε为介电常数，σ为电导率，为第i个源激发在rj处接收到的观测数据，rj为第j个接收器空间坐标向量；
(2)引入相对电导率σr＝σ/σ0，取参考介质的电导率σ0＝1/η0；引入无量纲比例因子β，将模型参数m设定为相对介电常数和相对电导率的线性组合(εr,σr/β)；
(3)加入先验模型约束的正则化项；在引入一个全变差正则化后，形成新的目标函数Φ(m)：
Φ(m)＝Φd(m)+λΦm(m)
式中Φd(m)为数据目标函数，λ为正则化因子，Φm(m)为模型参数目标函数且Φm(m)＝R(εr)+R(σr/β)，Ω为模型空间，为微分算子，R为全变差正则化算子且
5.根据权利要求4所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S7所述的对目标函数进行求解，具体为采用局部求导类的有限内存BFGS法求解，每次迭代需计算目标函数的Fréchet导数：
式中δm＝(δε(r),δσ(r),δμ(r))T，
6.根据权利要求5所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S8所述的引入伴随场，具体为引入伴随场定义算子L*为算子L的伴随算子，根据伴随作用：
<L*w,δu>＝<w,Lδu>
式中<,>表示内积；
目标函数对参数μ，ε和σ的梯度分别为gμ，gε和gσ：
重写尺度变化之后的梯度向量为：
式中，ε0与η0为常量，β为参数调整因子；
全变差正则化算子的梯度可表示为：
目标函数Φ(m)的梯度写为：
7.根据权利要求6所述的探地雷达多尺度全波形双参数反演方法，其特征在于步骤S9所述的计算更新步长，具体为采用非精确线搜索方法Wolfe准则计算更新步长。
CN201910477198.5A 2019-06-03 2019-06-03 探地雷达多尺度全波形双参数反演方法 Active CN110095773B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910477198.5A CN110095773B (zh) 2019-06-03 2019-06-03 探地雷达多尺度全波形双参数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910477198.5A CN110095773B (zh) 2019-06-03 2019-06-03 探地雷达多尺度全波形双参数反演方法

Publications (2)

Publication Number Publication Date
CN110095773A true CN110095773A (zh) 2019-08-06
CN110095773B CN110095773B (zh) 2022-11-01

Family

ID=67449991

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910477198.5A Active CN110095773B (zh) 2019-06-03 2019-06-03 探地雷达多尺度全波形双参数反演方法

Country Status (1)

Country Link
CN (1) CN110095773B (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111324972A (zh) * 2020-03-16 2020-06-23 郑州大学 一种基于gpu并行的探地雷达电磁波数值模拟计算方法
CN111830905A (zh) * 2020-08-10 2020-10-27 哈尔滨工业大学 一种基于简化牛顿法的多维系统轮廓误差估计方法
CN112434439A (zh) * 2020-12-01 2021-03-02 中国自然资源航空物探遥感中心 约束物性多区间分布的多元段反演方法及装置
CN112731379A (zh) * 2020-12-15 2021-04-30 郑州大学 粒子群遗传联合介电常数获取方法、雷达检测方法和系统
CN112748466A (zh) * 2019-10-30 2021-05-04 中国石油天然气集团有限公司 一种基于菲涅尔体的旅行时场数据处理方法及装置
CN112882022A (zh) * 2021-01-18 2021-06-01 云南航天工程物探检测股份有限公司 一种探地雷达数据时频组合全波形反演方法
CN113240791A (zh) * 2021-04-25 2021-08-10 北京航空航天大学 一种地下目标高精度成像探测方法
CN113569493A (zh) * 2021-09-26 2021-10-29 自然资源部第一海洋研究所 一种基于域分解的物理驱动深度学习反演方法
CN114545863A (zh) * 2022-03-07 2022-05-27 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法
CN116295158A (zh) * 2023-05-17 2023-06-23 四川华恒升科技发展有限公司 一种基于机载双频测深仪测量淤泥深度的仪器及方法
CN116327250A (zh) * 2023-02-13 2023-06-27 中国科学院地质与地球物理研究所 一种基于全波形反演技术的乳腺超声三维成像方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110238390A1 (en) * 2010-03-29 2011-09-29 Krebs Jerome R Full Wavefield Inversion Using Time Varying Filters
US20120314538A1 (en) * 2011-06-08 2012-12-13 Chevron U.S.A. Inc. System and method for seismic data inversion
CN106908835A (zh) * 2017-03-01 2017-06-30 吉林大学 带限格林函数滤波多尺度全波形反演方法
CN107450102A (zh) * 2017-07-28 2017-12-08 西安交通大学 基于分辨率可控包络生成算子的多尺度全波形反演方法
CN108254731A (zh) * 2018-04-25 2018-07-06 吉林大学 探地雷达数据的多尺度阶梯式层剥离全波形反演方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
US20180356548A1 (en) * 2017-06-12 2018-12-13 Institute Of Geology And Geophysics Chinese Academy Of Sciences Inversion velocity model, method for establishing the same and method for acquiring images of underground structure

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110238390A1 (en) * 2010-03-29 2011-09-29 Krebs Jerome R Full Wavefield Inversion Using Time Varying Filters
US20120314538A1 (en) * 2011-06-08 2012-12-13 Chevron U.S.A. Inc. System and method for seismic data inversion
CN106908835A (zh) * 2017-03-01 2017-06-30 吉林大学 带限格林函数滤波多尺度全波形反演方法
US20180356548A1 (en) * 2017-06-12 2018-12-13 Institute Of Geology And Geophysics Chinese Academy Of Sciences Inversion velocity model, method for establishing the same and method for acquiring images of underground structure
CN107450102A (zh) * 2017-07-28 2017-12-08 西安交通大学 基于分辨率可控包络生成算子的多尺度全波形反演方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
CN108254731A (zh) * 2018-04-25 2018-07-06 吉林大学 探地雷达数据的多尺度阶梯式层剥离全波形反演方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
冯德山等: "基于GPU并行的时间域全波形优化共轭梯度法快速GPR双参数反演", 《地球物理学报》 *
廖文等: "基于改进目标函数的频率域全波形反演研究", 《科技通报》 *
成景旺等: "频率域反射波全波形速度反演", 《地球科学(中国地质大学学报)》 *
李伦等: "基于正则化方法的高频地波雷达海浪方向谱反演", 《地球物理学报》 *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112748466A (zh) * 2019-10-30 2021-05-04 中国石油天然气集团有限公司 一种基于菲涅尔体的旅行时场数据处理方法及装置
CN112748466B (zh) * 2019-10-30 2024-03-26 中国石油天然气集团有限公司 一种基于菲涅尔体的旅行时场数据处理方法及装置
CN111324972A (zh) * 2020-03-16 2020-06-23 郑州大学 一种基于gpu并行的探地雷达电磁波数值模拟计算方法
CN111324972B (zh) * 2020-03-16 2023-02-24 郑州大学 一种基于gpu并行的探地雷达电磁波数值模拟计算方法
CN111830905A (zh) * 2020-08-10 2020-10-27 哈尔滨工业大学 一种基于简化牛顿法的多维系统轮廓误差估计方法
CN112434439A (zh) * 2020-12-01 2021-03-02 中国自然资源航空物探遥感中心 约束物性多区间分布的多元段反演方法及装置
CN112434439B (zh) * 2020-12-01 2024-02-27 中国自然资源航空物探遥感中心 约束物性多区间分布的多元段反演方法及装置
CN112731379A (zh) * 2020-12-15 2021-04-30 郑州大学 粒子群遗传联合介电常数获取方法、雷达检测方法和系统
CN112882022B (zh) * 2021-01-18 2023-08-11 云南航天工程物探检测股份有限公司 一种探地雷达数据时频组合全波形反演方法
CN112882022A (zh) * 2021-01-18 2021-06-01 云南航天工程物探检测股份有限公司 一种探地雷达数据时频组合全波形反演方法
CN113240791A (zh) * 2021-04-25 2021-08-10 北京航空航天大学 一种地下目标高精度成像探测方法
CN113569493A (zh) * 2021-09-26 2021-10-29 自然资源部第一海洋研究所 一种基于域分解的物理驱动深度学习反演方法
CN114545863A (zh) * 2022-03-07 2022-05-27 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法
CN114545863B (zh) * 2022-03-07 2024-02-13 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法
CN116327250B (zh) * 2023-02-13 2023-08-25 中国科学院地质与地球物理研究所 一种基于全波形反演技术的乳腺超声三维成像方法
CN116327250A (zh) * 2023-02-13 2023-06-27 中国科学院地质与地球物理研究所 一种基于全波形反演技术的乳腺超声三维成像方法
CN116295158B (zh) * 2023-05-17 2023-08-08 四川华恒升科技发展有限公司 一种基于机载双频测深仪测量淤泥深度的仪器及方法
CN116295158A (zh) * 2023-05-17 2023-06-23 四川华恒升科技发展有限公司 一种基于机载双频测深仪测量淤泥深度的仪器及方法

Also Published As

Publication number Publication date
CN110095773B (zh) 2022-11-01

Similar Documents

Publication Publication Date Title
CN110095773A (zh) 探地雷达多尺度全波形双参数反演方法
Müller-Petke et al. MRSmatlab—A software tool for processing, modeling, and inversion of magnetic resonance sounding data
CN113221393B (zh) 一种基于非结构有限元法的三维大地电磁各向异性反演方法
US11181653B2 (en) Reservoir characterization utilizing ReSampled seismic data
US8731987B2 (en) Method and apparatus to automatically recover well geometry from low frequency electromagnetic signal measurements
CN105388518A (zh) 一种质心频率与频谱比联合的井中地震品质因子反演方法
CN108267722A (zh) 地质雷达回波信号物性解构与探测目标数字重构智能化识取方法
CN103163554A (zh) 利用零偏vsp资料估计速度和q值的自适应波形的反演方法
CN103792573A (zh) 一种基于频谱融合的地震波阻抗反演方法
CN109061727A (zh) 一种频率域粘声介质全波形反演方法
CN106707334A (zh) 一种提高地震资料分辨率的方法
CN109459789A (zh) 基于振幅衰减与线性插值的时间域全波形反演方法
CN104730576A (zh) 基于Curvelet变换的地震信号去噪方法
US9702995B2 (en) Domain freezing in joint inversion
CN106842297A (zh) 井约束非稳态相位校正方法
CN115906559B (zh) 一种基于混合网格的大地电磁自适应有限元正演方法
Jin et al. A robust learning method for low-frequency extrapolation in gpr full waveform inversion
CN108363739A (zh) 一种基于稀疏采集的地震资料高低频拓展方法
Zhong et al. Electrical resistivity tomography with smooth sparse regularization
CN110865409B (zh) 一种基于波阻抗低秩正则化的地震波阻抗反演方法
CN110244383B (zh) 基于近地表数据的地质岩性综合模型创建方法
Xia et al. Automatic demarcation of sequence stratigraphy using the method of well logging multiscale data fusion
Su et al. Waveform and resistivity data fusion imaging method based on the reflection coefficient
CN104375171B (zh) 一种高分辨率地震反演方法
Hu et al. Fast inversion of array laterolog measurements in an axisymmetric medium

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20221110

Address after: Yuelu District City, Hunan province 410083 Changsha Lushan Road No. 932

Patentee after: CENTRAL SOUTH University

Address before: Yuelu District City, Hunan province 410083 Changsha Lushan Road No. 932

Patentee before: CENTRAL SOUTH University

Patentee before: ZHEJIANG HUADONG ENGINEERING SAFETY TECHNOLOGY CO.,LTD.