CN105467222B - 基于单基地测量的地表介质参数反演方法 - Google Patents

Abstract

基于单基地测量的地表介质参数反演方法，涉及探地雷达。利用单天线垂直入射地面的测量结果对大地表层介质介电常数、电导率进行反演。本发明首先基于矢量网络分析仪测试获得单个喇叭天线垂直照射土壤表面(或金属板)的反射信号S₁₁，根据测试结果引入Matrix Pencil(矩阵束)算法，进行不同界面反射信号的辨识分离，估计提取出由土壤表面(或金属板表面)造成的真正的反射信号，从而得到土壤表面的真正的反射系数，利用该方法能有效将各个界面反射分量加以区分，滤除非土壤表面反射信号及其他噪声信号，使得估计精度得到提高，对快速准确获得土壤表层电磁参数具有重要的意义。

Description

基于单基地测量的地表介质参数反演方法

技术领域

本发明涉及探地雷达(GPR)，尤其是涉及一种基于单基地测量的地表介质参数反演方法。

背景技术

反演地表的介电参数是GPR(探地雷达)的重要应用之一，近二十多年来，地表参数的反演已发展为多波段、多极化、多角度地表数据的反演，地表介电常数和粗糙度的反演取得了显著的发展。目标的介电常数反映了介质的基本特性,是微波遥感所研究的的重要参数之一。它与物质的结构、组成、密度等因素有关。测量介电常数的方法主要有波导法、谐振腔法、空间波法、探针法等[1]。每种测量方法都有一定的适用范围和应用特点。其中波导法、谐振腔法等属于实验室测量方法,测量通常是在室中进行,要求具有相应的样品采集技术,由于环境改变(如密度、压力和湿度改变)和对样品的加工处理势必会造成测量与真实值之间的偏差。空间波法、探针法等属于实地测量方法,探针法测量实验相对稳定,但是只能代表局部特征。自由空间法测量不需要采集样品,可以直接进行测量。与其他测量方法相比,自由空间法对测试样品没有严格的形状和工艺要求,只需被测样品的厚度均匀,且具有一定大的面积以避免边沿绕射,可以在较宽的频率范围内进行扫频测量。

利用自由空间法对样品的电磁参数进行测量的误差主要由系统误差与人为误差造成。自由空间法属于开放式测量，因此测量过程中会有很多因素影响测量结果，其主要系统误差有：

1.喇叭天线远场处的平面波近似导致的误差；

2.当电磁波照射到样品时由样品边缘产生的多次绕射及散射效应引起的误差；

3.喇叭天线与被测样品之间的电磁波多次反射带来的误差；

4.被测样品位置参考面的选择带来的误差。

在现有的GPR信号参数反演中，面临的主要困难是难以排除多径效应和收发信号间的直接耦合对实验精度的影响。在非理想暗室环境中天线测试时，暗室的侧壁与后墙、设备支架等造成的多径效应成为影响测量精度的主要原因。若受到客观条件限制无法在暗室中进行，在外场更是极易受到背景噪声电磁波干扰，从而很大程度上影响实验精度[2]，这也影响探地雷达在实际中的运用，Greg Hislop提出了一种方法反演大地介电常数的方法，通过建立介电常数与频率的数学关系模型，在不需了解天线或测试系统的情况，准确反演出大地电磁参数[3]。

为了有效减少多径干扰，国内外均采用两种方法：逆傅里叶变换法和矩阵束法[4][5]。由于信道频率响应的测量广泛用于多路径传播渠道的脉冲响应或时间响应，因而在时间域中，通过测量多径效应，可以用来估测暗室的脉冲响应。逆傅里叶变换的应用可以解决上述问题。然而，逆傅里叶变换只是单方面从时间域或频率域处理测试数据，并不能全面地描述信号的特征。矩阵束方法对于解决多径干扰具有较好的效果。通过将接收器接收到处于不同频率和不同方位角的波，等效成不同的指数函数，矩阵束法更方便的找出直接由发射器传送到接收器的波。

针对抗噪声能力和估计精度的进一步要求,矩阵束方法(Matrix Pencil Method)被广泛研究[6][7]。其基本思想是根据数据构造两个特殊的Hankel数据矩阵,根据数据矩阵间的关系求解它们的广义特征值,广义特征值包含了所要求解的信息(信号的极点),因此,求解衰减指数信号极点问题转化为求解矩阵束的广义特征值问题。在求解广义特征值过程中,为抑制噪声干扰,引入SVD分解和矩阵的低秩近似等方法,使算法具有好的鲁棒性。

中国专利201410299068.4公开一种基于Matrix Pencil的电力系统低频振荡模态辨识方法，该方法针对实际系统得到的测量数据通常受到现场环境等因素的影响，是带有一定信噪比的信号数据，提出利用旋转不变技术(ESPRIT)改进Matrix Pencil算法，直接以测量数据构成的数据矩阵为基础，将信号空间分解成信号子空间和噪声子空间，准确估计模型阶数，并检测出电力系统低频振荡信号不同振荡模态的振荡频率、衰减因子、振荡幅值和相位等信息，能够有效的提高计算效率和低频振荡辨识能力[8]。

中国专利201210352832.0公开了一种基于矩阵束的配电线路自适应电流速断保护方法，该发明提供了一种基于矩阵束的配电线路自适应电流速断保护方法，该发明将基于矩阵束算法的工频量提取方法引入到保护中，有效地避免了电气量中非周期分量的影响，使得工频量提取快速准确，明显改善了传统自适应电流速断保护的动作性能[9]。

参考文献：

[1]M.A.Saed.Measurement of Complex Permittivity of Low-Loss PlanarMicrowave Substrates Using Aperture-Coupled Microstrip Resonators.IEEE Trans,1993,41(8)1343-1348.

[2]唐东,张麟兮,呼斌,等.基于距离差分法消除天线测试多径干扰[J].现代电子技术,2014,37(11):101-103.

[3]Greg Hislop.Permittivity Estimation Using Coupling of CommercialGround Penetrating Radars[J].IEEE Transactions On Geoscience And RemoteSensing,vol.53,no.8,pp.4157–4164,Aug.2015.

[4]Fourestie B,Altman Z,Kanda M.Efficient detection of resonances inanechoic chambers using the matrix pencil method[J].ElectromagneticCompatibility,IEEE Transactions on,2000,42(1):1-5.

[5]FourestiéB,Altman Z,Kanda M.Anechoic chamber evaluation using thematrix pencil method[J].Electromagnetic Compatibility,IEEE Transactions on,1999,41(3):169-174.

[6]Y.Hua and T.K.Sarkar.Matrix pencil method and its performance[C].Proc.ICASSP-88,1988,4:2476-2479.

[7]Y.B.Hua and T.K.Sarkar.Matrix pencil method for estimatingparameters of exponentially damped/undamped sinusoids in noise[J].IEEE Transon Acoustics Speech and Signal Processing,1990,38(5):814-824.

[8]金涛,顾小兴,黄宴委,程远.基于Matrix Pencil的电力系统低频振荡模态辨识方法[P].中国专利:201410299068.4,2014-10-01.

[9]宋国兵,田小强,赵林平,王显峰,祁胜利,李德坤.基于矩阵束的配电线路自适应电流速断保护方法[P].中国专利:201210352832.0,2013-01-09.

发明内容

本发明的目的是为了解决现有技术中在利用反演土壤介电常数时，实际反演方法不准确，实验条件复杂的问题，提供可有效将各个界面反射分量加以区分，滤除非土壤表面反射信号及其他噪声信号，提高估计精度，便于实际运用的一种基于单基地测量的地表介质参数反演方法。

本发明包括如下步骤：

(1)在外场环境中预先将喇叭端口朝上，对准无障碍物的天空，使用矢量网络分析仪测量喇叭端口的反射信号若电缆为稳相电缆，则该反射信号也可以事先在暗室环境中测量；

(2)将喇叭垂直对准待测土壤表面，测得反射信号

(3)在待测土壤表面放置平面金属反射板，用于截获喇叭天线波束辐射的主要能量，测得反射信号

(4)将和分别减去喇叭自身反射信号得到土壤和金属板的反射信号Y₁和Y₂；

(5)分别对反射信号Y₁和Y₂进行矩阵束处理，滤除其他反射信号及噪声，得到土壤表面和金属板表面的反射信号和由于金属板反射系数已知，从而得到土壤表面的反射系数R^s＝-Y₁ ^s/Y₂ ^s；

其中η₀为空气波阻抗；另外，弱导电媒介的本征波阻抗可近似表述为

其中μ为自由空间中的磁导率，ε₀为自由空间中的介电常数，ε_r为表层土壤的相对介电常数，σ为表层土壤的电导率，ω为角频率，根据式(2)可得到表层土壤的电磁参数估算公式：

在步骤(5)中，所述分别对反射信号Y₁和Y₂进行矩阵束处理，可在基于衰减指数模型下，提取极点值，根据相位时延，估计出土壤表面和金属板表面的反射信号Y₁ ^s和Y₂ ^s；基于矩阵束方法的反射信号特征提取的算法流程如下：

步骤1：将测得的反射信号Y采样后表示为以下离散模型：

即由一组M个具有任意振幅、相位、频率和衰减因子的指数函数构成，其中R_i为复振幅，α_i为衰减因子，β_i为带有相位信息的角频率，j为虚数单位。

步骤2：利用信号样本y(n)构造一个(N-L)×(L+1)维矩阵Y,定义如下:

Y可分解为2个(N-L)×L维Hankel矩阵Y₁和Y₂,其中Y₁是由Y删除最后一列元素得到,Y₂是由Y删除第一列元素得到，即：

于是可将Y₁和Y₂表示如下：

Y₂＝Z₁RZ₀Z₂ (9)

Y₁＝Z₁RZ₂ (10)

Z₀＝diag{z₁,z₂,…,z_M} (13)

R＝diag{R₁,R₂,…,R_M} (14)

步骤3：构造矩阵束Y₂-λY₁＝Z₁R{Z₀-λI}Z₂ (15)

若λ＝z_i,i＝1,2,…,M，则参数z_i为矩阵对{Y₂:Y₁}的广义特征值，求解参数z_i等同于求解下式的本征特征值：

其中为Y₁的广义逆矩阵，其定义为：

步骤4：以上为不存在噪声等其他干扰的情况，实际运用中常常含有噪声,使真实的极点产生误差,并产生多余的虚假极点；为了减小噪声对极点提取的影响,在求Y₂-λY₁＝Z₁R{Z₀-λI}Z₂的广义特征值之前,先对Y₁和Y₂进行降秩处理。

用奇异值分解法(SVD)可将Y分解为Y＝UΣV^H,其中,U和V为酉阵,分别由YY^H和Y^HY的特征向量构成,Σ为对角阵,内为Y的奇异值(Y^HY特征值的算术平方根),并为正的非增序列。

步骤5：为得参数M，通过下式来确定，对一个奇异值σ_c,若满足

σ_c/σ_max≤10^-p (18)

则σ_c为临界奇异值(第M个奇异值),σ_max为最大奇异值,p代表数据的有效数字的位数。

步骤6：求解矩阵Y的低秩逼近矩阵，即滤除噪声后的矩阵,由矩阵V的M个主奇异向量构成的滤波矩阵:

V'＝[v₁,v₂,…,v_M] (19)

对应较小奇异值的奇异向量(v_M+1,v_M+2,…,v_L)被舍弃,因此,

其中V₁'和V₂ ^'是由V'分别删除最后一行元素和第一行元素得到,Σ'是由Σ的M个主奇异值对应的M个奇异向量组成。

步骤7：由矩阵束原理可知，信号的极点z_i就是Y₂相对于Y₁的广义特征值，求解该广义特征值的问题可以转化为求解如下矩阵的普通特征值问题：的特征值问题等效为求下式的特征值问题

在信噪比为25～30dB的情况下，以上估计方法依然适用。

步骤8：确定参数M和特征值z_i之后，幅度R可以通过求解下式的最小二乘解得到

本发明适用单喇叭天线垂直照射土壤表面，根据喇叭端口反射信号反演出地表介电常数与电导率。本发明将矩阵束方法运用到了大地电磁参数反演中，有效将各个界面反射分量加以区分，滤除非土壤表面反射信号及其他噪声信号，使得估计精度得到提高，也方便了实际运用。

本发明的优点和积极效果在于：本发明在非理想暗室环境或外场中信号测试时，能分离出不同的反射信号，并根据相位延时提取出土壤表面的反射信号，从而获得土壤表层的等效波阻抗。根据与介电常数及电导率的关系与之对应且通过数学关系转换，推导出了的反演公式，从而实现了土壤表层介电常数及电导率的反演。

在对单层和多层介质的测试验证后，证明可以进一步将本发明运用到多层介质结构中表层介质电磁参数的反演工作中，运算量小，易于工程实现。

附图说明

图1为本发明基于矢量网络分析仪测试单喇叭天线垂直照射土壤表面反射系数示意图。

图2为单个喇叭对单层介质土壤垂直照射模型图。

图3为对单层介质表面垂直照射后的反射信号再分别减去喇叭端口自身的反射信号。

图4为对金属板垂直照射后的反射信号再分别减去喇叭端口自身的反射信号。

图5为理论计算、测试获得以及经矩阵束算法处理后获得的介质层表面反射系数的比较。

图6为本发明方法反演得到的单层地表相对介电常数(真实值为3)。

图7为本发明方法反演得到的单层地表电导率(真实值为0.05S/m)。

图8为单喇叭对双层介质土壤垂直照射的模型图。

图9为本发明方法反演得到的双层介质中表层的相对介电常数(真实值为3)。

图10为本发明方法反演得到的双层介质中表层的电导率(真实值为0.05S/m)。

图11为单喇叭对3层介质土壤垂直照射的模型图。

图12为本发明方法反演得到的3层介质中表层的相对介电常数(真实值为3)。

图13为本发明方法反演得到的3层介质中表层的电导率(真实值为0.05S/m)。

具体实施方式

以下将结合附图和实施例对本发明作进一步的详细说明。

参见图1～4，本发明提供了基于单基地测量的地表介质参数反演方法，具体实现步骤如下：

(1)在外场环境中预先将喇叭1的端口朝上，对准无障碍物的天空，使用矢量网络分析仪2测量喇叭端口的反射信号若电缆为稳相电缆，则该反射信号也可以事先在暗室环境中测量；

(2)将喇叭1垂直对准待测土壤P表面，测得反射信号

(3)在待测土壤表面放置平面金属反射板(金属板具有一定的尺寸，以截获到喇叭天线波束辐射的主要能量)，测得反射信号

(4)和分别减去喇叭自身反射信号得到土壤和金属板的反射信号Y₁和Y₂

(5)分别对Y₁和Y₂进行矩阵束方法处理，在基于衰减指数模型下，提取极点值，根据相位时延，估计出土壤表面和金属板表面的反射信号和基于矩阵束方法的反射信号特征提取的算法流程如下：

步骤1：将测得的反射信号Y采样后表示为以下离散模型：

即由一组M个具有任意振幅、相位、频率和衰减因子的指数函数构成，其中R_i为复振幅，α_i为衰减因子，β_i为带有相位信息的角频率，j为虚数单位。

步骤2：利用信号样本y(n)构造一个(N-L)×(L+1)维矩阵Y,定义如下:

Y可分解为2个(N-L)×L维Hankel矩阵Y₁和Y₂,其中Y₁是由Y删除最后一列元素得到,Y₂是由Y删除第一列元素得到。即：

于是可将Y₁和Y₂表示如下：

Y₂＝Z₁RZ₀Z₂ (9)

Y₁＝Z₁RZ₂ (10)

Z₀＝diag{z₁,z₂,…,z_M} (13)

R＝diag{R₁,R₂,…,R_M} (14)

步骤3：构造矩阵束Y₂-λY₁＝Z₁R{Z₀-λI}Z₂ (15)

若λ＝z_i,i＝1,2,…,M，则参数z_i为矩阵对{Y₂:Y₁}的广义特征值。求解参数z_i等同于求解下式的本征特征值：

其中为Y₁的广义逆矩阵，其定义为：

步骤4：以上为不存在噪声等其他干扰的情况，实际运用中常常含有噪声,使真实的极点产生误差,并产生多余的虚假极点。为了减小噪声对极点提取的影响,在求Y₂-λY₁＝Z₁R{Z₀-λI}Z₂的广义特征值之前,先对Y₁和Y₂进行降秩处理。

用奇异值分解法(SVD)可将Y分解为Y＝UΣV^H,其中,U和V为酉阵,分别由YY^H和Y^HY的特征向量构成,Σ为对角阵,内为Y的奇异值(Y^HY特征值的算术平方根),并为正的非增序列。

步骤5：为得参数M，通过下式来确定，对一个奇异值σ_c,若满足

σ_c/σ_max≤10^-p (18)

则σ_c为临界奇异值(第M个奇异值),σ_max为最大奇异值,p代表数据的有效数字的位数。

步骤6：求解矩阵Y的低秩逼近矩阵，即滤除噪声后的矩阵,由矩阵V的M个主奇异向量构成的滤波矩阵:

V'＝[v₁,v₂,…,v_M] (19)

对应较小奇异值的奇异向量(v_M+1,v_M+2,…,v_L)被舍弃,因此,

其中V₁'和V₂'是由V'分别删除最后一行元素和第一行元素得到,Σ'是由Σ的M个主奇异值对应的M个奇异向量组成。

步骤7：由矩阵束原理可知，信号的极点z_i就是Y₂相对于Y₁的广义特征值，求解该广义特征值的问题可以转化为求解如下矩阵的普通特征值问题：的特征值问题等效为求下式的特征值问题

在信噪比为25～25dB的情况下，以上估计方法依然适用。

步骤8：确定了参数M和特征值z_i之后，幅度R可以通过求解下式的最小二乘解得到

在图2中，介质ε₁＝3，σ₁＝0.05s/m；喇叭口与介质上表面的距离为150mm；喇叭端口与介质上表面的距离为310mm；介质的高度为200mm。

(6)根据估计获得的反射信号，求出表层土壤的反射系数R^s＝-Y₁ ^s/Y₂ ^s，如图5为测试结果直接求得的介质表面反射系数与介质反射系数的估计值及介质反射系数的理论值的比较；

其中，η₀，η₁分别是真空中波阻抗和表层介质波阻抗。

弱导电媒介的本征波阻抗可近似为

其中μ为自由空间中的磁导率，ε₀为自由空间中的介电常数，ε_r为表层土壤的相对介电常数，σ为表层土壤的电导率，ω为角频率，根据式(18)可得到表层土壤的电磁参数估算公式：

结果如图6和7所示。

以上方法同样适用在多层介质结构的表层土壤电磁参数的估计。图8为双层模型(第1层介质ε₁＝3，σ₁＝0.05s/m；第2层介质ε₂＝10，σ₂＝0.06s/m；喇叭口与第1层介质上表面的距离为150mm；喇叭端口与第1层介质上表面的距离为310mm；第1层介质的高度为200mm，第2层介质的高度为150mm)，图9和10为反演结果，图11为3层模型(3层介质的高度均为100nm，第1层介质ε₁＝3,σ₁＝0.05s/m；第2层介质ε₂＝10，σ₂＝0.06s/m；第3层介质ε₃＝5，σ₂＝0.08s/m；喇叭口与第1层介质上表面的距离为150mm；喇叭端口与第1层介质上表面的距离为310mm)，图12和13为反演结果。

表1

表1给出三种模型土壤表层电磁参数估计结果(中心频率2.5GHz)。

Claims (2)

1.基于单基地测量的地表介质参数反演方法，其特征在于包括如下步骤：

(1)在外场环境中预先将喇叭端口朝上，对准无障碍物的天空，使用矢量网络分析仪测量喇叭端口的反射信号若电缆为稳相电缆，则该反射信号可事先在暗室环境中测量；

(2)将喇叭垂直对准待测土壤表面，测得反射信号

(3)在待测土壤表面放置平面金属反射板，用于截获喇叭天线波束辐射的主要能量，测得反射信号

(4)将和分别减去喇叭自身反射信号得到土壤和金属板的反射信号Y₁和Y₂；

(5)分别对反射信号Y₁和Y₂进行矩阵束处理，滤除其他反射信号及噪声，得到土壤表面和金属板表面的反射信号Y₁ ^s和Y₂ ^s；

(6)根据估计获得的反射信号，求出表层土壤的反射系数R^s＝-Y₁ ^s/Y₂ ^s，根据平面波对理想介质分界的垂直入射、以及反射系数与介质波阻抗的关系，获得地面表层介质的平均波阻抗

其中η₀为空气波阻抗；另外，弱导电媒介的本征波阻抗可近似表述为

其中μ为自由空间中的磁导率，ε₀为自由空间中的介电常数，ε_r为表层土壤的相对介电常数，σ为表层土壤的电导率，ω为角频率，根据式(2)可得到表层土壤的电磁参数估算公式：

相对介电常数

电导率

2.如权利要求1所述基于单基地测量的地表介质参数反演方法，其特征在于在步骤(5)中，所述分别对反射信号Y₁和Y₂进行矩阵束处理，在基于衰减指数模型下，提取极点值，根据相位时延，估计出土壤表面和金属板表面的反射信号Y₁ ^s和Y₂ ^s；基于矩阵束方法的反射信号特征提取的算法流程如下：

步骤1：将测得的反射信号Y采样后表示为以下离散模型：

即由一组M个具有任意振幅、相位、频率和衰减因子的指数函数构成，其中R_i为复振幅，α_i为衰减因子，β_i为带有相位信息的角频率，j为虚数单位；

步骤2：利用信号样本y(n)构造一个(N-L)×(L+1)维矩阵Y,定义如下:

Y可分解为2个(N-L)×L维Hankel矩阵Y₁和Y₂,其中Y₁是由Y删除最后一列元素得到,Y₂是由Y删除第一列元素得到，即：

于是可将Y₁和Y₂表示如下：

Y₂＝Z₁RZ₀Z₂ (9)

Y₁＝Z₁RZ₂ (10)

式中

Z₀＝diag{z₁,z₂,…,z_M} (13)

R＝diag{R₁,R₂,…,R_M} (14)

步骤3：构造矩阵束Y₂-λY₁＝Z₁R{Z₀-λI}Z₂ (15)

若λ＝z_i,i＝1,2,…,M，则参数z_i为矩阵对{Y₂:Y₁}的广义特征值，求解参数z_i等同于求解下式的本征特征值：

其中为Y₁的广义逆矩阵，其定义为：

步骤4：以上为不存在噪声等其他干扰的情况，实际运用中常常含有噪声,使真实的极点产生误差,并产生多余的虚假极点；为了减小噪声对极点提取的影响,在求Y₂-λY₁＝Z₁R{Z₀-λI}Z₂的广义特征值之前,先对Y₁和Y₂进行降秩处理；

用奇异值分解法(SVD)可将Y分解为Y＝U∑V^H,其中,U和V为酉阵,分别由YY^H和Y^HY的特征向量构成,∑为对角阵,内为Y的奇异值(Y^HY特征值的算术平方根),并为正的非增序列；

步骤5：为得参数M，通过下式来确定，对一个奇异值σ_c,若满足

σ_c/σ_max≤10^-p (18)

则σ_c为临界奇异值(第M个奇异值),σ_max为最大奇异值,p代表数据的有效数字的位数；

步骤6：求解矩阵Y的低秩逼近矩阵，即滤除噪声后的矩阵,由矩阵V的M个主奇异向量构成的滤波矩阵:

V'＝[v₁,v₂,…,v_M] (19)

对应较小奇异值的奇异向量(v_M+1,v_M+2,…,v_L)被舍弃,因此,

其中V′₁和V′₂是由V'分别删除最后一行元素和第一行元素得到,∑'是由∑的M个主奇异值对应的M个奇异向量组成；

步骤7：由矩阵束原理可知，信号的极点z_i就是Y₂相对于Y₁的广义特征值，求解该广义特征值的问题可以转化为求解如下矩阵的普通特征值问题：的特征值问题等效为求下式的特征值问题

在信噪比为25～30dB的情况下，以上估计方法依然适用；

步骤8：确定参数M和特征值z_i之后，幅度R通过求解下式的最小二乘解得到