一种探地雷达数据时频组合全波形反演方法技术

本发明专利技术提供一种探地雷达数据时频组合全波形反演方法，包括以下步骤：输入参数：发射信号子波波形、预处理后的GPR观测数据、研究区域的介电常数模型和研究区域的电导率模型；设定双参数反演的反演频率；频率域和时间域联合反演；达到迭代总结束条件，输出最终的电导率模型和介电常数模型。本发明专利技术提供的探地雷达数据时频组合全波形反演方法具有以下优点：本发明专利技术提供的探地雷达数据时频组合全波形反演方法，该反演策略简单有效，充分利用了信号时间域和频率域的信息，可以使FWI在抗参数失衡(trade off)和计算效率上达到一个最优化状态。

【技术实现步骤摘要】
一种探地雷达数据时频组合全波形反演方法
本专利技术属于雷达数据处理
，具体涉及一种探地雷达数据时频组合全波形反演方法。
技术介绍
全波形反演(FullWaveformInversion,简称FWI)本质上是一个非线性的病态问题，其目标是得到观测数据和模拟数据之间残差的最小二乘解。该方法首先由地震勘探领域的学者提出，并分别被扩展到了时间域算法和频率域算法中。随着该技术的不断发展，FWI已经逐渐被应用到实际生产项目中。探地雷达(GroundPenetratingRadar，简称GPR)是一种高效、无损的可以对目标介质实现远程探测的地球物理方法。GPR利用发射出的高频电磁脉冲在地下介质中经过反射、折射以及吸收作用之后所测量到的接收信号来实现对探测区域的电性参数成像，并主要用于探测分析地下介质结构，物性参数的改变以及孔隙或者裂缝的出现，可应用于隧道质量检测、路面密实度检测、混凝土构件非破坏性检测、水文环境检测评定等多方面。由于GPR与地震勘探在物理本质上的差异，导致GPR数据的FWI距离实际应用还有较大距离。这主要是由于：1)GPR探测深度较浅，增加了反演问题的不适定性；2)小偏移距测量模式导致GPR的低波数覆盖，使得GPR数据的FWI对初始模型的精度要求更高；3)介电常数与电导率在数量级上的巨大差异，使得它们具有不同的敏感度，给双参数反演带来了极大的挑战。一般而言，GPR数据的FWI对于相位信息(介电常数)更加敏感，因此当FWI反演过程中对相位信息进行拟合时会造成对振幅信息(电导率)的过拟合现象(overshoot)。如何有效提高GPR数据的FWI的精度，是目前急需解决的问题。
技术实现思路
针对现有技术存在的缺陷，本专利技术提供一种探地雷达数据时频组合全波形反演方法，可有效解决上述问题。本专利技术采用的技术方案如下：本专利技术提供一种探地雷达数据时频组合全波形反演方法，包括以下步骤：步骤1，获取研究区域的GPR原始观测数据；对所述GPR原始观测数据进行预处理，得到预处理后的GPR观测数据；步骤2，从预处理后的GPR观测数据中提取发射信号子波波形；根据预处理后的GPR观测数据，获得研究区域的介电常数模型和研究区域的电导率模型；其中，所述研究区域的介电常数模型是指：将研究区域格网化后，每个格网的位置坐标与格网的介电常数构成位置-介电常数数据对；研究区域的介电常数模型由研究区域的各个格网的位置-介电常数数据对组成；初始时，每个格网的介电常数通过对所述预处理后的GPR观测数据进行初级反演方法获得；所述研究区域的电导率模型是指：将研究区域格网化后，每个格网的位置坐标与格网的电导率构成位置-电导率数据对；研究区域的电导率模型由研究区域的各个格网的位置-电导率数据对组成；初始时，每个格网的电导率通过对所述预处理后的GPR观测数据进行初级反演方法获得；步骤3，设定双参数反演的反演频率，按从低频向高频方向，分别为：f1,f2,...,fm，f1,f2,...,fm分别为：第1个反演频率,第2个反演频率,...,第m个反演频率；设定每个反演频率的迭代次数均相等，为n；其中，反演频率f1,f2,...,fm通过以下方式确定：分析获得GPR原始观测数据包含的电磁波频率中的频率最大值和频率最小值；在频率最小值和频率最大值之间，按低频时反演频率采样密，随着频率的升高，反演频率采样逐渐变稀疏的原则，在频率最小值和频率最大值范围内，选取反演频率f1,f2,...,fm；对于每个反演频率，均构造一个对应的低通滤波器；因此，构造与反演频率f1,f2,...,fm对应的低通滤波器filter1,filter2,...,filterm；对于步骤1得到的预处理后的GPR观测数据，采用低通滤波器filter1进行滤波，生成对应的时间域观测数据data1；采用低通滤波器filter2进行滤波，生成对应的时间域观测数据data2；依此类推，采用低通滤波器filterm进行滤波，生成对应的时间域观测数据datam；对于步骤1得到的预处理后的GPR观测数据，进行傅立叶变换，得到与反演频率f1,f2,...,fm对应的频率域观测数据frequency1,frequency2,...,frequencym；步骤4，令i＝1；步骤5，当前采用的反演频率为反演频率fi；步骤6，频率域反演过程，方法为：步骤6.1，构建频率域目标函数S1(ε,σ)为：其中：ε是介电常数；σ是电导率；Ns是GPR发射点数量；Nr是每个发射点对应的GPR接收点数量；s表示发射点，以每个发射点为一组，一组内包含的接收点数量为r；表示模拟的频率域电场强度；表示观测的频率域电场强度；表示频率域介电常数扰动下的拟合残差；T表示转置矩阵；步骤6.2，电导率σ保持不变，根据介电常数模型，计算与反演频率fi对应的模拟的频率域电场强度根据频率域观测数据frequencyi，计算与反演频率fi对应的将和代入步骤6.1建立的频率域目标函数S1(ε,σ)，得到频率域目标函数值S1；步骤6.3，基于伴随矩阵方法，计算频率域目标函数值S1对介电常数的梯度方向；步骤6.4，根据抛物线原理，计算介电常数最优步长Δεi；步骤7，时间域反演过程，方法为：步骤7.1，构建时间域目标函数S2(ε,σ)为：其中：表示模拟的时间域电场强度；表示观测的时间域电场强度；表示时间域介电常数扰动下的拟合残差；δ表示时间变化系数；步骤7.2，介电常数ε保持不变，根据电导率模型，计算与反演频率fi对应的模拟的时间域电场强度根据时间域观测数据datai，计算与反演频率fi对应的将和代入步骤7.1建立的时间域目标函数S2(ε,σ)，得到时间域目标函数值S2；步骤7.3，基于伴随矩阵方法，计算时间域目标函数值S2对电导率的梯度方向；步骤7.4，根据抛物线原理，计算电导率最优步长Δσi；步骤8，更新双参数模型，即：更新介电常数模型和电导率模型，具体方法为：对于介电常数模型，其每个格网的介电常数，按步骤6.3确定的介电常数的梯度方向，增加介电常数最优步长Δεi，得到格网的更新后的介电常数，从而得到更新的介电常数模型；对于电导率模型，其每个格网的电导率，按步骤7.3确定的电导率的梯度方向，增加电导率最优步长Δσi，得到格网的更新后的电导率，从而得到更新的电导率模型；步骤9，判断是否达到迭代总结束条件；如果达到，则执行步骤11；如果没有达到，则执行步骤10；步骤10，判断当前反演频率fi是否已迭代n次；如果是，则令i＝i+1；返回步骤5，采用更新后的电导率模型和更新后的介电常数模型，进行下一个反演频率的全波形反演；如果否，则返回步骤5，采用更新后的电导本文档来自技高网...

【技术保护点】
1.一种探地雷达数据时频组合全波形反演方法，其特征在于，包括以下步骤：/n步骤1，获取研究区域的GPR原始观测数据；对所述GPR原始观测数据进行预处理，得到预处理后的GPR观测数据；/n步骤2，从预处理后的GPR观测数据中提取发射信号子波波形；/n根据预处理后的GPR观测数据，获得研究区域的介电常数模型和研究区域的电导率模型；/n其中，所述研究区域的介电常数模型是指：将研究区域格网化后，每个格网的位置坐标与格网的介电常数构成位置-介电常数数据对；研究区域的介电常数模型由研究区域的各个格网的位置-介电常数数据对组成；初始时，每个格网的介电常数通过对所述预处理后的GPR观测数据进行初级反演方法获得；/n所述研究区域的电导率模型是指：将研究区域格网化后，每个格网的位置坐标与格网的电导率构成位置-电导率数据对；研究区域的电导率模型由研究区域的各个格网的位置-电导率数据对组成；初始时，每个格网的电导率通过对所述预处理后的GPR观测数据进行初级反演方法获得；/n步骤3，设定双参数反演的反演频率，按从低频向高频方向，分别为：f

【技术特征摘要】
1.一种探地雷达数据时频组合全波形反演方法，其特征在于，包括以下步骤：
步骤1，获取研究区域的GPR原始观测数据；对所述GPR原始观测数据进行预处理，得到预处理后的GPR观测数据；
步骤2，从预处理后的GPR观测数据中提取发射信号子波波形；
根据预处理后的GPR观测数据，获得研究区域的介电常数模型和研究区域的电导率模型；
其中，所述研究区域的介电常数模型是指：将研究区域格网化后，每个格网的位置坐标与格网的介电常数构成位置-介电常数数据对；研究区域的介电常数模型由研究区域的各个格网的位置-介电常数数据对组成；初始时，每个格网的介电常数通过对所述预处理后的GPR观测数据进行初级反演方法获得；
所述研究区域的电导率模型是指：将研究区域格网化后，每个格网的位置坐标与格网的电导率构成位置-电导率数据对；研究区域的电导率模型由研究区域的各个格网的位置-电导率数据对组成；初始时，每个格网的电导率通过对所述预处理后的GPR观测数据进行初级反演方法获得；
步骤3，设定双参数反演的反演频率，按从低频向高频方向，分别为：f1,f2,...,fm，f1,f2,...,fm分别为：第1个反演频率,第2个反演频率,...,第m个反演频率；设定每个反演频率的迭代次数均相等，为n；
其中，反演频率f1,f2,...,fm通过以下方式确定：分析获得GPR原始观测数据包含的电磁波频率中的频率最大值和频率最小值；在频率最小值和频率最大值之间，按低频时反演频率采样密，随着频率的升高，反演频率采样逐渐变稀疏的原则，在频率最小值和频率最大值范围内，选取反演频率f1,f2,...,fm；
对于每个反演频率，均构造一个对应的低通滤波器；因此，构造与反演频率f1,f2,...,fm对应的低通滤波器filter1,filter2,...,filterm；
对于步骤1得到的预处理后的GPR观测数据，采用低通滤波器filter1进行滤波，生成对应的时间域观测数据data1；采用低通滤波器filter2进行滤波，生成对应的时间域观测数据data2；依此类推，采用低通滤波器filterm进行滤波，生成对应的时间域观测数据datam；
对于步骤1得到的预处理后的GPR观测数据，进行傅立叶变换，得到与反演频率f1,f2,...,fm对应的频率域观测数据frequency1,frequency2,...,frequencym；
步骤4，令i＝1；
步骤5，当前采用的反演频率为反演频率fi；
步骤6，频率域反演过程，方法为：
步骤6.1，构建频率域目标函数S1(ε,σ)为：



其中：
ε是介电常数；
σ是电导率；
Ns是GPR发射点数量；
Nr是每个发射点对应的GPR接收点数量；
s表示发射点，以每个发射点为一组，一组内包含的接收点数量为r；

表示模拟的频率域电场强度；

表示观...

【专利技术属性】
技术研发人员：杨曦，
申请(专利权)人：云南航天工程物探检测股份有限公司，
类型：发明
国别省市：云南;53