CN105928987A - 基于探地雷达的盐碱地电导率测定方法 - Google Patents

Abstract

本发明公开了一种基于探地雷达的盐碱地电导率测定方法，包括如下步骤：步骤一、使用探地雷达对土壤进行探测，获取探地雷达回波相关系数ρ_ts；步骤二、使用公式计算土壤电导率σ，其中，a_i(i＝1,2,…,n)分别为i阶系数。本发明通过在盐碱地选择不同类型的典型样地进行雷达回波探测，并将雷达回波与雷达发射机的发射子波进行相关运算，得到相关系数。然后，根据样点的实测电导率与雷达回波相关系数的对应关系，得到电导率与雷达回波相关系数关系，当根据待测盐碱地的雷达回波计算出相关系数后，就可以反演出探测点的盐碱土电导率。因此，本方法可快速的测量出探测点的电导率，进而分析出盐碱程度，具有较高的准确性。

Description

基于探地雷达的盐碱地电导率测定方法

技术领域

本发明涉及雷达探测技术领域，特别涉及一种基于探地雷达的盐碱地电导率测定方法。

背景技术

目前，我国一些地区的盐碱化与沙化范围扩大速度正在加快，已引起各界人士的关注。某盐碱地改良需要得到盐碱地的盐渍化程度信息，目前一般是通过测量盐碱土电导率实现。盐碱地含盐量高，在水的溶解作用下，溶解盐离子会提高盐碱土的导电能力，因此盐碱土的电导率比一般土壤大很多。地质雷达具有解决盐碱土分布状态和分布厚度的能力，可通过地质雷达来了解盐碱土壤的含盐量、PH值、有机物含量、含水量等信息。

中国专利申请号200910089204.6公开了一种土壤盐碱化的雷达遥感监测方法，基于全极化SAR书记计算极化参数，进行裸地与植被覆盖区域的划分，生成掩膜图像，划定试验区，利用测量数据建立含水含盐土壤介电模型，利用区域复介电常数和区域含水量，进行盐土和碱土的区分。上述方法虽然能够获取土壤的含盐量，但过程复杂，准确性差。因此，需要提供一种能够快速、准确的测量出土壤含盐量的方法。

发明内容

本发明设计开发了一种基于探地雷达的盐碱地电导率测定方法，目的是通过雷达回波相关系数反演出探测点的盐碱土电导率，以得到土壤盐碱化程度。

本发明提供的技术方案为：

基于探地雷达的盐碱地电导率测定方法，包括如下步骤：

步骤一、使用探地雷达对土壤进行探测，获取探地雷达回波相关系数ρ_ts；

步骤二、计算土壤电导率σ

$\mathrm{\σ}=a_0+a_1{\mathrm{\ρ}}_{ts}+a_2{\mathrm{\ρ}}_{ts}^2+...+a_i{\mathrm{\ρ}}_{ts}^i+...+a_n{\mathrm{\ρ}}_{ts}^n,$

其中，a_i(i＝1,2,…,n)分别为i阶系数。

优选的是，步骤二中，n＝3。

优选的是，至少采集四个试验点，分别测量导电率σ(1),σ(2),…,σ(n)，并获取相应的试验点的雷达回波相关系数ρ_ts(1),ρ_ts(2),…,ρ_ts(n)，代入方程

$\mathrm{\σ}\left(i\right)=a_0+a_1{\mathrm{\ρ}}_{ts}\left(i\right)+a_2{\mathrm{\ρ}}_{ts}^2\left(i\right)+a_3{\mathrm{\ρ}}_{ts}^3\left(i\right),i=1,2,...,n$

获得a₀、a₁、a₂、a₃的解。

优选的是，采用最小二乘法计算a₀、a₁、a₂、a₃的解。

优选的是，使用土壤水分温度电导率测试仪测量试验点土壤的电导率。

本发明的有益效果是：本发明提供了一种基于探地雷达的盐碱地电导率测定方法，通过在盐碱地选择不同类型的典型样地进行雷达回波探测，并将雷达回波与雷达发射机的发射子波进行相关运算，得到相关系数。然后，根据样点的实测电导率与雷达回波相关系数的对应关系，得到电导率与雷达回波相关系数的拟合曲线，并求得以回波相关系数作为自变量的电导率计算公式。得到先验公式以后，当根据待测盐碱地的雷达回波计算出相关系数后，就可以反演出探测点的盐碱土电导率。因此，本方法可快速的测量出探测点的电导率，进而分析出盐碱程度，具有较高的准确性。

附图说明

图1为初始相位不同的同频正弦波相关系数曲线。

图2为低电导率土壤介质回波波形图。

图3为高电导率土壤介质回波波形图。

图4为反射系数曲线示意图。

图5为介电常数虚部与回波相关系数关系曲线。

图6为测线1俯视图。

图7为测线1雷达回波相关系数与电导率的关系示意图。

图8为测线2俯视图。

图9为测线2雷达回波相关系数与电导率的关系示意图。

图10为测线1,2雷达回波相关系数与电导率关系散点图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本发明提供了一种基于探地雷达的盐碱地电导率测定方法，分析雷达回波相关系数与导电率的对应关系，通过获取回波相关系数得到探测点的盐碱土电导率。

考虑一个周期时间窗的某一频率信号w₁表示为：

w₁＝A₁ cos(ω₁t)

令一同频率信号w₂表示为：

如图1所示，当取值从0～2π变化时，w₁与w₂的相关系数是相位差的余弦函数，即相关系数可以用来衡量波形之间的相位差异。

盐碱土的高损耗主要是由于土壤中溶解盐离子浓度较高所导致的电导率升高，从而使土壤的介电常数虚部增大，当介电常数虚部与实部相比不能忽略时，土壤介质的波阻抗变成复数，进而土壤界面的反射系数Г也变为复数，复反射系数会造成回波的相位偏移(这与低损耗介质回波的无相位偏移有本质的不同)，时域上体现为回波波形的畸变，如图2、图3所示。

考虑具有较高含水量与含盐量的土壤介质，由于电导率较高，介电常数虚部不能忽略，令某一频率f_n下土壤介电常数为ε_n，相对介电常数为ε_rn，电导率为σ，对由电导率表征的有损介质有介电常数表达式：

ε_n＝ε′_n-jδ/ω_n

对一般损耗介质(介电常数虚部不能由低频电导率直接得到)作线性近似：

ε″_rn＝a_1nδ_n+a_0n

其中a_1n，a_0n为常数，a_0n表征除电导率以外因素导致的介质损耗，则有：

ε_n＝ε_rn·ε₀＝(ε′_rn-jε″_rn)·ε₀＝(ε′_rn-j(a_1nδ_n+a_0n))·ε₀

其中ε₀为真空介电常数，考虑实际高电导率土壤介电常数的实部与虚部：

ε′_rn＞0，ε″_rn＞0

土壤介质在频率f_n下的波阻抗为：

$Z_n=\sqrt{u_1/{\mathrm{\ϵ}}_n}=\sqrt{\frac{{\mathrm{\μ}}_1}{{\mathrm{\ϵ}}_0}}\sqrt{\frac{1}{{\mathrm{\ϵ}}_{rn}}}=Z_0\sqrt{\frac{1}{{\mathrm{\ϵ}}_{rn}}}=Z_{nr}+{\mathrm{jZ}}_{ni}$

其中Z₀为真空波阻抗，将Z_n归一化，可写为：

$z_n=Z_n/Z_0=\sqrt{1/{\mathrm{\ϵ}}_{rn}}=z_{nr}+{\mathrm{jz}}_{ni}$

${\left(z_n\right)}^2={(z_{nr}+{\mathrm{jz}}_{ni})}^2=\left({\left(z_{nr}\right)}^2-{\left(z_{ni}\right)}^2\right)+2z_{nr}{z}_{ni}\·j=1/{\mathrm{\ϵ}}_{rn}=\frac{1}{{\left({\mathrm{\ϵ}}_{rn}^{\′}\right)}^2+{\left({\mathrm{\ϵ}}_{rn}^{\′\′}\right)}^2}({\mathrm{\ϵ}}_{rn}^{\′}+{\mathrm{j\ϵ}}_{rn}^{\′\′})$

考虑土壤介电常数实部与虚部的符号，可得：

$z_{nr}={(\frac{1}{{\left({\mathrm{\ϵ}}_{rn}^{\′}\right)}^2+{\left({\mathrm{\ϵ}}_{rn}^{\′\′}\right)}^2}\·\frac{\sqrt{{\left({\mathrm{\ϵ}}_{rn}^{\′}\right)}^2+{\left({\mathrm{\ϵ}}_{rn}^{\′\′}\right)}^2}+{\mathrm{\ϵ}}_{rn}^{\′}}{2})}^{\frac{1}{2}}>z_{ni}>0$

$z_{ni}={(\frac{1}{{\left({\mathrm{\ϵ}}_{rn}^{\′}\right)}^2+{\left({\mathrm{\ϵ}}_{rn}^{\′\′}\right)}^2}\·\frac{\sqrt{{\left({\mathrm{\ϵ}}_{rn}^{\′}\right)}^2+{\left({\mathrm{\ϵ}}_{rn}^{\′\′}\right)}^2}-{\mathrm{\ϵ}}_{rn}^{\′}}{2})}^{\frac{1}{2}}>0$

高损耗介质条件下电磁波穿透能力非常有限，假设地下无明显反射体，则可忽略土壤深层产生的反射波，而仅以土壤表层反射回波参与计算。令土壤表层在频率f_n的反射系数Г_n：

${\mathrm{\Γ}}_n=\frac{Z_n-Z_0}{Z_n+Z_0}=\frac{z_n-1}{z_n+1}=\frac{z_{nr}+{\mathrm{jz}}_{ni}-1}{z_{nr}+{\mathrm{jz}}_{ni}+1}={\mathrm{\Γ}}_{nr}+{\mathrm{j\Γ}}_{ni}$

其中，Z₀为空气波阻抗(≈377Ω)。

雷达回波(不考虑深层反射)包括直达波与地表反射波，因此频率为f_n的雷达回波幅度与相位可由一等效反射系数R_n确定：

R_n＝A_dn+Γ_n＝(A_dn+Γ_nr)+jΓ_ni

其中A_dn为直达波系数(令雷达接收的直达波与发射子波相比仅在幅度上有衰减，相位关系相同)，图4为加入直达波系数的雷达回波等效反射系数曲线，可以看到曲线上每一反射系数点沿逆时针的方向变化时，相位的变化是非线性的。

对每一反射系数的相位取余弦再取反，得到取反后的频率f_n雷达回波与发射子波的相关系数(-ρ_tsn)，该回波相关系数与介电常数虚部的关系如图5所示，介电常数虚部与回波相关系数呈单调上升的非线性关系(该曲线可用3次多项式拟合)，考虑之前的假设——介电常数虚部与电导率呈线性关系：

ε″_rn＝a_1nδ_n+a_0n (14)

故可得电导率σ与雷达回波相关系数-ρ_tsn呈与图5中曲线形式相近的单调关系。

进一步，由傅式变换可将雷达发射的脉冲信号表示为通频带内各频率分量的线性组合：

则雷达回波可表示为各频率分量反射波的线性组合：

其中各频率分量f_n的回波相关系数ρ_tsn与该频率下介电常数虚部关系曲线都与图5类似，将整体的雷达回波与雷达子波的相关系数近似表示为各频率分量相关系数的线性组合：

ρ_ts＝C₁·ρ_ts1+C₂·ρ_ts2+…+C_n·ρ_tsn

则雷达回波与雷达发射子波的总相关系数ρ_tsn与电导率σ也呈近似图5的非线性关系。

令雷达发射子波为w_t(m)，高电导率土壤介质回波为w_sk(m)，取相关系数表征回波的相位偏移：

${\mathrm{\ρ}}_{ts}\left(k\right)=\frac{\underset{m=1}{\overset{M}{\Σ}}\left(w_t\right(m)-\stackrel{\‾}{w_t})\left(w_{sk}\right(m)-\stackrel{\‾}{w_s})}{\sqrt{\underset{m=1}{\overset{M}{\Σ}}{\left(w_t\right(m)-\stackrel{\‾}{w_t})}^2\·\underset{m=1}{\overset{M}{\Σ}}{\left(w_{sk}\right(m)-\stackrel{\‾}{w_{sk}})}^2}}$

经盐碱地实验发现相关系数与介质电导率有较好的对应关系，令不同电导率土壤的相关系数集合为{-ρ_ts(1),...,-ρ_ts(k),...}(相关系数取反来获得与电导率一致的对应关系)，对应的实测土壤电导率集合为{σ(1),...,σ(k),...}，令

$\left\{\begin{array}{c}\mathrm{\δ}\left(1\right)=a_0+a_1(-{\mathrm{\ρ}}_{ts}(1\left)\right)+a_2{(-{\mathrm{\ρ}}_{ts}(1\left)\right)}^2+a_3{(-{\mathrm{\ρ}}_{ts}(1\left)\right)}^3\\ .\\ .\\ .\\ \mathrm{\δ}\left(k\right)=a_0+a_1(-{\mathrm{\ρ}}_{ts}(k\left)\right)+a_2{(-{\mathrm{\ρ}}_{ts}(k\left)\right)}^2+a_3{(-{\mathrm{\ρ}}_{ts}(k\left)\right)}^3\\ .\\ .\\ .\end{array}\right.$

解超定方程得到系数a₀～a₃的估计值，则可以根据雷达回波的相关系数ρ_ts(k)反演土壤电导率σ(k)。

首先使用探地雷达对土壤进行探测，获取至少四个试验点处的探地雷达回波相关系数ρ_ts(1),ρ_ts(2),…,ρ_ts(m)，并使用土壤水分温度电导率测试仪测量对应的试验点处的电导率σ(1),σ(2),…,σ(m)。其中m为大于或等于四的整数。

然后将相关系数ρ_ts(1),ρ_ts(2),…,ρ_ts(m)和电导率σ(1),σ(2),…,σ(m)带入到下式中

$\mathrm{\σ}\left(i\right)=a_0+a_1{\mathrm{\ρ}}_{ts}\left(i\right)+a_2{\mathrm{\ρ}}_{ts}^2\left(i\right)+a_3{\mathrm{\ρ}}_{ts}^3\left(i\right),i=1,2,...,m$

使用最小二乘法计算a₀、a₁、a₂、a₃的解。

最终获得了相关系数与电导率的对应关系：

$\mathrm{\σ}=a_0+a_1{\mathrm{\ρ}}_{ts}+a_2{\mathrm{\ρ}}_{ts}^2+a_3{\mathrm{\ρ}}_{ts}^3.$

利用上式，只需使用探地雷达对土壤进行探测，获取其中一个探测点的相关系数后即可得到该探测点的电导率。通过电导率即可分析出该处土壤的盐碱程度。

例如在每一次在探地雷达实验中，选取具有代表性的两条测线说明雷达回波与样地水盐特性之间的关系。

如图6所示，测线1样地，测线长度90m，该样地地表特征多样，测线最左侧是一片地表开裂的干涸小水塘，由于地表开裂严重，初步估计在雨水淋溶与冲刷作用下，盐分很大一部分被水从裂缝处带到深层的地下，因此浅表土层的含盐量应该较低；水塘右侧经过4米左右的裸土地带后，是一片宽度2米左右的有稀疏植被的条带；条带右侧又是长度为20多米的裸土地带；再右侧是宽度大约3米左右的具有细密植被的草地条带；草地右侧经过另外一片裸土后，是最右侧的一片干涸的水塘，该水塘地表特征与最左侧的水塘不同，地表无裂纹，土壤疏松，与重度盐碱地的致密土质有很大不同，疏松的土壤在雨水淋溶渗透作用下，盐分上下输送通畅，因此初步估计含盐量也不高。

根据雷达回波成像，可以清晰分辨干涸的水塘、裸土、稀疏植被、草地所在的区域界限。在测线雷达经过的每一点，使用WET土壤水分温度电导率测试仪测量了土壤体积含水量、介电常数实部与电导率。

图7为电导率曲线与雷达回波相关系数曲线对比图，图中可以看到两条曲线的变化趋势高度契合，体现了雷达回波相位变化与盐碱土电导率的相关关系。由于WET单点测量范围大约数平方厘米，探地雷达波束覆盖范围大约1平方米，2种测量手段的尺度不同，因此造成了在地表介电特征急剧变化的位置，2者的趋势虽然一致性好，但测量值绝对差异较大，契合度相对较低。

如图8所示，测线2样地，测线长度30m，该样地最右侧紧邻一个较大的水塘，地势从草地右侧开始呈明显的左高右低(降水时，由于盐碱土透水性差，右侧地表的很大一部分积水会流向右边的水塘，同时冲刷沿途的地面，造成盐度降低)。基于测线1的先验知识，根据雷达回波成像结果，可初步分析测线15m以前(左侧)电导率较高，15米以后(右侧靠近水塘的部分)电导率急剧降低，而含水量变化不大(虽然目测地表特征仍为典型的重度盐碱地)。

图9为WET电导率曲线与雷达回波相关系数曲线对比图，与测线1的结果类似，WET电导率曲线与雷达回波相关系数曲线趋势的契合度仍然较好，电导率的变化趋势结果与雷达回波成像结果也大致吻合。

图10为测线1,2所有WET电导率与雷达回波相关系数的对应关系散点图，由研究方案部分得到的电导率与雷达回波相关系数的对应关系，这里采用3次多项式进行曲线拟合：

y＝9.084·x³+20.02·x²+19.84·x+8.81

如图10中曲线所示。利用该公式可在一定含水量条件下，根据雷达回波快速确定当前盐碱土的电导率。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (5)

1.基于探地雷达的盐碱地电导率测定方法，其特征在于，包括如下步骤：

步骤一、使用探地雷达对土壤进行探测，获取探地雷达回波相关系数ρ_ts；

步骤二、计算土壤电导率σ

$\mathrm{\σ}=a_0+a_1{\mathrm{\ρ}}_{ts}+a_2{\mathrm{\ρ}}_{ts}^2+...+a_i{\mathrm{\ρ}}_{ts}^i+...+a_n{\mathrm{\ρ}}_{ts}^n,$

其中，a_i(i＝1,2,…,n)分别为i阶系数，n为阶数。

2.根据权利要求1所述的基于探地雷达的盐碱地电导率测定方法，其特征在于，步骤二中，n＝3。

3.根据权利要求2所述的基于探地雷达的盐碱地电导率测定方法，其特征在于，至少采集四个试验点，分别测量导电率σ(1),σ(2),…,σ(m)，并获取相应的试验点的雷达回波相关系数ρ_ts(1),ρ_ts(2),…,ρ_ts(m)，代入方程

$\mathrm{\σ}\left(i\right)=a_0+a_1{\mathrm{\ρ}}_{ts}\left(i\right)+a_2{\mathrm{\ρ}}_{ts}^2\left(i\right)+a_3{\mathrm{\ρ}}_{ts}^3\left(i\right),i=1,2,...,m$

获得a₀、a₁、a₂、a₃的解。

4.根据权利要求3所述的基于探地雷达的盐碱地电导率测定方法，其特征在于，采用最小二乘法计算a₀、a₁、a₂、a₃的解。

5.根据权利要求4所述的基于探地雷达的盐碱地电导率测定方法，其特征在于，使用土壤水分温度电导率测试仪测量试验点土壤的电导率。