一种平面三层介质中电磁波传播时延的计算方法

摘要

本发明公开了一种平面三层介质中电磁波传播时延的计算方法，通过双曲线波前近似理论获得电磁波在第一二层介质中由天线位置到二三层介质交界面上某点C的传播时延。随着点C在二三层介质交界面上的移动，总的单程传播时延也随之变化。即可通过一维搜索或数学解析的方法获得传播时延的最小值作为真实传播时延的估计值，该方法避免了传统的二维搜索或求解二元四次方程组的计算方法，大大降低了运算量，保证了传播时延的计算精度，适用于探地雷达实时成像处理。

说明书

技术领域

本发明属于探地雷达成像技术领域，涉及一种平面三层介质中电磁波传播时 延的计算方法。

背景技术

探地雷达是一种有效的无损探测技术。它通过空域扫描向探测区域发射电 磁波并接收散射回波，可实现对未知区域内部的成像处理，获得未知区域中的隐 蔽目标参数，即目标分布信息和散射强度信息，有效应用于市政工程、考古、地 雷探测、反恐等多种场合。

探地雷达的成像方法有多种，基于“延时-累加”处理的后向投影成像算 法适用于非等间距采样下对分层介质中点散射型目标的成像处理，因而广泛应用 于探地雷达信号处理中。后向投影成像算法的关键处理步骤是传播时延的计算。 对成像区域中的每一网格位置(x_o，z_o)，根据探测扫描场景分别计算各测点对应 的发射天线x_T和接收天线x_R到成像点(x_o，z_o)的电磁波单程传播时延。本专利适 用于平面三层介质环境下第一层介质(空气)中各天线测点到第三层介质中各成 像点的单程传播时延的计算。

对三层平面介质结构的扫描区域而言，探地雷达天线位于第一层介质，即 空气中，用于对整个探测区域进行扫描。整个成像区域可分为两部分，即第二层 介质部分和第三层介质部分。成像处理便是获得第二层介质中各点和第三层介质 中成像区域覆盖的各点的散射强度。在进行后向投影成像运算时，对第二层介质 中的各点，需要计算电磁波经过一次折射后的传播时延，常用的方法有搜索法、 线性近似法【参考文献：E.M.Johansson，J.E.Mast，Three dimensional ground  penetrating radar imaging using a synthetic aperture time-domain focusing， in：Proceedings of the SPIE Conference on Advanced Microwave and  Millimeter Wave Detectors，2275，1994，pp.205-214.】和双曲线波前近似 法【参考文献：Rappaport，C.，A novel，non-iterative，analytic method to  find the surface refraction point for air-coupled ground penetrating  radar，Proceedings of the 5th European Conference on Antennas and  Propagation(EUCAP)，1786-1789，2011】。对第三层介质中成像区域的各点， 需要计算各成像点到第一层介质(空气)中探地雷达天线各测点的传播时延。根 据电磁波传播规律，电磁波在探地雷达天线和第三层介质中成像点间的传播要经 过两次折射，即在第一二层介质交界面处产生一次折射，在第二三层介质交界面 处产生第二次折射，如图1所示。

传统的时延计算方法是通过两维搜索获得单程传播时延的极小值，即

$\left(\begin{array}{c}\min \tau =\frac{\sqrt{{ϵ}_1}·\sqrt{z_a^2+{(x_a-x_{r1})}^2}+\sqrt{{ϵ}_2}·\sqrt{d^2+{(x_{r2}-x_{r1})}^2}+\sqrt{{ϵ}_3}·\sqrt{{(d-z_o)}^2+{(x_o-x_{r2})}^2}}{c}\\ x_{r1}\in [x_a,x_o]\\ x_{r2}\in [x_a,x_o]\end{array}\right)$

其中，x_r1和x_r2分别为第一二层介质交界面和第二三层介质交界面上的点，ε₁、 ε₂和ε₃分别为第一、二、三层介质的相对介电常数值，第一层介质中发射天线 (测点)的坐标为(x_a，z_a)，第三层介质中成像点的坐标为(x_o，z_o)，第二三层介 质交界面的纵向坐标为z＝d。

或者是求解二元四次方程组 $\left(\begin{array}{c}{ϵ}_1{(x_{r1}-x_a)}^2[d^2+{(x_{r2}-x_{r1})}^2]={ϵ}_2{(x_{r2}-x_{r1})}^2[z_a^2+{(x_a-x_{r1})}^2]\\ {ϵ}_2{(x_{r2}-x_{r1})}^2[{(d-z_o)}^2+{(x_o-x_{r2})}^2]={ϵ}_3{(x_o-x_{r2})}^2[d^2+{(x_{r2}-x_{r1})}^2]\end{array}\right)$获得两个折 射点的横向坐标x_r1和x_r2，再通过公式 $\tau =\frac{\sqrt{{ϵ}_1}·\sqrt{z_a^2+{(x_a-x_{r1})}^2}+\sqrt{{ϵ}_2}·\sqrt{d^2+{(x_{r2}-x_{r1})}^2}+\sqrt{{ϵ}_3}·\sqrt{{(d-z_o)}^2+{(x_o-x_{r2})}^2}}{c}$计 算获得单程传播时延。这种时延值的计算方法计算量很大，不适用于探地雷达实 时成像处理。

发明内容

本发明所要解决的技术问题是提供一种平面三层介质中电磁波传播时延的 计算方法，该计算方法计算量小(本专利提出的计算方法的计算量要低于常规计 算方法的计算量)，能提高探地雷达的成像处理效率。

发明的技术解决方案如下：

一种平面三层介质中电磁波传播时延的计算方法，其特征在于，包括以下步 骤：

首先设定第二层与第三层介质交界面处折射点的估计区间[x_a，x_o]，其中x_a和 x_o分别为第一层介质中天线测点A的横向坐标和第一层介质中天线测点A在第 三层介质中的成像点D的横向坐标；

然后在估计区间中的每一点均对应一个传播时延值；

估计区间中任一点x_r′对应的传播时延值计算公式为 $\tau \left(x_r^{\prime }\right)=\frac{\sqrt{{ϵ}_2}{l}_1\left(x_r^{\prime }\right)+\sqrt{{ϵ}_3}\sqrt{{(x_o-x_r^{\prime })}^2+{(d-z_o)}^2}}{c},$其中ε₂和ε₃分别为第二介质和第 三层介质的相对介电常数值，x_o和z_o分别为第三层介质中成像点D的横向和纵 向坐标，x_r′为估计区间中任一点的横向坐标，d为第二介质和第三层介质交界 面的纵向坐标值，c＝3×10⁸m/s为光速， $l_1\left(x_r^{\prime }\right)=\frac{-b+\sqrt{b^2-4a[d^2-k^2{(x_r^{\prime }-x_a)}^2-{\left(\frac{z_a}{n_1}\right)}^2]}}{2a},$$a=2\sqrt{1+k^2}-1,$$b=2[(d-\frac{z_a}{n_1})\sqrt{1+k^2}+\frac{z_a}{n_1}],$$k=\frac{1}{\sqrt{n_1^2-1}},$$\sqrt{\frac{{ϵ}_2}{{ϵ}_1}}=n_1,$ε₁为第一层介质的相对介 电常数值，z_a为第一层介质中天线的纵向坐标；

最后求取最小的传播时延值作为真实电磁波传播时延的估计值。

一种方法为：通过一维搜索方法选取最小的传播时延值作为真实电磁波传播 时延的估计值。

另一种方法为：通过数学解析方法求解最小的传播时延值作为真实电磁波传 播时延的估计值；

所述的数学解析方法为：求得所述传播时延值计算公式的一阶导数，令该一 阶导数为零，求出对应的x_r′值，然后再将x_r′代入传播时延值计算公式计算得 到的传播时延值即为最小的传播时延值。

所述的一维搜索方法为等距搜索或最速下降法。

所述的等距搜索方法如下：(根据传播时延的计算精度要求)将第二层与第三 层介质交界面处折射点的估计区间[x_a，x_o]等分为N份，该N份对应于N+1个点； 在每个点处，都计算获得一个传播时延值，从而获得N+1个传播时延值，再通 过逐一比较的方法获得这N+1个传播时延的最小值作为真实传播时延的估计 值。

本发明通过双曲线波前近似理论获得电磁波在第一二层介质中由天线位置 到二三层介质交界面上某点C的传播时延。随着点C在二三层介质交界面上的 移动，总的单程传播时延也随之变化。即可通过一维搜索或数学解析的方法获得 传播时延的最小值作为真实传播时延的估计值，该方法避免了传统的二维搜索或 求解二元四次方程组的计算方法，大大降低了运算量，保证了传播时延的计算精 度，适用于探地雷达实时成像处理。

公式推导过程如下：

总的传播时延可表示为第一层介质中源点到第二三层介质交界面处某点的 传播时延和第二三层介质交界面处某点到第三层介质中成像点时延的累加形式。 下面首先计算第一层介质中源点到第二三层介质交界面处某点的传播时延。根据 双曲线波前近似理论【参考文献：Rappaport，C.，A novel，non-iterative， analytic method to find the surface refraction point for air-coupled  ground penetrating radar，Proceedings of the 5th European Conference on  Antennas and Propagation(EUCAP)，1786-1789，2011】，当源点位于第一 层介质中的某点(x_a，z_a)向外辐射电磁波时，若假定第一层介质与第二层介质相 同，则辐射场的波前可等效为双曲线形式。波前的双曲线方程为 其中ε₁和ε₂分别为第一、二层介质的相对 介电常数值。根据惠更斯原理，该双曲线波前上的任一点可以作为二次波源向外 辐射电磁波，辐射方向沿该点的法线方向。根据费马原理，相同的时间内，波前 上任一点沿各自法线方向的传播距离是相同的，这就又形成了对应时刻的波前。 考虑到波前是向无穷远处延伸，无穷远处两个波前的距离应该是不变的。因此， 两个双曲线波前的渐近线应该是平行的。根据以上三点准则，即可获得任一时刻 (或任一传播距离)时的双曲线波前形式。

设该双曲线经过l个距离的平移后，过第二三层介质交界面上的折射点 (x′_r，d)。渐近线的顶点平移后为双曲线的顶点平移 后为从而，新的双曲线波前方程为 ${(z+l\sqrt{1+k^2})}^2-k^2{(x-x_a)}^2={(\frac{z_a}{n_1}-l+l\sqrt{1+k^2})}^2.$该双曲线过点(x′_r，d)，从而有 ${(d+l\sqrt{1+k^2})}^2-k^2{(x_r^{\prime }-x_a)}^2={(\frac{z_a}{n_1}-l+l\sqrt{1+k^2})}^2.$整理后可得 $l^2(2\sqrt{1+k^2}-1)+2l[(d-\frac{z_a}{n_1})\sqrt{1+k^2}+\frac{z_a}{n_1}]+d^2-k^2{(x_r^{\prime }-x_a)}^2-{\left(\frac{z_a}{n_1}\right)}^2=0.$该方程为l 的一元二次方程，求解获得其中的一个根为其中a、b、 c分别为该一元二次方程的二次项系数、一次项系数和常数项。则第一层介质中 源点(x_a，z_a)到第二三层介质交界面处某点(x′_r，d)的传播距离为第二 三层介质交界面处折射点(x′_r，d)到第三层介质中成像点(x_o，z_o)的传播路径为直 线，传播距离为从而，第一层介质中源点(x_a，z_a)到 第三层介质中成像点(x_o，z_o)的总传播时延可表示为 $\tau \left(x_r^{\prime }\right)=\frac{\sqrt{{ϵ}_2}{l}_1\left(x_r^{\prime }\right)+\sqrt{{ϵ}_3}\sqrt{{(x_o-x_r^{\prime })}^2+{(d-z_o)}^2}}{c},$其中ε₂和ε₃分别为第二三层介质 的相对介电常数值，x_o和z_o为第三层介质中成像点的横向和纵向坐标，x_r′为估 计区间中任一点的横向坐标，d为二三层介质交界面的纵向坐标值， c＝3×10⁸m/s为光速，$l_1\left(x_r^{\prime }\right)=\frac{-b+\sqrt{b^2-4a[d^2-k^2{(x_r^{\prime }-x_a)}^2-{\left(\frac{z_a}{n_1}\right)}^2]}}{2a}.$$a=2\sqrt{1+k^2}-1,$$b=2[(d-\frac{z_a}{n_1})\sqrt{1+k^2}+\frac{z_a}{n_1}],$$k=\frac{1}{\sqrt{n_1^2-1}},$$\sqrt{\frac{{ϵ}_2}{{ϵ}_1}}=n_1,$ε₁第一 层介质的相对介电常数值，z_a为第一层介质中天线的纵向坐标。

最后对整个估计区间中获得的所有时延值进行一维搜索，获得最小时延值 作为真实时延值的估计。

或者通过解析方法(即解方程的方法)求取最小时延值；

两种方法的计算结果是相同的，计算所需时间上略有差异而已。

有益效果：

本发明的平面三层介质中电磁波传播时延的计算方法，基于双曲线波前近似 方法，通过一维搜索或解析求解即可获得三层平面介质中的电磁波传播时延，避 免了传统的二维搜索或求解二元四次方程组的计算方法，大大降低了运算量，保 证了传播时延的计算精度，适用于探地雷达实时成像处理。

附图说明

图1示出了三层平面介质中电磁波传播路径示意图。

图2示出了第三层介质中各成像点到源点的真实单程传播时延。具体参数如 下：ε₁＝ε₀，ε₂＝8ε₀，ε₃＝12ε₀，ε₀＝8.85419×10^-12F/m为常数，表示了真空 中的介电常数值，d＝-0.5m，第一层介质中源点(天线)的坐标为$\left(\begin{array}{c}{x}_a=1m\\ z_a=0.4m\end{array}\right).$第三层介质中的成像区域设置为$\left(\begin{array}{c}\mathrm{xo}\in [-\mathrm{1,3}]m\\ \mathrm{zo}\in [-0.6,1.6]m\end{array}\right).$

图3示出了运用本专利方法的单程时延计算结果，具体参数与图2计算的参 数相同。

图4示出了图2和图3对应的单程传播时延真实值和近似值的相对误差。

图5示出了第三层介质中各成像点到源点的真实单程传播时延。具体参数如 下：ε₁＝ε₀，ε₂＝8ε₀，ε₃＝4ε₀，d＝-0.8m，第一层介质中源点(天线)的坐标 为$\left(\begin{array}{c}{x}_a=1m\\ z_a=0.2m\end{array}\right).$第三层介质中的成像区域设置为$\left(\begin{array}{c}\mathrm{xo}\in [-\mathrm{1,3}]m\\ \mathrm{zo}\in [-\mathrm{1,3}]m\end{array}\right).$

图6示出了运用本专利方法的单程时延计算结果，具体参数与图5计算的参 数相同。

图7示出了图5和图6对应的单程传播时延真实值和近似值的相对误差。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

本实例中，ε₁＝ε₀，ε₂＝8ε₀，ε₃＝12ε₀，ε₀为常数，表示了真空中的介电 常数值，d＝-0.5m，第一层介质中源点(天线测点)的坐标为$\left(\begin{array}{c}{x}_a=1m\\ z_a=0.4m\end{array}\right).$第 三层介质中的成像区域设置为$\left(\begin{array}{c}\mathrm{xo}\in [-\mathrm{1,3}]m\\ \mathrm{zo}\in [-0.6,1.6]m\end{array}\right).$电磁波从该成像区域中任一点 到源点的传播都是经历了两次折射。本实例计算电磁波在第一层介质中固定的源 点(天线测点)与第三层介质中各成像点间的单程传播时延。首先运用传统的处 理方法，通过二维寻优，计算获得成像区域中各点对应的单程传播时延如图2 所示，记为τ_real，这里是很多时延值，而且是真实值。此处计算真实值的目的是 与本专利的计算结果相对比。第三层介质中的二维区间中每个点都对应一个传播 时延。此处的传播时延是真实值，每个点的传播时延都是通过大运算量的二维寻 优得到的。在CPU为Pentium(R)Dual-Core 2.00GHz，内存为2.00GB的PC上 运行，所需的时间为7.3937秒。通过本专利所提的一维等距搜索方法计算获得 的成像区域中各点对应的单程传播时延如图3所示，记为τ_app。相同的PC机上， 计算所需的时间为0.3239秒。定义相对误差时延值为则二者的 相对误差如图4所示。从图4可见，二者的相对误差沿着x_a点的位置向左右两 侧增大，相对误差最大值为1.48％。在某些成像点上，相对误差为零。本专利所 提的处理方法的计算时间不到传统的时延计算方法处理时间的

在本实施例中，对第三层介质成像区域中的每一点，分别运用本专利所提的 一维等距搜索时延计算方法，获得了成像区域中各点对应的传播时延。

实施例2：

本实例中，ε₁＝ε₀，ε₂＝8ε₀，ε₃＝4ε₀，d＝-0.8m，第一层介质中源点(天 线测点)的坐标为$\left(\begin{array}{c}{x}_a=1m\\ z_a=0.2m\end{array}\right).$第三层介质中的成像区域设置为$\left(\begin{array}{c}\mathrm{xo}\in [-\mathrm{1,3}]m\\ \mathrm{zo}\in [-\mathrm{1,3}]m\end{array}\right).$传统的时延计算方法结果如图5所示，计算时间为8.4228秒。本专利所提的数学 解析方法的计算结果如图6所示，计算时间为0.3408秒。二者的相对误差如图7 所示。从图7可见，二者的相对误差沿着x_a点的位置向左右两侧增大，相对误差 最大值为1.47％，最小值为零。本专利所提的处理方法的计算时间不到传统的时 延计算方法处理时间的

在本实施例中，对第三层介质成像区域中的每一点，分别运用本专利所提的 数学解析时延计算方法，获得了成像区域中各点对应的传播时延。