CN113970732A - 一种三维频率域探地雷达双参数同步反演方法 - Google Patents

一种三维频率域探地雷达双参数同步反演方法 Download PDF

Info

Publication number
CN113970732A
CN113970732A CN202111243358.3A CN202111243358A CN113970732A CN 113970732 A CN113970732 A CN 113970732A CN 202111243358 A CN202111243358 A CN 202111243358A CN 113970732 A CN113970732 A CN 113970732A
Authority
CN
China
Prior art keywords
representing
model
current
iteration
dielectric constant
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.)
Pending
Application number
CN202111243358.3A
Other languages
English (en)
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.)
Third Engineering Co Ltd Of Cccc Third Highway Engineering Co ltd
Central South University
Guangzhou Municipal Engineering Design & Research Institute Co Ltd
Original Assignee
Third Engineering Co Ltd Of Cccc Third Highway Engineering Co ltd
Central South University
Guangzhou Municipal Engineering Design & Research Institute Co Ltd
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 Third Engineering Co Ltd Of Cccc Third Highway Engineering Co ltd, Central South University, Guangzhou Municipal Engineering Design & Research Institute Co Ltd filed Critical Third Engineering Co Ltd Of Cccc Third Highway Engineering Co ltd
Priority to CN202111243358.3A priority Critical patent/CN113970732A/zh
Publication of CN113970732A publication Critical patent/CN113970732A/zh
Pending legal-status Critical Current

Links

Images

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
    • 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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种三维频率域探地雷达双参数同步反演方法，包括输入观测数据，设置三维初始模型，并初始化参数；对建立的三维初始模型正演；构建反演目标函数，并分解为两个交替的子问题；初始化介电常数和电导率，得到初始先导模型向量；根据先导模型向量求解子问题1，并更新模型介质参数向量；根据得到的模型介质参数向量求解子问题2，并更新先导模型向量；同步交替进行更新，直到满足反演终止条件，输出反演结果。本发明能够同步更新介电常数和电导率值，以提高反演的效率；改进全变差模型约束方法到全波形反演过程中，增加算法的鲁棒性，结合频率权重因子，采用并行方式实现了3D的GPR定量成像，使成像过程更加高效精确。

Description

一种三维频率域探地雷达双参数同步反演方法
技术领域
本发明属于地球物理探测领域，具体涉及一种三维频率域探地雷达双参数同步反演方法。
背景技术
探地雷达(Ground penetrating radar,GPR)是一种对地下异常体、结构和特性或物体内部不可见目标体进行定位的电磁无损探测技术。实际GPR探测中仅仅依靠获取的雷达剖面，即可大致推断探测目标体的位置、大小与尺寸。但是，探底雷达多局限于粗略估计，尚不能给出异常体的介电常数、电导率等参数，实现目标体属性的精准界定。全波形反演(Full waveform inversion,FWI)是新兴起的重要成像手段。GPR全波形反演能够获取介质的介电常数、电导率等物性参数，有助于岩性的精确描述和异常体预测，是重建近地表复杂地质结构的有力工具。
目前FWI已在GPR成像中得到广泛关注与研究，可分为时间域FWI及频率域FWI。在频率域中，通过将数据从低频率反转到高频率，可以在计算期间有效地减少问题的非线性。通过合理地选取少数频率分量数据进行反演，即可达到与时间域反演一致的效果。
鉴于电磁波横向波的3D散射比2D纵波的散射更复杂，模拟数据与模型参数之间非线性更强，而且3D模型中会有更多的地下参数拟合观测数据(多解性更强)。再加上解决3D反演问题需要更大的存储成本，并且其涉及到的正演问题多次计算成本远高于2D，因此现有技术多是关于2D探地雷达数据FWI的研究。然而对于3D变化的物体，2D成像方法限制了FWI的适用性和准确性。
考虑到探地雷达全波形反演是一个高度非线性、不适定问题。Lavoué等在2D频率域GPR全波形反演中引入了Tikhonov正则化，改善了电导率模型中的高振幅振荡现象，然而Tikhonov正则化的反演倾向于产生平滑模型。Watson使用近似的全变差正则化约束介电常数模型，使其反演结果具有尖锐界面，这种近似全变差方法可能不稳定并且反演的收敛性不能保证。因此需要寻找新的方法来进行GPR数据全波形反演。
发明内容
本发明的目的在于提供一种三维频率域探地雷达双参数同步反演方法，该方法能够提高探地雷达反演精度。
本发明提供的这种三维频率域探地雷达双参数同步反演方法，包括如下步骤：
S1.输入观测数据，设置三维初始模型，并初始化参数；
S2.对建立的三维初始模型正演；
S3.构建反演目标函数，并分解为两个交替的子问题，包括子问题1和子问题2，子问题1为求解模型介质参数向量，子问题2为求解先导模型向量；
S4.初始化介电常数和电导率，得到初始先导模型向量；
S5.根据先导模型向量求解子问题1，并更新模型介质参数向量；
S6.根据得到的模型介质参数向量求解子问题2，并更新先导模型向量；
S7.重复步骤S5-S6，同步交替进行更新，直到满足反演终止条件，结束反演并输出反演结果。
所述的步骤S2包括，利用有限单元法对建立的三维初始模型进行正演，具体为采用线性的矢量基函数的有限单元法进行模拟区域的离散，剖分网格为矩形块单元，截断边界处采用各向异性完全匹配层；三维GPR满足的频率域麦克斯韦方程表示为：
Figure BDA0003319979520000021
其中，
Figure BDA0003319979520000022
表示电场强度；ε表示当前介电常数；μ表示磁导率；σ表示当前电导率；ω表示角频率，ω＝2πf，f表示频率；j表示虚数单位，
Figure BDA0003319979520000023
为电流密度。
求得多源***的线性方程组：
Figure BDA0003319979520000024
其中，
Figure BDA0003319979520000025
表示与频率相关的所有源公用的系数矩阵；
Figure BDA0003319979520000026
表示离散电场值；
Figure BDA0003319979520000027
表示离散化激励源项；ω表示角频率；ε表示介电常数；σ表示电导率；s表示当前位置的激励源；与频率相关的所有源公用的系数矩阵
Figure BDA0003319979520000031
由所有单元矩阵
Figure BDA0003319979520000032
合成：
Figure BDA0003319979520000033
其中，
Figure BDA0003319979520000034
表示当前单元e内的体积积分；
Figure BDA0003319979520000035
和
Figure BDA0003319979520000036
为矢量基函数；使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；对矩阵
Figure BDA0003319979520000037
进行一次LU分解；对于单源单个频率的正演算子描述为：
Figure BDA0003319979520000038
其中，
Figure BDA0003319979520000039
表示角频率ω和当前激励源s下的模拟数据；
Figure BDA00033199795200000310
表示当前源的接收点的对应位置矩阵；
Figure BDA00033199795200000311
表示离散电场值；
Figure BDA00033199795200000312
表示与频率相关的所有源公用的系数矩阵；
Figure BDA00033199795200000313
表示离散化激励源项。
所述的步骤S3，利用已知的实测数据重构地下介质中的模型参数，包括介电常数ε和电导率σ的空间分布；并引入全变差正则化方法，反演目标函数
Figure BDA00033199795200000314
定义为：
Figure BDA00033199795200000315
其中，
Figure BDA00033199795200000316
表示模型介质参数向量，
Figure BDA00033199795200000317
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000318
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200000319
表示先验模型向量，
Figure BDA00033199795200000320
表示当前介电常数对应的先验模型向量，
Figure BDA00033199795200000321
表示当前电导率对应的先验模型向量；
Figure BDA00033199795200000322
表示数据目标函数；ε表示当前介电常数；σ表示当前电导率；αε、ασ、βε和βσ表示值为正的正则化参数；||·||TV表示全变差算子；将反演目标函数
Figure BDA00033199795200000323
进一步分解为子问题1和子问题2，子问题1为求解模型介质参数向量，子问题2为求解先验模型向量。
所述的子问题1具体为：
Figure BDA00033199795200000324
其中，
Figure BDA00033199795200000325
表示模型介质参数向量，
Figure BDA00033199795200000326
表示当前介电常数对应的模型介质参数向量；
Figure BDA0003319979520000041
表示当前电导率对应的模型介质参数向量；
Figure BDA0003319979520000042
表示第k次迭代后的模型介质参数向量；
Figure BDA0003319979520000043
表示关于模型介质参数向量的反演目标函数；
Figure BDA0003319979520000044
表示数据目标函数；ε表示当前介电常数；σ表示当前电导率；
Figure BDA0003319979520000045
和
Figure BDA0003319979520000046
表示第k次迭代、值为正的正则化参数；
Figure BDA0003319979520000047
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure BDA0003319979520000048
表示第k-1次迭代后当前电导率对应的先导模型向量；其中，第k次迭代的正则化参数
Figure BDA0003319979520000049
和
Figure BDA00033199795200000410
分别为：
Figure BDA00033199795200000411
和
Figure BDA00033199795200000412
其中，
Figure BDA00033199795200000413
和
Figure BDA00033199795200000414
为无量纲数；
Figure BDA00033199795200000415
表示数据目标函数；
Figure BDA00033199795200000416
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000417
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200000418
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure BDA00033199795200000419
表示第k-1次迭代后当前电导率对应的先导模型向量。
所述的子问题2包括，子问题2：
Figure BDA00033199795200000420
其中，ε表示当前介电常数；σ表示当前电导率；
Figure BDA00033199795200000421
βσ和βε表示值为正的正则化参数；||·||TV表示全变差算子；
Figure BDA00033199795200000422
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200000423
表示当前介电常数对应的先导模型向量；
Figure BDA00033199795200000424
表示关于当前介电常数对应的先导模型向量的目标函数；
Figure BDA00033199795200000425
表示第k次迭代后的当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000426
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200000427
表示当前介电常数对应的先导模型向量；
Figure BDA00033199795200000428
表示关于当前电导率对应的先导模型向量的目标函数；
Figure BDA0003319979520000051
表示第k次迭代后的当前介电常数对应的模型介质参数向量；其中，第k次迭代的正则化参数
Figure BDA0003319979520000052
和
Figure BDA0003319979520000053
分别为：
Figure BDA0003319979520000054
和
Figure BDA0003319979520000055
其中，
Figure BDA0003319979520000056
和
Figure BDA0003319979520000057
为无量纲数；
Figure BDA0003319979520000058
表示数据目标函数；
Figure BDA0003319979520000059
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000510
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200000511
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure BDA00033199795200000512
表示第k-1次迭代后当前电导率对应的先导模型向量。
所述的步骤S5具体包括：利用有限内存的BFGS方法对子问题1进行求解，对于第k+1次迭代后的模型介质参数向量
Figure BDA00033199795200000513
Figure BDA00033199795200000514
其中，
Figure BDA00033199795200000515
表示第k次迭代后的模型介质参数向量；a(k)为第k次迭代的搜索步长；
Figure BDA00033199795200000516
表示第k次迭代的搜索方向，
Figure BDA00033199795200000517
表示逆海森矩阵的近似；
Figure BDA00033199795200000518
表示第k次迭代的，关于模型介质参数向量的目标函数
Figure BDA00033199795200000519
的梯度；
第k次迭代的，关于模型介质参数向量的目标函数
Figure BDA00033199795200000520
的梯度
Figure BDA00033199795200000521
具体为：
Figure BDA00033199795200000522
其中，
Figure BDA00033199795200000523
表示关于模型介质参数向量的目标函数的梯度；
Figure BDA00033199795200000524
表示数据目标函数；
Figure BDA00033199795200000525
和
Figure BDA00033199795200000526
表示第k次迭代的值为正的正则化参数；
Figure BDA00033199795200000527
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000528
表示当前电导率对应的模型介质参数向量；
Figure BDA0003319979520000061
表示第k-1次迭代后当前介电常数对应的先验模型向量；
Figure BDA0003319979520000062
表示第k-1次迭代后当前电导率对应的先验模型向量；
Figure BDA0003319979520000063
表示关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度；
Figure BDA0003319979520000064
表示关于当前电导率对应的模型介质参数向量的数据目标函数的梯度；
对于地下介质单元e，关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度
Figure BDA0003319979520000065
采用伴随状态法求得：
Figure BDA0003319979520000066
其中，
Figure BDA0003319979520000067
为模型介质参数向量；Nω表示角频率的最大数量；Ns表示激励源的最大数量；i表示频率编号；r表示激励源编号；
Figure BDA0003319979520000068
为观测数据与模拟数据之间的差值，具体为
Figure BDA0003319979520000069
表示观测数据；
Figure BDA00033199795200000610
表示模拟数据；T表示转置，*表示共轭；
Figure BDA00033199795200000611
表示实部算子；
Figure BDA00033199795200000612
表示离散的电场值；
Figure BDA00033199795200000613
表示阻抗矩阵；
Figure BDA00033199795200000614
表示当前源的接收点的对应位置矩阵；
Figure BDA00033199795200000615
为一个稀疏矩阵，
Figure BDA00033199795200000616
表示阻抗矩阵
Figure BDA00033199795200000617
对当前单元e对应的模型介质参数me的敏感性，
Figure BDA00033199795200000618
为一个对角矩阵，当前单元e对应的系数为非零元素，通过
Figure BDA00033199795200000619
计算得到；其模型介质参数向量me包括介电常数和电导率；
求解单元矩阵
Figure BDA00033199795200000620
对于相应的当前单元e中当前介电常数ε对应的模型介质参数
Figure BDA00033199795200000621
的敏感性
Figure BDA00033199795200000622
求解单元矩阵
Figure BDA00033199795200000623
对于相应的当前单元e中当前电导率σ对应的模型介质参数
Figure BDA00033199795200000624
的敏感性
Figure BDA00033199795200000625
Figure BDA00033199795200000626
和
Figure BDA00033199795200000627
其中，
Figure BDA00033199795200000628
表示当前单元e中当前介电常数ε对应的模型介质参数；ω表示角频率；j表示虚数单位，
Figure BDA0003319979520000071
ε0为自由空间介电常数；
Figure BDA0003319979520000072
为矢量基函数，使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；
Figure BDA0003319979520000073
表示当前单元e内的体积积分；
Figure BDA0003319979520000074
表示当前单元e中当前电导率σ对应的模型介质参数；β表示值为正的正则化参数；σ0为模型背景电导率。
所述的步骤S6具体包括，利用Split-Bregman迭代方法求解子问题2，并进行更新：
子问题2采用统一TV模型表示：
Figure BDA0003319979520000075
其中，
Figure BDA0003319979520000076
表示先验模型向量，
Figure BDA0003319979520000077
表示当前介电常数对应的先验模型向量，
Figure BDA0003319979520000078
表示当前电导率对应的先验模型向量；
Figure BDA0003319979520000079
表示第k次迭代后的模型介质参数向量；
Figure BDA00033199795200000710
表示第k次迭代后当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200000711
表示第k次迭代后当前电导率对应的模型介质参数向量；||·||TV表示全变差算子；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
在全变差算子||·||TV中，对于一个3D各向同性的正则化算子
Figure BDA00033199795200000712
表示为：
Figure BDA00033199795200000713
其中，e表示当前单元；
Figure BDA00033199795200000714
表示先验模型向量在x轴上的梯度；
Figure BDA00033199795200000715
表示先验模型向量在y轴上的梯度；
Figure BDA00033199795200000716
表示先验模型向量在z轴上的梯度；设置辅助变量
Figure BDA00033199795200000717
和
Figure BDA00033199795200000718
分割问题的L1约束项和L2惩罚项分量，得到统一TV模型***的Bregman公式的最优性条件为：
Figure BDA00033199795200000719
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
Figure BDA00033199795200000720
表示单位矩阵；
Figure BDA00033199795200000721
表示第k次迭代后的模型介质参数向量；
Figure BDA00033199795200000722
表示三维下的x,y,z三个方向；
Figure BDA0003319979520000081
表示第k次迭代的、
Figure BDA0003319979520000082
方向上的辅助变量；
Figure BDA0003319979520000083
表示第k次迭代的、
Figure BDA0003319979520000084
方向上的Bregman系数，其迭代公式为
Figure BDA0003319979520000085
Figure BDA0003319979520000086
表示第k+1次迭代、
Figure BDA0003319979520000087
方向上的Bregman系数；
Figure BDA0003319979520000088
表示第k+1次迭代的、
Figure BDA0003319979520000089
方向上的辅助变量；
Figure BDA00033199795200000810
表示
Figure BDA00033199795200000811
与
Figure BDA00033199795200000812
的向量差在
Figure BDA00033199795200000813
方向上的梯度的转置矩阵；Δ为拉普拉斯算子；
Figure BDA00033199795200000814
表示第k+1次迭代后的先验模型向量；
设
Figure BDA00033199795200000815
辅助变量
Figure BDA00033199795200000816
使用广义收缩公显式求解：
Figure BDA00033199795200000817
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
Figure BDA00033199795200000818
表示收缩系数，
Figure BDA00033199795200000819
表示第k次迭代后的先验模型向量在
Figure BDA00033199795200000820
方向上的梯度；
Figure BDA00033199795200000821
表示第k次迭代后的先验模型向量在x方向上的梯度；
Figure BDA00033199795200000822
表示第k次迭代后的先验模型向量在y方向上的梯度；
Figure BDA00033199795200000823
表示第k次迭代后的先验模型向量在z方向上的梯度；
Figure BDA00033199795200000824
表示第k次迭代的、
Figure BDA00033199795200000825
方向上的Bregman系数；
Figure BDA00033199795200000826
表示第k次迭代的、x方向上的Bregman系数；
Figure BDA00033199795200000827
表示第k次迭代的、y方向上的Bregman系数；
Figure BDA00033199795200000828
表示第k次迭代的、z方向上的Bregman系数。
所述的步骤S7具体为，重复步骤S5-S6，对子问题1和子问题2同步交替进行更新，直到满足设定反演终止条件，结束反演并输出反演结果：设定反演终止条件包括达到设定迭代精度或达到最大的迭代次数；修正的全变差问题交替求解子问题1和子问题2；在第k次迭代中，根据步骤S4使用变量
Figure BDA00033199795200000829
和
Figure BDA00033199795200000830
求解第k次迭代后的模型介质参数向量
Figure BDA00033199795200000831
和
Figure BDA00033199795200000832
表示第k-1次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200000833
表示第k-1次迭代后当前电导率对应的先导模型向量；然后根据步骤S5使用第k次迭代后的模型介质参数向量
Figure BDA00033199795200000834
和
Figure BDA00033199795200000835
求解变量
Figure BDA0003319979520000091
和
Figure BDA0003319979520000092
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA0003319979520000093
表示第k次迭代后当前电导率对应的先导模型向量；得到的变量
Figure BDA0003319979520000094
和
Figure BDA0003319979520000095
在下一个迭代步骤中使用；对于第一个迭代步骤，初始当前介电常数对应的先导模型向量
Figure BDA0003319979520000096
等于初始当前介电常数对应的模型介质参数向量
Figure BDA0003319979520000097
初始当前电导率对应的先导模型向量
Figure BDA0003319979520000098
等于初始当前电导率对应的模型介质参数向量
Figure BDA0003319979520000099
得到反演序列：
Figure BDA00033199795200000910
Figure BDA00033199795200000911
本发明提供的这种三维频率域探地雷达双参数同步反演方法，采用有限元求解器来模拟GPR正演问题；使用拟牛顿法中的L-BFGS法求解反演问题，避免了Hessian矩阵的计算；通过参数变换和参数调整因子，便于在相同的下降步骤内能够同步更新介电常数和电导率值，以提高反演的效率；将改进全变差模型约束方法引入到GPR全波形反演过程中，增加算法的鲁棒性，结合频率权重因子，采用并行方式实现了3D的GPR数据的定量成像，使成像过程更加高效精确。
附图说明
图1为本发明方法的流程示意图。
图2a和图2b为本发明实施例的模型1的三维模型示意图。
图3a和图3b为本发明实施例的模型1的全波形反演结果示意图。
图4a和图4b为本发明实施例的模型1的相对数据残差和重构误差示意图。
图5a和图5b为本发明实施例的模型2的三维模型示意图。
图6a、图6b、图6c和图6d为本发明实施例的模型2的不同正则化策略下的全波形反演结果示意图。
具体实施方式
如图1为本发明方法的流程示意图：本发明提供的这种三维频率域探地雷达双参数同步反演方法，包括如下步骤：
S1.输入观测数据，设置三维初始模型，并初始化参数；
S2.对建立的三维初始模型正演；
S3.构建反演目标函数，并分解为两个交替的子问题，包括子问题1和子问题2，子问题1为求解模型介质参数向量，子问题2为求解先导模型向量；
S4.初始化介电常数和电导率，得到初始先导模型向量；
S5.根据先导模型向量求解子问题1，并更新模型介质参数向量；
S6.根据得到的模型介质参数向量求解子问题2，并更新先导模型向量；
S7.重复步骤S5-S6，同步交替进行更新，直到满足反演终止条件，结束反演并输出反演结果。
所述的步骤S2包括，利用有限单元法对建立的三维初始模型进行正演，具体为采用线性的矢量基函数的有限单元法进行模拟区域的离散，剖分网格为矩形块单元，截断边界处采用各向异性完全匹配层：三维GPR满足的频率域麦克斯韦方程表示为：
Figure BDA0003319979520000101
其中，
Figure BDA0003319979520000102
表示电场强度；
Figure BDA0003319979520000103
表示磁场强度；ε表示当前介电常数；μ表示磁导率；σ表示当前电导率；ω表示角频率，ω＝2πf，f表示频率；j表示虚数单位，
Figure BDA0003319979520000104
表示电流密度；本方程的时谐因子为ejωt；
选用线性的矢量基函数的有限单元法进行离散，截断边界处采用各向异性完全匹配层，经过总体合并得到多源***的线性方程组：
Figure BDA0003319979520000105
其中，
Figure BDA0003319979520000106
表示与频率相关的所有源公用的系数矩阵；
Figure BDA0003319979520000107
表示离散电场值；
Figure BDA0003319979520000108
表示离散化激励源项；ω表示角频率；ε表示介电常数；σ表示电导率；s表示当前位置的激励源；与频率相关的所有源公用的系数矩阵
Figure BDA0003319979520000109
由所有单元矩阵
Figure BDA00033199795200001010
合成：
Figure BDA00033199795200001011
其中，
Figure BDA00033199795200001012
表示当前单元e内的体积积分；
Figure BDA00033199795200001013
和
Figure BDA00033199795200001014
为矢量基函数；使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；此处采用了容易实现的矩形块单元。线性***使用直接求解器MUMPS求解，对于给定的介质和频率，仅需要对矩阵
Figure BDA00033199795200001015
进行一次LU分解，有效地实现多个右端项问题求解。
对于单源单个频率的正演算子可描述为：
Figure BDA0003319979520000111
其中，
Figure BDA0003319979520000112
表示角频率ω和当前激励源s下的模拟数据；
Figure BDA0003319979520000113
表示当前源的接收点的对应位置矩阵；
Figure BDA0003319979520000114
表示离散电场值；
Figure BDA0003319979520000115
表示与频率相关的所有源公用的系数矩阵；
Figure BDA0003319979520000116
表示离散化激励源项。
所述的步骤S3具体包括，探地雷达全波形反演实质是利用已知的实测数据来重构地下介质中的模型参数，包括介电常数ε和电导率σ的空间分布；根据正演模拟数据和实测数据之间的拟合最优，数据目标函数
Figure BDA0003319979520000117
定义为：
Figure BDA0003319979520000118
其中，
Figure BDA0003319979520000119
为模型介质参数向量；Nω表示角频率的数量；Ns表示激励源的数量；ηi表示多频率数据加权因子；ωi表示第i个频率；
Figure BDA00033199795200001110
表示第j个激励源；
Figure BDA00033199795200001111
为观测数据与模拟数据之间的差值，具体为
Figure BDA00033199795200001112
表示观测数据；
Figure BDA00033199795200001113
表示模拟数据；T表示转置，*表示共轭；ηi表示多频率数据加权因子：
Figure BDA00033199795200001114
由于介电常数、电导率在数量级上相差很大，给反演计算带来了诸多不便；为了在每次迭代时同时更新介电常数和电导率模型，需考虑如何平衡不同参数的影响；变化之后，模型介质参数向量
Figure BDA00033199795200001115
为：
Figure BDA00033199795200001116
其中，
Figure BDA00033199795200001117
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200001118
表示当前电导率对应的模型介质参数向量；εr表示相对介电常数，εr＝ε/ε0，ε0为自由空间介电常数，ε为当前介电常数；σr表示相对电导率，σr＝σ/σ0，σ0为模型背景电导率，模型背景电导率取设定参考介质的电导率，σ为当前电导率；к表示参数调整因子，用于在优化过程中通过控制相对电导率σr对相对介电常数εr的权重，避免由相对电导率与相对介电常数定义不准确引起反演过程的不确定性。
由于GPR的FEI属于不适定性问题，为了提高FWI的鲁棒性，在反演过程中常引入改进的全变差正则化方法，则反演目标函数
Figure BDA0003319979520000121
定义为：
Figure BDA0003319979520000122
其中，
Figure BDA0003319979520000123
表示模型介质参数向量；
Figure BDA0003319979520000124
表示当前介电常数对应的模型介质参数向量；
Figure BDA0003319979520000125
表示当前电导率对应的模型介质参数向量；
Figure BDA0003319979520000126
表示先验模型向量，
Figure BDA0003319979520000127
表示当前介电常数对应的先验模型向量，
Figure BDA0003319979520000128
表示当前电导率对应的先验模型向量；
Figure BDA0003319979520000129
表示数据目标函数；
Figure BDA00033199795200001210
表示当前介电常数对应的模型介质参数向量；ε表示当前介电常数；σ表示当前电导率；αε、ασ、βε和βσ表示值为正的正则化参数；||·||TV表示全变差算子；将反演目标函数
Figure BDA00033199795200001211
进一步分解为2个交替最小化子问题：
子问题1：
Figure BDA00033199795200001212
子问题2：
Figure BDA00033199795200001213
其中，
Figure BDA00033199795200001214
表示模型介质参数向量；
Figure BDA00033199795200001215
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200001216
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200001217
表示第k次迭代后的模型介质参数向量；
Figure BDA00033199795200001218
表示关于模型介质参数向量的目标函数；
Figure BDA00033199795200001219
表示数据目标函数；ε表示当前介电常数；σ表示当前电导率；αε、ασ、αε,σ、βσ和βε表示值为正的正则化参数；
Figure BDA00033199795200001220
表示第k-1次迭代后当前介电常数对应的先验模型向量，
Figure BDA00033199795200001221
表示第k-1次迭代后当前电导率对应的先验模型向量；||·||TV表示全变差算子；
Figure BDA0003319979520000131
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA0003319979520000132
表示当前介电常数对应的先导模型向量；
Figure BDA0003319979520000133
表示关于当前介电常数对应的先导模型向量的目标函数；
Figure BDA0003319979520000134
表示第k次迭代后的当前介电常数对应的模型介质参数向量；
Figure BDA0003319979520000135
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA0003319979520000136
表示当前介电常数对应的先导模型向量；
Figure BDA0003319979520000137
表示关于当前电导率对应的先导模型向量的目标函数；
Figure BDA0003319979520000138
表示第k次迭代后的当前介电常数对应的模型介质参数向量。
子问题1是使用基于Tikhonov正则化和多参数的先验模型，即第k-1次迭代后当前介电常数和当前电导率对应的先导模型向量
Figure BDA0003319979520000139
求解第k次迭代后的模型介质参数向量
Figure BDA00033199795200001310
的全波形反演问题，Tikhonov正则化反演能保证反演的稳定性；
子问题2是进行2次标准的L2-全变差最小化方法求解不同参数的、第k次迭代后当前介电常数对应的模型介质参数向量
Figure BDA00033199795200001311
和当前电导率对应的先验模型向量
Figure BDA00033199795200001312
以保持反演结果，反演结果即第k次迭代后当前介电常数对应的模型介质参数向量
Figure BDA00033199795200001313
和当前电导率对应的模型介质参数向量
Figure BDA00033199795200001314
的分界面清晰；交替迭代求解方程不仅能增强界面的清晰度，同时能保证全波形反演更加稳定。
其中，第k次迭代的正则化参数
Figure BDA00033199795200001315
和
Figure BDA00033199795200001316
分别为：
Figure BDA00033199795200001317
和
Figure BDA00033199795200001318
其中，
Figure BDA00033199795200001319
和
Figure BDA00033199795200001320
为无量纲数，在本实施例中，
Figure BDA00033199795200001321
表示数据目标函数；
Figure BDA00033199795200001322
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200001323
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200001324
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure BDA00033199795200001325
表示第k-1次迭代后当前电导率对应的先导模型向量；第k次迭代的正则化参数
Figure BDA0003319979520000141
和
Figure BDA0003319979520000142
的迭代方程可以实现正则化因子的自适应选取。
所述的步骤S5具体包括：利用有限内存的BFGS方法对子问题1进行求解，同步更新介电常数模型和电导率模型；
采用L-BFGS算法求解模型介质参数向量
Figure BDA0003319979520000143
的反问题，对于第k+1次迭代后的模型介质参数向量
Figure BDA0003319979520000144
Figure BDA0003319979520000145
其中，
Figure BDA0003319979520000146
表示第k次迭代后的模型介质参数向量；a(k)为第k次迭代的搜索步长，采用强Wolfe准则的非精确线搜素方法计算；
Figure BDA0003319979520000147
表示第k次迭代的搜索方向，
Figure BDA0003319979520000148
表示逆海森矩阵的近似，可以用两个循环递归算法计算；
Figure BDA0003319979520000149
表示第k次迭代的，关于模型介质参数向量的目标函数
Figure BDA00033199795200001410
的梯度：
Figure BDA00033199795200001411
其中，
Figure BDA00033199795200001412
表示关于模型介质参数向量的目标函数的梯度；
Figure BDA00033199795200001413
表示数据目标函数；
Figure BDA00033199795200001414
和
Figure BDA00033199795200001415
表示第k次迭代的值为正的正则化参数；
Figure BDA00033199795200001416
表示当前介电常数对应的模型介质参数向量；
Figure BDA00033199795200001417
表示当前电导率对应的模型介质参数向量；
Figure BDA00033199795200001418
表示第k-1次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200001419
表示第k-1次迭代后当前电导率对应的先导模型向量；
Figure BDA00033199795200001420
表示关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度；
Figure BDA00033199795200001421
表示关于当前电导率对应的模型介质参数向量的数据目标函数的梯度。
对于地下介质单元e，关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度
Figure BDA00033199795200001422
采用伴随状态法求得：
Figure BDA0003319979520000151
其中，
Figure BDA0003319979520000152
为模型介质参数向量；Nω表示角频率的最大数量；Ns表示激励源的最大数量；i表示频率编号；r表示激励源编号；
Figure BDA0003319979520000153
为观测数据与模拟数据之间的差值，具体为
Figure BDA0003319979520000154
表示观测数据；
Figure BDA0003319979520000155
表示模拟数据；T表示转置，*表示共轭；
Figure BDA0003319979520000156
表示实部算子；
Figure BDA0003319979520000157
表示离散的电场值；
Figure BDA0003319979520000158
表示阻抗矩阵；
Figure BDA0003319979520000159
表示当前源的接收点的对应位置矩阵；
Figure BDA00033199795200001510
为一个稀疏矩阵，
Figure BDA00033199795200001511
表示阻抗矩阵
Figure BDA00033199795200001512
对当前单元e对应的模型介质参数me的敏感性，
Figure BDA00033199795200001513
为一个对角矩阵，只有当前单元e对应的系数为非零元素，可以通过
Figure BDA00033199795200001514
计算得到。其模型介质参数向量me包括介电常数和电导率。
求解单元矩阵
Figure BDA00033199795200001515
对于相应的当前单元e中当前介电常数ε对应的模型介质参数
Figure BDA00033199795200001516
的敏感性
Figure BDA00033199795200001517
求解单元矩阵
Figure BDA00033199795200001518
对于相应的当前单元e中当前电导率σ对应的模型介质参数
Figure BDA00033199795200001519
的敏感性
Figure BDA00033199795200001520
Figure BDA00033199795200001521
其中，
Figure BDA00033199795200001522
表示当前单元e中当前介电常数ε对应的模型介质参数；ω表示角频率；j表示虚数单位，
Figure BDA00033199795200001523
ε0为自由空间介电常数；
Figure BDA00033199795200001524
为矢量基函数，使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；
Figure BDA00033199795200001525
表示当前单元e内的体积积分；
Figure BDA00033199795200001526
表示当前单元e中当前电导率σ对应的模型介质参数；β表示值为正的正则化参数；σ0为模型背景电导率。
所述的步骤S6具体为，利用Split-Bregman迭代方法求解子问题2，即两个参数的辅助变量，并对其进行更新，可以采用如下步骤：
子问题2采用如下统一TV模型表示：
Figure BDA0003319979520000161
其中，
Figure BDA0003319979520000162
表示先验模型向量，
Figure BDA0003319979520000163
表示当前介电常数对应的先验模型向量，
Figure BDA0003319979520000164
表示当前电导率对应的先验模型向量；
Figure BDA0003319979520000165
表示第k次迭代后的模型介质参数向量；
Figure BDA0003319979520000166
表示第k次迭代后当前介电常数对应的模型介质参数向量；
Figure BDA0003319979520000167
表示第k次迭代后当前电导率对应的模型介质参数向量；||·||TV表示全变差算子；μTV＝μTV,ε,μTV,σ，μTV,ε和μTV,σ表示正则化参数；
对于一个3D各向同性的正则化算子
Figure BDA0003319979520000168
表示为：
Figure BDA0003319979520000169
其中，e表示当前单元；
Figure BDA00033199795200001610
表示先验模型向量在x轴上的梯度；
Figure BDA00033199795200001611
表示先验模型向量在y轴上的梯度；
Figure BDA00033199795200001612
表示先验模型向量在z轴上的梯度；子问题2不能直接求解
Figure BDA00033199795200001613
因此需要通过设置辅助变量
Figure BDA00033199795200001614
和
Figure BDA00033199795200001615
分割问题的L1约束项和L2惩罚项分量，得到统一TV模型***的Bregman公式的最优性条件为：
Figure BDA00033199795200001616
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV,ε和μTV,σ表示正则化参数；
Figure BDA00033199795200001617
表示单位矩阵；
Figure BDA00033199795200001618
表示第k次迭代后的模型介质参数向量；
Figure BDA00033199795200001619
表示三维下的x,y,z三个方向；
Figure BDA00033199795200001620
表示第k次迭代的，
Figure BDA00033199795200001621
方向上的辅助变量；
Figure BDA00033199795200001622
表示第k次迭代的，
Figure BDA00033199795200001623
方向上的Bregman系数，其迭代公式为
Figure BDA00033199795200001624
Figure BDA00033199795200001625
表示第k+1次迭代，
Figure BDA00033199795200001626
方向上的Bregman系数；
Figure BDA00033199795200001627
表示第k+1次迭代的，
Figure BDA00033199795200001628
方向上的辅助变量；
Figure BDA00033199795200001629
表示
Figure BDA00033199795200001630
与
Figure BDA00033199795200001631
的向量差在
Figure BDA00033199795200001632
方向上的梯度的转置矩阵；Δ为拉普拉斯算子；
Figure BDA00033199795200001633
表示第k+1次迭代后的先验模型向量；
设
Figure BDA00033199795200001634
辅助变量
Figure BDA00033199795200001635
使用广义收缩公显式求解：
Figure BDA0003319979520000171
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
Figure BDA0003319979520000172
表示收缩系数，
Figure BDA0003319979520000173
表示第k次迭代后的先验模型向量在
Figure BDA0003319979520000174
即x、y或z，方向上的梯度；
Figure BDA0003319979520000175
表示第k次迭代后的先验模型向量在x方向上的梯度；
Figure BDA0003319979520000176
表示第k次迭代后的先验模型向量在y方向上的梯度；
Figure BDA0003319979520000177
表示第k次迭代后的先验模型向量在z方向上的梯度；
Figure BDA0003319979520000178
表示第k次迭代的，
Figure BDA0003319979520000179
方向上的Bregman系数；
Figure BDA00033199795200001710
表示第k次迭代的、x方向上的Bregman系数；
Figure BDA00033199795200001711
表示第k次迭代的、y方向上的Bregman系数；
Figure BDA00033199795200001712
表示第k次迭代的，z方向上的Bregman系数。
所述的步骤S7具体为，重复步骤S5-S6，对子问题1和子问题2同步交替进行更新，直到满足设定反演终止条件，结束反演并输出反演结果。设定反演终止条件包括满足迭代精度或达到最大的迭代次数。修正的全变差问题交替求解子问题1和子问题2；在第k次迭代中，根据步骤S4使用变量
Figure BDA00033199795200001713
和
Figure BDA00033199795200001714
求解第k次迭代后的模型介质参数向量
Figure BDA00033199795200001715
和
Figure BDA00033199795200001716
表示第k-1次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200001717
表示第k-1次迭代后当前电导率对应的先导模型向量；然后根据步骤S5使用第k次迭代后的模型介质参数向量
Figure BDA00033199795200001718
和
Figure BDA00033199795200001719
求解变量
Figure BDA00033199795200001720
和
Figure BDA00033199795200001721
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure BDA00033199795200001722
表示第k次迭代后当前电导率对应的先导模型向量；得到的变量
Figure BDA00033199795200001723
和
Figure BDA00033199795200001724
在下一个迭代步骤中使用；对于第一个迭代步骤，初始当前介电常数对应的先导模型向量
Figure BDA00033199795200001725
等于初始当前介电常数对应的模型介质参数向量
Figure BDA00033199795200001726
初始当前电导率对应的先导模型向量
Figure BDA00033199795200001727
等于初始当前电导率对应的模型介质参数向量
Figure BDA00033199795200001728
得到反演序列：
Figure BDA00033199795200001729
Figure BDA0003319979520000181
以下结合一个实施案例，对本发明方法进行进一步说明：
如图2a和图2b为本发明实施例的模型1的三维模型示意图。图2a表示本发明实施例的模型1的介电常数εr的三维模型示意图；图2b表示本发明实施例的模型1的电导率的三维模型示意图。地表上方存在0.48m厚的空气层，地下背景介质的介电常数为εr＝4和σ＝3mS/m。地下0.72m深有4个大小为0.72m×0.72m×0.72m的异常体，两个蓝色异常体的电性参数为εr＝1和σ＝0；红色异常体为εr＝8和σ＝10mS/m，计算域由50×50×50网格组成，网格尺寸为0.12m×0.12m×0.12m，PML厚度为0.48m。正演模拟的激励源采用主频为200MHz的Ricker子波，在据地表0.12m的空气层中平均分布121个源和420个接收器(分别用黑色和白色“·”表示)，间距分别为0.48m和0.24m，每个源的信号将由所有接收器记录。
该异常体具有较大的对比度，4个异常体之间产生的多重散射以及地面观测方式均增加了反演的复杂性。将5％的高斯噪声添加到观测数据中，反演初始模型为背景均匀介质，反演过程中限定模型中空气层的参数，同时避免源和接收器位置处的奇异性，待求的地下介电常数和电导率的参数为210000个。
该全波形反演试验在Dell Precision T7920工作站上执行，具有双Intel(R)Xeon(R)Platinum 8168CPU，@2.70GHz，物理内存512GB，Windows 10操作***。反演过程中的正演模拟和梯度计算模块在频率上进行了CPU并行加速。选择10个频率数据(即50、60、70、80、100、120、150、180、200和220MHz)进行同步多频反演，经过100次迭代后的相对介电常数和电导率的模型结果如图3所示。反演过程中采用TV正则化，设置参数调节因子κ＝1.2，保证反演过程由介电常数起主要引导作用。反演耗时约16.87小时，在整个算法中使用的最大内存为300GB。
如图3a和图3b为本发明实施例的模型1的全波形反演结果示意图。图3a表示本发明实施例的模型1的介电常数全波形反演结果；图3b表示本发明实施例的模型1的电导率的全波形反演结果。如图3a和图3b所示，4个异常体的近似形状和物性参数均被成功地重构。对比图2a与图3a：介电常数的反演结果界面比较清晰，低介电常数异常体的反演结果与真实值基本一致，高介电常数异常体与真实值结果有一定差异；4个异常体上方均出现了相反的伪像。对比图2b与图3b：电导率异常体在体积上有一定的变化，中心位置发生了上移。相对而言，介电常数的重构比电导率的重构更为精确，特别是异常体的界面较为清晰。
如图4a和图4b为本发明实施例的模型1的相对数据残差和重构误差示意图。图4a表示本发明实施例的模型1的相对数据残差；图4b表示本发明实施例的模型1的重构误差。在图4a中，相对数据残差在迭代初期迅速下降，在迭代的中后期朝着最小值下降，并且收敛速度更稳定。图4b中蓝色实线表示的介电常数重构误差曲线较红色虚线表示的电导率重构误差曲线整体要小。根据两个收敛曲线的趋势，在更多的迭代之后可以实现更精确地重构。然而，由于4个异常体的物性参数已得到了较好地重构，更多的迭代次数将会耗费更多的时间。反演结果表明，频率域加权方法的同步多频策略，能有效地重构地下介质的三维分布。
本发明方法的优势如下：
通过上述案例反演结果表明，本发明方法对地面多偏移量噪声数据具有良好的适应性。不仅可以同步反演介电常数和电导率，提高反演精度，并且在异常体界面的清晰度上有所提高。以下结合一个案例进行进一步说明：
如图5a和图5b为本发明实施例的模型2的三维模型示意图。图5a表示本发明实施例的模型2的介电常数的三维模型示意图；图5b表示本发明实施例的模型2的电导率的三维模型示意图。为了观察此算法对于不同深度异常体的反演效果，在图2a和图2b所示的模型1的基础上，改变了其中两个异常体的埋深至1.08m，其余两个保持不变，并且进一步减小四个异常体的距离，3D模型2的示意图如图5所示。其他剖分和测量参数与3D模型1保持一致。在总场合成数据增加SNR＝30dB的高斯噪声作为观测数据，其他反演参数与策略与3D模型1保持一致。该算法利用MTV算法提高反演，因此在此基础上开展了标准Tikhonov正则化与MTV正则化反演实验，以对比不同正则化方法在含噪数据的效果。
如图6a、图6b、图6c和图6d为本发明实施例的模型2的不同正则化策略下的全波形反演结果示意图。图6a表示本发明实施例的模型2使用MTV正则化的介电常数的三维模型示意图；图6b表示本发明实施例的模型2使用标准Tikhonov正则化的介电常数的三维模型示意图；图6c表示本发明实施例的模型2使用MTV正则化的电导率的三维模型示意图；图6d表示本发明实施例的模型2使用标准Tikhonov正则化的电导率的三维模型示意图。反演中MTV正则化反演耗时约为15.91小时，标准Tikhonov正则化反演的反演耗时约为21.57小时(比MTV正则化方法长36％)，在整个算法中使用的最大内存为300GB。
观察图6可以发现，图6a-图6d中的异常体界面不如图5中的真实模型的界面清晰。图6a及图6c中的界面明显比图6b及图6d中的界面清晰，这表明了FWI和MTV正则化在保持尖锐界面方面比标准Tikhonov正则化更具优势。比较图6a、图6b和图6c、图6d中的结果可发现，相对而言，介电常数的重构比电导率的重构更为精确，这是由于数据对介电常数的灵敏度比电导率高。将反演结果图6c、图6d与图5b的真实电导率模型进行比较，发现电导率异常体在体积上有一定的变化，中心位置发生了上移，此外浅部的两个异常体反演参数更接近真实值。
值得指出的是，与标准Tikhonov正则化相比，MTV正则化反演结果的背景值更接近真实模型。总的来说，反演结果中没有明显的人为噪声，验证了该方法的鲁棒性。为了进一步定量评估反演结果的精度，使用了两种不同的重构误差指标：
Figure BDA0003319979520000201
Figure BDA0003319979520000202
其中，
Figure BDA0003319979520000203
表示标准重构误差；
Figure BDA0003319979520000204
表示相对重构误差；
Figure BDA0003319979520000205
表示真实模型中的介质参数向量；
Figure BDA0003319979520000206
表示最终反演得到的介质参数向量。
表1不同正则化策略下FWI的迭代次数及模型误差参数
方法 迭代次数 MSE(ε<sub>r</sub>) MSE(σ) Err<sub>m</sub>(ε<sub>r</sub>) Err<sub>m</sub>(σ)
MTV 72 0.058 0.149 0.117 0.240
Tikhonov 71 0.061 0.172 0.252 0.595
在表1中，列出了两种正则化方案的迭代次数及模型误差参数。结果表明，两种正则化方案的介电常数重构误差均小于电导率重构误差；与标准Tikhonov正则化相比，MTV正则化的FWI具有更低的重构误差。
两种正则化策略下的反演结果表明，基于MTV正则化的方法能够提高介电常数和电导率反演的精度，在异常体界面的清晰度方面优于经典的Tikhonov正则化方法。

Claims (8)

1.一种三维频率域探地雷达双参数同步反演方法，其特征在于包括如下步骤：
S1.输入观测数据，设置三维初始模型，并初始化参数；
S2.对建立的三维初始模型正演；
S3.构建反演目标函数，并分解为两个交替的子问题，包括子问题1和子问题2，子问题1为求解模型介质参数向量，子问题2为求解先导模型向量；
S4.初始化介电常数和电导率，得到初始先导模型向量；
S5.根据先导模型向量求解子问题1，并更新模型介质参数向量；
S6.根据得到的模型介质参数向量求解子问题2，并更新先导模型向量；
S7.重复步骤S5-S6，同步交替进行更新，直到满足反演终止条件，结束反演并输出反演结果。
2.根据权利要求1所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的步骤S2包括，利用有限单元法对建立的三维初始模型进行正演，具体为采用线性的矢量基函数的有限单元法进行模拟区域的离散，剖分网格为矩形块单元，截断边界处采用各向异性完全匹配层；三维GPR满足的频率域麦克斯韦方程表示为：
Figure FDA0003319979510000011
其中，
Figure FDA0003319979510000012
表示电场强度；ε表示当前介电常数；μ表示磁导率；σ表示当前电导率；ω表示角频率，ω＝2πf，f表示频率；j表示虚数单位，
Figure FDA0003319979510000013
Figure FDA0003319979510000014
为电流密度。
求得多源***的线性方程组：
Figure FDA0003319979510000015
其中，
Figure FDA0003319979510000016
表示与频率相关的所有源公用的系数矩阵；
Figure FDA0003319979510000017
表示离散电场值；
Figure FDA0003319979510000018
表示离散化激励源项；ω表示角频率；ε表示介电常数；σ表示电导率；s表示当前位置的激励源；与频率相关的所有源公用的系数矩阵
Figure FDA0003319979510000019
由所有单元矩阵
Figure FDA00033199795100000110
合成：
Figure FDA00033199795100000111
其中，<·,·>Ωe表示当前单元e内的体积积分；
Figure FDA0003319979510000021
和
Figure FDA0003319979510000022
为矢量基函数；使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；对矩阵
Figure FDA0003319979510000023
进行一次LU分解；对于单源单个频率的正演算子描述为：
Figure FDA0003319979510000024
其中，
Figure FDA0003319979510000025
表示角频率ω和当前激励源s下的模拟数据；
Figure FDA0003319979510000026
表示当前源的接收点的对应位置矩阵；
Figure FDA0003319979510000027
表示离散电场值；
Figure FDA0003319979510000028
表示与频率相关的所有源公用的系数矩阵；
Figure FDA0003319979510000029
表示离散化激励源项。
3.根据权利要求2所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的步骤S3，利用已知的实测数据重构地下介质中的模型参数，包括介电常数ε和电导率σ的空间分布；并引入全变差正则化方法，反演目标函数
Figure FDA00033199795100000210
定义为：
Figure FDA00033199795100000211
其中，
Figure FDA00033199795100000212
表示模型介质参数向量，
Figure FDA00033199795100000213
Figure FDA00033199795100000214
表示当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000215
表示当前电导率对应的模型介质参数向量；
Figure FDA00033199795100000216
表示先验模型向量，
Figure FDA00033199795100000217
Figure FDA00033199795100000218
表示当前介电常数对应的先验模型向量，
Figure FDA00033199795100000219
表示当前电导率对应的先验模型向量；
Figure FDA00033199795100000220
表示数据目标函数；ε表示当前介电常数；σ表示当前电导率；αε、ασ、βε和βσ表示值为正的正则化参数；||·||TV表示全变差算子；将反演目标函数
Figure FDA00033199795100000221
进一步分解为子问题1和子问题2，子问题1为求解模型介质参数向量，子问题2为求解先验模型向量。
4.根据权利要求3所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的子问题1具体为：
Figure FDA00033199795100000222
其中，
Figure FDA00033199795100000223
表示模型介质参数向量，
Figure FDA00033199795100000224
Figure FDA00033199795100000225
表示当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000226
表示当前电导率对应的模型介质参数向量；
Figure FDA00033199795100000227
表示第k次迭代后的模型介质参数向量；
Figure FDA0003319979510000031
表示关于模型介质参数向量的反演目标函数；
Figure FDA0003319979510000032
表示数据目标函数；ε表示当前介电常数；σ表示当前电导率；
Figure FDA0003319979510000033
和
Figure FDA0003319979510000034
表示第k次迭代、值为正的正则化参数；
Figure FDA0003319979510000035
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure FDA0003319979510000036
表示第k-1次迭代后当前电导率对应的先导模型向量；其中，第k次迭代的正则化参数
Figure FDA0003319979510000037
和
Figure FDA0003319979510000038
分别为：
Figure FDA0003319979510000039
其中，
Figure FDA00033199795100000310
和
Figure FDA00033199795100000311
为无量纲数；
Figure FDA00033199795100000312
表示数据目标函数；
Figure FDA00033199795100000313
表示当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000314
表示当前电导率对应的模型介质参数向量；
Figure FDA00033199795100000315
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure FDA00033199795100000316
表示第k-1次迭代后当前电导率对应的先导模型向量。
5.根据权利要求3所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的子问题2包括，子问题2：
Figure FDA00033199795100000317
其中，ε表示当前介电常数；σ表示当前电导率；
Figure FDA00033199795100000318
βσ和βε表示值为正的正则化参数；||·||TV表示全变差算子；
Figure FDA00033199795100000319
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure FDA00033199795100000320
表示当前介电常数对应的先导模型向量；
Figure FDA00033199795100000321
表示关于当前介电常数对应的先导模型向量的目标函数；
Figure FDA00033199795100000322
表示第k次迭代后的当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000323
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure FDA00033199795100000324
表示当前介电常数对应的先导模型向量；
Figure FDA00033199795100000325
表示关于当前电导率对应的先导模型向量的目标函数；
Figure FDA00033199795100000326
表示第k次迭代后的当前介电常数对应的模型介质参数向量；其中，第k次迭代的正则化参数
Figure FDA0003319979510000041
和
Figure FDA0003319979510000042
分别为：
Figure FDA0003319979510000043
其中，
Figure FDA0003319979510000044
和
Figure FDA0003319979510000045
为无量纲数；
Figure FDA0003319979510000046
表示数据目标函数；
Figure FDA0003319979510000047
表示当前介电常数对应的模型介质参数向量；
Figure FDA0003319979510000048
表示当前电导率对应的模型介质参数向量；
Figure FDA0003319979510000049
表示第k-1次迭代后当前介电常数对应的先导模型向量，
Figure FDA00033199795100000410
表示第k-1次迭代后当前电导率对应的先导模型向量。
6.根据权利要求5所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的步骤S5具体包括：利用有限内存的BFGS方法对子问题1进行求解，对于第k+1次迭代后的模型介质参数向量
Figure FDA00033199795100000411
Figure FDA00033199795100000412
其中，
Figure FDA00033199795100000413
表示第k次迭代后的模型介质参数向量；a(k)为第k次迭代的搜索步长；
Figure FDA00033199795100000414
表示第k次迭代的搜索方向，
Figure FDA00033199795100000415
Figure FDA00033199795100000416
表示逆海森矩阵的近似；
Figure FDA00033199795100000417
表示第k次迭代的，关于模型介质参数向量的目标函数
Figure FDA00033199795100000418
的梯度；
第k次迭代的，关于模型介质参数向量的目标函数
Figure FDA00033199795100000419
的梯度
Figure FDA00033199795100000420
具体为：
Figure FDA00033199795100000421
其中，
Figure FDA00033199795100000422
表示关于模型介质参数向量的目标函数的梯度；
Figure FDA00033199795100000423
表示数据目标函数；
Figure FDA00033199795100000424
和
Figure FDA00033199795100000425
表示第k次迭代的值为正的正则化参数；
Figure FDA00033199795100000426
表示当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000427
表示当前电导率对应的模型介质参数向量；
Figure FDA0003319979510000051
表示第k-1次迭代后当前介电常数对应的先验模型向量；
Figure FDA0003319979510000052
表示第k-1次迭代后当前电导率对应的先验模型向量；
Figure FDA0003319979510000053
表示关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度；
Figure FDA0003319979510000054
表示关于当前电导率对应的模型介质参数向量的数据目标函数的梯度；
对于地下介质单元e，关于当前介电常数对应的模型介质参数向量的数据目标函数的梯度
Figure FDA0003319979510000055
采用伴随状态法求得：
Figure FDA0003319979510000056
其中，
Figure FDA0003319979510000057
为模型介质参数向量；Nω表示角频率的最大数量；Ns表示激励源的最大数量；i表示频率编号；r表示激励源编号；
Figure FDA0003319979510000058
为观测数据与模拟数据之间的差值，具体为
Figure FDA0003319979510000059
Figure FDA00033199795100000510
表示观测数据；
Figure FDA00033199795100000511
表示模拟数据；T表示转置，*表示共轭；
Figure FDA00033199795100000512
表示实部算子；
Figure FDA00033199795100000513
表示离散的电场值；
Figure FDA00033199795100000514
表示阻抗矩阵；
Figure FDA00033199795100000515
表示当前源的接收点的对应位置矩阵；
Figure FDA00033199795100000516
为一个稀疏矩阵，
Figure FDA00033199795100000517
表示阻抗矩阵
Figure FDA00033199795100000518
对当前单元e对应的模型介质参数me的敏感性，
Figure FDA00033199795100000519
为一个对角矩阵，当前单元e对应的系数为非零元素，通过
Figure FDA00033199795100000520
计算得到；其模型介质参数向量me包括介电常数和电导率；
求解单元矩阵
Figure FDA00033199795100000521
对于相应的当前单元e中当前介电常数ε对应的模型介质参数
Figure FDA00033199795100000522
的敏感性
Figure FDA00033199795100000523
求解单元矩阵
Figure FDA00033199795100000524
对于相应的当前单元e中当前电导率σ对应的模型介质参数
Figure FDA00033199795100000525
的敏感性
Figure FDA00033199795100000526
Figure FDA00033199795100000527
其中，
Figure FDA00033199795100000528
表示当前单元e中当前介电常数ε对应的模型介质参数；ω表示角频率；j表示虚数单位，
Figure FDA0003319979510000061
ε0为自由空间介电常数；
Figure FDA0003319979510000062
为矢量基函数，使用六面体单元对计算区域进行离散时：l,n＝1,2,…,12；
Figure FDA0003319979510000063
表示当前单元e内的体积积分；
Figure FDA0003319979510000064
表示当前单元e中当前电导率σ对应的模型介质参数；β表示值为正的正则化参数；σ0为模型背景电导率。
7.根据权利要求6所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的步骤S6具体包括，利用Split-Bregman迭代方法求解子问题2，并进行更新：
子问题2采用统一TV模型表示：
Figure FDA0003319979510000065
其中，
Figure FDA0003319979510000066
表示先验模型向量，
Figure FDA0003319979510000067
Figure FDA0003319979510000068
表示当前介电常数对应的先验模型向量，
Figure FDA0003319979510000069
表示当前电导率对应的先验模型向量；
Figure FDA00033199795100000610
表示第k次迭代后的模型介质参数向量；
Figure FDA00033199795100000611
Figure FDA00033199795100000612
表示第k次迭代后当前介电常数对应的模型介质参数向量；
Figure FDA00033199795100000613
表示第k次迭代后当前电导率对应的模型介质参数向量；||·||TV表示全变差算子；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
在全变差算子||·||TV中，对于一个3D各向同性的正则化算子
Figure FDA00033199795100000614
表示为：
Figure FDA00033199795100000615
其中，e表示当前单元；
Figure FDA00033199795100000616
表示先验模型向量在x轴上的梯度；
Figure FDA00033199795100000617
表示先验模型向量在y轴上的梯度；
Figure FDA00033199795100000618
表示先验模型向量在z轴上的梯度；设置辅助变量
Figure FDA00033199795100000619
和
Figure FDA00033199795100000620
分割问题的L1约束项和L2惩罚项分量，得到统一TV模型***的Bregman公式的最优性条件为：
Figure FDA00033199795100000621
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
Figure FDA00033199795100000622
表示单位矩阵；
Figure FDA0003319979510000071
表示第k次迭代后的模型介质参数向量；
Figure FDA00033199795100000729
z表示三维下的x,y,z三个方向；
Figure FDA0003319979510000072
表示第k次迭代的、
Figure FDA00033199795100000732
方向上的辅助变量；
Figure FDA0003319979510000073
表示第k次迭代的、
Figure FDA00033199795100000730
方向上的Bregman系数，其迭代公式为
Figure FDA0003319979510000074
Figure FDA0003319979510000075
Figure FDA0003319979510000076
表示第k+1次迭代、
Figure FDA00033199795100000731
方向上的Bregman系数；
Figure FDA0003319979510000077
表示第k+1次迭代的、
Figure FDA00033199795100000734
方向上的辅助变量；
Figure FDA0003319979510000078
表示
Figure FDA0003319979510000079
与
Figure FDA00033199795100000710
的向量差在
Figure FDA00033199795100000733
方向上的梯度的转置矩阵；Δ为拉普拉斯算子；
Figure FDA00033199795100000711
表示第k+1次迭代后的先验模型向量；
设
Figure FDA00033199795100000735
辅助变量
Figure FDA00033199795100000712
使用广义收缩公显式求解：
Figure FDA00033199795100000713
其中，λ＝2μTV；μTV＝μTV,ε,μTV,σ，同时μTV、μTV,ε和μTV,σ表示正则化参数；
Figure FDA00033199795100000714
表示收缩系数，
Figure FDA00033199795100000715
Figure FDA00033199795100000716
表示第k次迭代后的先验模型向量在
Figure FDA00033199795100000738
方向上的梯度；
Figure FDA00033199795100000717
表示第k次迭代后的先验模型向量在x方向上的梯度；
Figure FDA00033199795100000718
表示第k次迭代后的先验模型向量在y方向上的梯度；
Figure FDA00033199795100000719
表示第k次迭代后的先验模型向量在z方向上的梯度；
Figure FDA00033199795100000720
表示第k次迭代的、
Figure FDA00033199795100000737
方向上的Bregman系数；
Figure FDA00033199795100000721
表示第k次迭代的、x方向上的Bregman系数；
Figure FDA00033199795100000722
表示第k次迭代的、y方向上的Bregman系数；
Figure FDA00033199795100000723
表示第k次迭代的、z方向上的Bregman系数。
8.根据权利要求7所述的三维频率域探地雷达双参数同步反演方法，其特征在于所述的步骤S7具体为，重复步骤S5-S6，对子问题1和子问题2同步交替进行更新，直到满足设定反演终止条件，结束反演并输出反演结果：设定反演终止条件包括达到设定迭代精度或达到最大的迭代次数；修正的全变差问题交替求解子问题1和子问题2；在第k次迭代中，根据步骤S4使用变量
Figure FDA00033199795100000724
和
Figure FDA00033199795100000725
求解第k次迭代后的模型介质参数向量
Figure FDA00033199795100000726
和
Figure FDA00033199795100000727
Figure FDA00033199795100000728
表示第k-1次迭代后当前介电常数对应的先导模型向量；
Figure FDA0003319979510000081
表示第k-1次迭代后当前电导率对应的先导模型向量；然后根据步骤S5使用第k次迭代后的模型介质参数向量
Figure FDA0003319979510000082
和
Figure FDA0003319979510000083
求解变量
Figure FDA0003319979510000084
和
Figure FDA0003319979510000085
Figure FDA0003319979510000086
表示第k次迭代后当前介电常数对应的先导模型向量；
Figure FDA0003319979510000087
表示第k次迭代后当前电导率对应的先导模型向量；得到的变量
Figure FDA0003319979510000088
和
Figure FDA0003319979510000089
在下一个迭代步骤中使用；对于第一个迭代步骤，初始当前介电常数对应的先导模型向量
Figure FDA00033199795100000810
等于初始当前介电常数对应的模型介质参数向量
Figure FDA00033199795100000811
初始当前电导率对应的先导模型向量
Figure FDA00033199795100000812
等于初始当前电导率对应的模型介质参数向量
Figure FDA00033199795100000813
得到反演序列：
Figure FDA00033199795100000814
Figure FDA00033199795100000815
CN202111243358.3A 2021-10-25 2021-10-25 一种三维频率域探地雷达双参数同步反演方法 Pending CN113970732A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111243358.3A CN113970732A (zh) 2021-10-25 2021-10-25 一种三维频率域探地雷达双参数同步反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111243358.3A CN113970732A (zh) 2021-10-25 2021-10-25 一种三维频率域探地雷达双参数同步反演方法

Publications (1)

Publication Number Publication Date
CN113970732A true CN113970732A (zh) 2022-01-25

Family

ID=79588240

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111243358.3A Pending CN113970732A (zh) 2021-10-25 2021-10-25 一种三维频率域探地雷达双参数同步反演方法

Country Status (1)

Country Link
CN (1) CN113970732A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115598714A (zh) * 2022-12-14 2023-01-13 西南交通大学（Cn） 基于时空耦合神经网络的探地雷达电磁波阻抗反演方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115598714A (zh) * 2022-12-14 2023-01-13 西南交通大学（Cn） 基于时空耦合神经网络的探地雷达电磁波阻抗反演方法

Similar Documents

Publication Publication Date Title
Grayver et al. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver
Rücker et al. Three-dimensional modelling and inversion of dc resistivity data incorporating topography—I. Modelling
CN110968826A (zh) 一种基于空间映射技术的大地电磁深度神经网络反演方法
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
WO2011139413A1 (en) Artifact reduction in method of iterative inversion of geophysical data
Wang et al. Multiparameter full-waveform inversion of 3-D on-ground GPR with a modified total variation regularization scheme
CN110852025B (zh) 一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法
CN112327374B (zh) Gpu探地雷达复杂介质dgtd正演方法
Song et al. Fast three-dimensional electromagnetic nonlinear inversion in layered media with a novel scattering approximation
CN113970732A (zh) 一种三维频率域探地雷达双参数同步反演方法
Vatankhah et al. Improving the use of the randomized singular value decomposition for the inversion of gravity and magnetic data
CN113569447B (zh) 一种基于舒尔补方法的瞬变电磁三维快速正演方法
Guo et al. Application of supervised descent method for transient EM data inversion
CN112346139B (zh) 一种重力数据多层等效源延拓与数据转换方法
Mannseth et al. Assimilating spatially dense data for subsurface applications—balancing information and degrees of freedom
Scheunert et al. A cut-&-paste strategy for the 3-D inversion of helicopter-borne electromagnetic data—I. 3-D inversion using the explicit Jacobian and a tensor-based formulation
CN115563791B (zh) 基于压缩感知重构的大地电磁数据反演方法
Varentsov Modern trends in the solution of forward and inverse 3D electromagnetic induction problems
CN116088048A (zh) 包含不确定性分析的各向异性介质多参数全波形反演方法
Tsourlos et al. Efficient 2D inversion of long ERT sections
CN109655910A (zh) 基于相位校正的探地雷达双参数全波形反演方法
Zhong et al. Electrical resistivity tomography with smooth sparse regularization
Pirani et al. Multi-frequency identification of defects in conducting media
CN113325482A (zh) 一种时间域电磁数据反演成像方法
Gebre et al. L0-norm gravity inversion with new depth weighting function and bound constraints

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