本发明涉及钻孔雷达目标成像技术领域，具体涉及一种基于非均匀频谱估计的钻孔雷达成像方法。
背景技术：
作为一种新型的深地层地球物理探测技术，钻孔雷达能够感知地下几千米深处的目标(如裂缝和空洞等异常地质结构)，获取目标的距离、方位、形状及电磁参数等信息。钻孔雷达已广泛应用于空洞和裂缝探测，水文地质学研究，温室气体跟踪，矿藏勘探，盐丘调查和测井等诸多领域，受到了全世界的广泛关注。
钻孔雷达在沿井孔运动过程中不断接收来自井孔周围介质的回波，并在后期通过雷达成像等技术获取井周介质的分布信息。在钻孔雷达的工作过程中，基于时间驱动的雷达系统在沿井孔运动时不可避免地会与井壁发生摩擦，造成系统沿井孔方向空间采样不均匀，进而对目标空间的高效率成像产生不利影响。
目前，钻孔雷达成像都是利用现有的雷达成像算法实现，包括后向投影算法，Kirchhoff偏移算法，逆时偏移算法和Stolt偏移算法等经典成像算法。然而，现有的雷达成像算法不能满足钻孔雷达非均匀采样条件下的高效率成像需求。
基于几何绕射理论，后向投影算法能够在时域中修正来自不同探测位置处的雷达回波，使得来自同一目标的回波能量能够在真实目标空间位置积累。然而，后向投影算法的计算复杂度会随着采样规模的增加而迅速增加。Kirchhoff偏移算法是基于波动方程的Huygen原理表达。与后向投影算法相比，Kirchhoff偏移算法能够取得更加精确的目标空间成像结果，但计算复杂度更高。逆时偏移算法能够有效地解决复杂介质环境中的目标成像问题，但该算法严重依赖于探测环境的先验信息，且计算复杂度很高。后向投影算法，Kirchhoff偏移算法和逆时偏移算法虽然都可以处理非均匀采样数据获取目标空间图像，但成像效率很低。
Stolt偏移算法可以借助于快速傅里叶变换及快速傅里叶逆变换等运算大大提高成像速度，在处理均匀采样数据时具有很高的计算效率。然而对于非均匀空间采样的钻孔雷达数据，Stolt偏移算法只能利用离散傅里叶变换来获取采样数据的频率波数谱，其计算复杂度很高，成像效率降低。
现有的雷达成像算法在钻孔雷达实际工程应用中的性能较差，特别是在处理非均匀空间采样数据时成像效率低下，制约了钻孔雷达在实际工程中的应用。
技术实现要素：
本发明所要解决的技术问题是满足钻孔雷达非均匀采样条件下的高效率成像需求，目的在于提供一种基于非均匀频谱估计的钻孔雷达成像方法，在保证钻孔雷达成像质量的前提下提高非均匀采样数据的高效率精确成像，提升钻孔雷达的实际应用价值。
本发明通过下述技术方案实现：
一种基于非均匀频谱估计的钻孔雷达成像方法，包括如下步骤：
步骤A、初始化钻孔雷达采样数据并建立位移数据向量表；
步骤B、利用频谱估计算法快速获取非均匀空间采样数据的频谱；
步骤C、对上一步骤中得到的采样数据的估计频谱进行插值运算获取目标空间频谱；
步骤D、得到目标空间图像。
本发明相比于传统后向投影算法，Kirchhoff偏移算法和逆时偏移算法，基于波动方程在频率波数域利用插值运算将非均匀采样数据矩阵转换为目标空间图像，在计算非均匀采样数据矩阵频谱时利用了快速傅里叶变换等运算，具有更高的运算效率；与传统Stolt偏移算法相比，将实际的系统位移信息考虑在内，在处理非均匀采样数据时具有更高的成像精度。
本发明在对钻孔雷达采样数据不做特殊要求的情况下，利用频谱估计算法快速获取非均匀空间采样数据的频谱，实现对钻孔雷达非均匀采样数据的高效率精确成像，该方法更加贴近钻孔雷达工作的实际情况，避免了在处理非均匀采样数据时带来的大规模计算量，可以在成像分辨率和算法复杂度之间取得很好的平衡，算法实用性强且成像性能稳定，大大提高了对非均匀钻孔雷达回波数据的成像效率，实现钻孔雷达的高效精确成像。同时本发明的思想和方法也可应用其他类似的雷达成像场景。
优选的，步骤A中初始化钻孔雷达的采样数据包括如下步骤：
A1.利用钻孔雷达系统对井周介质的雷达回波进行采样，建立若干维度的非均匀采样数据矩阵u，u的表达式如下
利用滑动滤波去除非均匀采样数据矩阵中的直达波信号及其他串扰信号，同时利用定位系统测量井下系统的位置并记录位移信息，建立相应维度的位移数据向量表[z x]，包括均匀采样维度集z和非均匀采样维度集x，其中z和x均为单一维度集，z＝[t1，t2，…，tL]，x＝[x1，x2，…，xK]，L和K为正整数；
A2.将钻孔雷达的采样数据矩阵u，在均匀采样维度集z和非均匀采样维度集x上补零至必要长度。
优选的，步骤B中将采样数据利用频谱估计算法快速获取非均匀空间采样数据的频谱的具体步骤如下：
B1.建立采样数据矩阵u的频率波数谱表达U，其中
为矩阵U的元素，Δxn为x维度上的实际采样间隔，Δx为x维度上的平均采样间隔，n′＝-N/2+1,…,N/2，l′＝-L/2+1,…,L/2，在采样数据矩阵的所有均匀采样维度集z上进行快速傅里叶变换，得到均匀采样维度下的部分频谱表达
B2.综合考虑采样数据矩阵频率波数谱的估计误差及计算复杂度，选取合适的过采样系数m和估算指数q；
B3.选取某一非均匀采样维度xi(i＝1，2，…，K)，将部分频谱表达U1在该维度下的每个非均匀指数因子都用若干加权的均匀指数因子近似估计其中ε为最大估计误差，f(m,n′)为加权系数，并根据最小二乘准则计算该非均匀采样维度下的加权因子向量ρ；
B4.利用步骤B3中的加权估计关系将该非均匀采样维度均匀化，并利用U1构造均匀的虚拟采样集合
B5.对该非均匀采样维度下的虚拟采样集合τ进行快速傅里叶变换，并选取[-π,π]内的结果按加权处理，加权结果作为该非均匀采样维度下的部分频谱估计表达
B6.将估计频谱U2赋给U1，重复步骤B3-B5，直至所有非均匀采样维度下的频谱估算完毕，将最终估算得到的频谱U2作为频率波数谱表达U的值，记为U’。通过该方法即可快速获取非均匀空间采样数据的频谱。
优选的，步骤B2中过采样系数m>1，估算指数q为大于等于8的偶数。
优选的，过采样系数m和估算指数q的取值范围为1.4<m<4，8≤q≤16。
优选的，步骤C中对上一步骤中得到的采样数据的估计频谱进行插值运算获取目标空间频谱的具体内容为确定钻孔雷达发射信号在介质中的波速及系统的时空采样间隔，初始化采样波数，由波数频率约束关系对采样数据矩阵的估计频谱U’做角频率插值运算，生成目标空间的频率波数谱。
优选的，步骤D中通过对目标空间的频率波数谱做快速傅里叶逆变换得到目标空间图像。
本发明与现有技术相比，具有如下的优点和有益效果：
1、本发明相比于传统后向投影算法，Kirchhoff偏移算法和逆时偏移算法，基于波动方程在频率波数域利用插值运算将非均匀采样数据矩阵转换为目标空间图像，在计算非均匀采样数据矩阵频谱时利用了快速傅里叶变换等运算，具有更高的运算效率；与传统Stolt偏移算法相比，将实际的系统位移信息考虑在内，在处理非均匀采样数据时具有更高的成像精度。
2、本发明在对钻孔雷达采样数据不做特殊要求的情况下，利用频谱估计算法快速获取非均匀空间采样数据的频谱，实现对钻孔雷达非均匀采样数据的高效率精确成像，该方法更加贴近钻孔雷达工作的实际情况，避免了在处理非均匀采样数据时带来的大规模计算量，可以在成像分辨率和算法复杂度之间取得很好的平衡，算法实用性强且成像性能稳定，大大提高了对非均匀钻孔雷达回波数据的成像效率，实现钻孔雷达的高效精确成像。
附图说明
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：
图1为本发明的总体流程图；
图2为仿真场景示意图；
图3为钻孔雷达系统二维非均匀采样数据成像图；
图4为基于图3中的仿真数据，本发明所提成像方法、Stolt偏移方法、后向投影方法和Kirchhoff偏移方法所得到的成像结果；
图5为本发明所提成像方法、后向投影方法和Kirchhoff偏移方法在不同非均匀采样规模下成像所需时间的仿真结果。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。
实施例1：
如图1所示，本发明包括一种基于非均匀频谱估计的钻孔雷达成像方法，包括如下步骤：
步骤A、初始化钻孔雷达采样数据并建立位移数据向量表；
步骤B、利用频谱估计算法快速获取非均匀空间采样数据的频谱；
步骤C、对上一步骤中得到的采样数据的估计频谱进行插值运算获取目标空间频谱；
步骤D、得到目标空间图像。
由于在钻孔雷达的工作过程中，基于时间驱动的雷达系统在沿井孔运动时不可避免地会与井壁发生摩擦，造成系统沿井孔方向空间采样不均匀，进而对目标空间的高效率成像产生不利影响。目前，钻孔雷达成像都是利用现有的雷达成像算法实现，包括后向投影算法，Kirchhoff偏移算法，逆时偏移算法和Stolt偏移算法等经典成像算法。然而，现有的雷达成像算法不能满足钻孔雷达非均匀采样条件下的高效率成像需求。基于几何绕射理论，后向投影算法能够在时域中修正来自不同探测位置处的雷达回波，使得来自同一目标的回波能量能够在真实目标空间位置积累。然而，后向投影算法的计算复杂度会随着采样规模的增加而迅速增加。Kirchhoff偏移算法是基于波动方程的Huygen原理表达。与后向投影算法相比，Kirchhoff偏移算法能够取得更加精确的目标空间成像结果，但计算复杂度更高。逆时偏移算法能够有效地解决复杂介质环境中的目标成像问题，但该算法严重依赖于探测环境的先验信息，且计算复杂度很高。后向投影算法，Kirchhoff偏移算法和逆时偏移算法虽然都可以处理非均匀采样数据获取目标空间图像，但成像效率很低。Stolt偏移算法可以借助于快速傅里叶变换及快速傅里叶逆变换等运算大大提高成像速度，在处理均匀采样数据时具有很高的计算效率。然而对于非均匀空间采样的钻孔雷达数据，Stolt偏移算法只能利用离散傅里叶变换来获取采样数据的频率波数谱，其计算复杂度很高，成像效率降低。现有的雷达成像算法在钻孔雷达实际工程应用中的性能较差，特别是在处理非均匀空间采样数据时成像效率低下，制约了钻孔雷达在实际工程中的应用。
本发明相比于传统后向投影算法，Kirchhoff偏移算法和逆时偏移算法，基于波动方程在频率波数域利用插值运算将非均匀采样数据矩阵转换为目标空间图像，在计算非均匀采样数据矩阵频谱时利用了快速傅里叶变换等运算，具有更高的运算效率；与传统Stolt偏移算法相比，将实际的系统位移信息考虑在内，在处理非均匀采样数据时具有更高的成像精度。
本发明在对钻孔雷达采样数据不做特殊要求的情况下，利用频谱估计算法快速获取非均匀空间采样数据的频谱，实现对钻孔雷达非均匀采样数据的高效率精确成像，该方法更加贴近钻孔雷达工作的实际情况，避免了在处理非均匀采样数据时带来的大规模计算量，可以在成像分辨率和算法复杂度之间取得很好的平衡，算法实用性强且成像性能稳定，大大提高了对非均匀钻孔雷达回波数据的成像效率，实现钻孔雷达的高效精确成像。同时本发明的思想和方法也可应用其他类似的雷达成像场景。
实施例2：
本实施例在实施例1的基础上以二维钻孔雷达探测为例说明具体实施步骤，优选如下：步骤A中获取钻孔雷达的采样数据包括如下步骤：
A1.利用钻孔雷达系统对井周介质的雷达回波进行采样，建立若干维度的非均匀采样数据矩阵u，u的表达式如下
利用滑动滤波去除非均匀采样数据矩阵中的直达波信号及其他串扰信号，同时利用定位系统测量井下系统的位置并记录位移信息，建立相应维度的位移数据向量表[z x]，包括均匀采样维度集z和非均匀采样维度集x，其中z和x均为单一维度集，z＝[t1，t2，…，tL]，x＝[x1，x2，…，xK]，L和K为正整数；
A2.将钻孔雷达的采样数据矩阵u，在均匀采样维度集z和非均匀采样维度集x上补零至必要长度。该方法为快速获取钻孔雷达的采样数据包的一种具体实现方式，通过该方法能够快速得到不做特殊处理后的原始数据，供后面步骤使用。
实施例3：
本实施例在上述实施例的基础上优选如下：步骤B中将采样数据利用频谱估计算法快速获取非均匀空间采样数据的频谱的具体步骤如下：
B1.建立采样数据矩阵u的频率波数谱表达U，其中
为矩阵U的元素，Δxn为x维度上的实际采样间隔，Δx为x维度上的平均采样间隔，n′＝-N/2+1,…,N/2，l′＝-L/2+1,…,L/2，在采样数据矩阵的所有均匀采样维度集z上进行快速傅里叶变换，得到均匀采样维度下的部分频谱表达
B2.综合考虑采样数据矩阵频率波数谱的估计误差及计算复杂度，选取合适的过采样系数m和估算指数q；
B3.选取某一非均匀采样维度xi(i＝1，2，…，K)，将部分频谱表达U1在该维度下的每个非均匀指数因子都用若干加权的均匀指数因子近似估计其中ε为最大估计误差，f(m,n′)为加权系数，并根据最小二乘准则计算该非均匀采样维度下的加权因子向量ρ；
B4.利用步骤B3中的加权估计关系将该非均匀采样维度均匀化，并利用U1构造均匀的虚拟采样集合
B5.对该非均匀采样维度下的虚拟采样集合τ进行快速傅里叶变换，并选取[-π,π]内的结果按加权处理，加权结果作为该非均匀采样维度下的部分频谱估计表达
B6.将估计频谱U2赋给U1，重复步骤B3-B5，直至所有非均匀采样维度下的频谱估算完毕，将最终估算得到的频谱U2作为频率波数谱表达U的值，记为U’。通过该方法即可快速获取非均匀空间采样数据的频谱。
实施例4：
本实施例在上述实施例的基础上优选如下：步骤B2中过采样系数m>1，估算指数q为大于等于8的偶数。
过采样系数m和估算指数q的取值范围为1.4<m<4，8≤q≤16。
步骤C中对上一步骤中得到的采样数据的估计频谱进行插值运算获取目标空间频谱的具体内容为确定钻孔雷达发射信号在介质中的波速及系统的时空采样间隔，初始化采样波数，由波数频率约束关系对采样数据矩阵的估计频谱U’做角频率插值运算，生成目标空间的频率波数谱。
步骤D中通过对目标空间的频率波数谱做快速傅里叶逆变换得到目标空间图像。
该方法更加贴近钻孔雷达工作的实际情况，进一步避免了在处理非均匀采样数据时带来的大规模计算量，可以在成像分辨率和算法复杂度之间取得很好的平衡，算法实用性强且成像性能稳定，大大提高了对非均匀钻孔雷达回波数据的成像效率，实现钻孔雷达的高效精确成像。
本发明的效果通过以下仿真对比试验进一步说明：
如图2所示，考虑二维单孔钻孔雷达探测场景，钻孔直径为0.2m，系统水平最大探测距离为2m，垂直最大探测深度为8m。井周介质为均匀的花岗岩，其相对介电常数为9。使用单收单发双基钻孔雷达系统对井周介质中直径为0.3m的空洞进行探测，系统收发天线均视为理想点源天线。钻孔雷达收发天线之间的间隔为0.1m。仿真使用中心频率为300MHz的一阶Balckman-Harris脉冲作为发射信号，发射信号的最小波长为0.17m。系统的时间采样间隔约为0.1ns，沿井孔方向空间采样不均匀。接收天线在200个非均匀探测深度上采集了来自0.5m至6.5m井深的回波信号。在每个探测深度上系统采样350次。得到系统采样数据如图3所示。
如图4所示成像结果，是基于图3中的仿真数据，采用本发明所提成像方法、Stolt偏移方法、后向投影方法和Kirchhoff偏移方法所得到的成像结果。本发明所提成像方法使用的过采样系数m和估算指数q的值分别为m＝2和q＝8。由图4中各成像方法结果对比可知，本发明所提成像方法的成像结果中目标能量聚焦效果最好，目标定位准确，成像结果中的旁瓣水平为-24.48dB。Stolt偏移方法、后向投影方法和Kirchhoff偏移方法成像结果的旁瓣水平分别为-19.96dB、-22.53dB和-23.11dB，均高于本发明所提成像方法，其中Stolt偏移方法的成像结果中目标定位有明显偏差，成像结果有误。此外，在本次仿真中，本发明所提成像方法计算时间为0.2s，后向投影方法和Kirchhoff偏移方法计算时间分别为4s和5.6s，均远高于本发明所提成像方法，其成像效率均低于本发明所提成像方法。
本发明所提成像方法、后向投影方法和Kirchhoff偏移方法在不同非均匀采样规模下成像所需时间仿真结果如图5所示。由图5可以看出，随着采样规模的增大，后向投影方法和Kirchhoff偏移方法成像所需的时间迅速增加，本发明所提成像方法的高成像效率优势更加明显。
通过本发明的具体实施可以看出，在钻孔雷达非均匀采样场景中，本发明提供的成像方法更加适合于钻孔雷达工作的实际情况，避免了在处理非均匀采样数据时带来的大规模计算量，可以在成像分辨率和算法复杂度之间取得很好的平衡，算法实用性强且成像性能稳定，大大提高了对非均匀钻孔雷达回波数据的成像效率，实现钻孔雷达的高效精确成像。
本发明未详细说明部分，为本领域内人员公知常识。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。