一种土壤盐碱化的雷达遥感监测方法

摘要

本发明公开了一种土壤盐碱化的雷达遥感监测方法，基于全极化SAR数据计算极化参数，进行裸地与植被覆盖区的划分，生成掩膜图像；划定试验区，利用测量数据建立含水含盐土壤介电模型；利用全极化SAR数据反演复介电常数图并进行RD地理编码、数据信息分离与转化；基于介电模型及复介电常数，利用遗传算法进行反演，得到区域含水量和含盐量，综合地理坐标信息和植被覆盖掩膜图，生成含水量图和含盐量图；利用区域复介电常数和含水量，进行盐土和碱土的区分。采用了本发明的技术方案，能够获取盐渍区土壤的含水量和含盐量，而且对地面辅助数据的依赖程度减小，同时可以很好地刻画含水含盐土壤的介电行为。

说明书

技术领域

本发明涉及微波遥感应用技术领域，尤其涉及一种土壤盐碱化的雷达遥感监测方法。

背景技术

微波遥感技术是一门较为年轻的学科，我国在该领域的研究起步于上个世纪70年代。研究工作从过去的单波段、单极化的机载合成孔径雷达(Synthetic Aperture Radar，SAR)数据定性分析已经发展到了现阶段的多波段、多极化、多模式的星载SAR数据定量分析。随着高分辨率的全极化SAR卫星的发射，无论从理论上，还是从技术手段上，微波遥感技术都将进入一个全新的发展时期，这需要开展大量的基础性研究。目前土壤盐碱化是世界上面临的一个严重问题，加强土壤盐碱化的监测具有极其重要的价值。传统的野外定位观测法不仅费时、费力，而且观测点少，代表性差，无法实现大面积动态监测。光学遥感不论是目视判读还是计算机自动分类，都要依赖于盐碱土的光谱响应特征，而单纯通过光谱响应来提取土壤盐碱化信息，具有一定局限性。而且，光学数据对于大气情况的依赖较大，高质量的数据获取困难。研究表明，微波遥感是电磁波谱中唯一能够真实定量化估计土壤湿度的电磁波谱频段。另外，在微波波段土壤盐分变化会影响土壤的导电性，进而影响土壤复介电常数，以及土壤表面后向散射系数。这就为通过雷达遥感监测土壤含盐量的变化提供了可能。

全极化SAR技术的理论研究由来已久，以往由于数据源稀少，使得该项技术的开展受到了限制。随着最近国内外一些全极化SAR卫星的成功发射，数据源方面得到了保证，该项技术受到了越来越多的关注。全极化技术可以从电磁波传播机理上对地物表面的粗糙程度和地物的电学性质进行分别表征，并可以计算出代表地物某方面性质的指标量，实现地物区分和定量分析。

由于地物的介电性质和几何形态是SAR数据的两个主要影响因素，因此开展地物的介电性质研究对于利用微波遥感技术反演地物参数是十分必要的。目前，土壤的介电性质研究大多针对于含水土壤。Dobson等人提出了一个含水土壤的介电模型，用于分析含水土壤的介电性质，进而作为反演土壤水分的理论基础。

第一步：测量含水土壤的容重(ρ_b)和粒度分布(粘土C、砂土S)。

第二步：根据经验公式确定模型中间辅助参量：形状参量(α)、粒度相关参量(β)和有效电导率(σ_eff)：

α＝0.65

β′＝1.2748-0.519S-0.152C

β″＝1.33797-0.603S-0.166C

σ_eff＝-1.645+1.939ρ_b-2.25622S+1.594C

根据理论公式计算自由水(fw)复介电常数：

${ϵ}_{\mathrm{fw}}^{\prime }={ϵ}_{w\infty }+\frac{{ϵ}_{w0}-{ϵ}_{w\infty }}{1+{\left(2\mathrm{\pi f}{\tau }_w\right)}^2}$

${ϵ}_{\mathrm{fw}}^{\prime \prime }=\frac{2\mathrm{\pi f}{\tau }_w({ϵ}_{w0}-{ϵ}_{w\infty })}{1+{\left(2\mathrm{\pi f}{\tau }_w\right)}^2}+\frac{{\sigma }_{\mathrm{eff}}}{2\pi {ϵ}_{0}f}\frac{({\rho }_s-{\rho }_b)}{{\rho }_s{m}_v}$

其中，ε₀是自由空间电容率，τ_w是水的弛豫时间，ε_w0是水的静态介电常数，ε_w∞是自由水(fw)的介电常数高频限值，ρ_s是土粒密度，以上各参量均为经验常量。f是频率(赫兹)，m_v是体积含水量。

第三步：建立含水土壤介电模型：

${ϵ}^{\prime }={[1+\frac{{\rho }_b}{{\rho }_s}({ϵ}_s^{\alpha }-1)+m_v^{{\beta }^{\prime }}{ϵ}_{\mathrm{fw}}^{\prime \alpha }-m_v]}^{1/\alpha }$

${ϵ}^{\prime \prime }={\left[m_v^{{\beta }^{\prime \prime }}{ϵ}_{\mathrm{fw}}^{\prime \prime \alpha }\right]}^{1/\alpha }$

其中，ε′和ε″分别为含水土壤的复介电常数的实部和虚部，ε_s为土粒介电常数，是经验常量。

第四步：以SAR数据反演出的区域复介电常数图作为输入，基于建立的模型进行单参数非线性反演，得到区域土壤含水量图。

该技术方法存在以下四点不足：

1、需要土壤粒度数据作为先验基础数据，阻碍大范围反演土壤含水量的实施；

2、第二步中需要计算的中间辅助参量都是基于经验公式的，存在局域适用性，不易于大范围推广；

3、对土壤复介电常数虚部的影响因素考虑得不充分，使得该技术只适用于含水土壤，对于广泛存在的盐渍化土壤，该技术无法表达其介电特征；

4、传统的SAR技术反演复介电常数存在一定的局限性，通常需要多角度、多极化数据才能进行较好的反演，这使得数据获取难以保证，监测周期随之加长。

发明内容

本发明的目的在于提出一种土壤盐碱化的雷达遥感监测方法，能够获取盐渍区土壤的含水量和含盐量，而且对地面辅助数据的依赖程度减小，同时可以很好地刻画含水含盐土壤的介电行为。

为达此目的，本发明采用以下技术方案：

一种土壤盐碱化的雷达遥感监测方法，包括以下步骤：

A、基于全极化合成孔径雷达数据进行同极化相关分析，获得相关系数a，相位差phi和极化度p；

B、根据相关系数a，相位差phi和极化度p划分土壤是裸地还是植被覆盖区域，对裸地生成掩膜图像；

C、根据裸地掩膜图划定试验区，在试验区中均匀采集不少于10个试验点土壤样本，测定所述土壤的含水量、含盐量和复介电常数；

D、利用试验点土壤样本数据建立适用于试验区含水含盐土壤介电模型；

E、针对试验区，基于全极化合成孔径雷达数据反演出区域复介电常数，并基于距离-多普勒(Range-Doppler，RD)算法进行地理编码，生成图件；

F、将复介电常数和地理坐标信息分离，分别存储，复介电常数转化为文本形式，并保存映射关系；

G、根据试验区含水含盐土壤介电模型和复介电常数，利用遗传算法(Genetic Algorithm，GA)进行多参数非线性反演，获得试验区区域含水量和含盐量，并利用映射关系进行反变换，形成区域含水量图和含盐量图；

H、利用地理坐标信息对区域含水量图和含盐量图进行地理编码；

I、利用区域复介电常数和含水量，进行盐土和碱土的区分。

步骤B中，如果相关系数a大于0.7、相位差phi接近于0并且极化度p大于0.7，划分为裸地，相反划分为植被覆盖区域。当然这些取值均为经验值，在实际操作过程中需要根据情况调整。

还包括以下步骤：

生成植被覆盖区域掩膜图；

基于植被覆盖区域掩膜图，对步骤G中的区域含水量图和含盐量图进行无效数据的标注。

步骤D进一步包括以下步骤：

利用试验点土壤样本数据训练含水含盐土壤介电模型实部，公式为：

${ϵ}^{\prime }={[1+\frac{{\rho }_b}{{\rho }_s}({ϵ}_s^{\alpha }-1)+m_v^{{\beta }^{\prime }}{ϵ}_{\mathrm{sw}}^{{\prime \alpha }^{\prime }}-m_v]}^{1/\alpha };$

利用试验点土壤样本数据训练含水含盐土壤介电模型虚部，其中，低频公式为：

${ϵ}_{\mathrm{low}}^{\prime \prime }=m_v^{{\beta }^{\prime \prime }/\alpha }[B^{{\alpha }^{\prime \prime }/\alpha }+B^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha B}}A+O\left(A^2\right)]=m_v^{{\beta }^{\prime \prime }/\alpha }{B}^{{\alpha }^{\prime \prime }/\alpha }+m_v^{{\beta }^{\prime \prime }/\alpha }{B}^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha B}}A,$

高频公式为：

${ϵ}_{\mathrm{high}}^{\prime \prime }=m_v^{{\beta }^{\prime \prime }/\alpha }[A^{{\alpha }^{\prime \prime }/\alpha }+A^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha A}}B+O\left(B^2\right)]=m_v^{{\beta }^{\prime \prime }/\alpha }{A}^{{\alpha }^{\prime \prime }/\alpha }+m_v^{{\beta }^{\prime \prime }/\alpha }{A}^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha A}}B;$

其中，α′，α″，β′和β″为具有实际物理意义的拟合参量，$A=\frac{2\mathrm{\pi f}{\tau }_{\mathrm{sw}}(T,0)b(N,T)[{ϵ}_{\mathrm{sw}0}(T,0)a\left(N\right)-{ϵ}_{\mathrm{sw}\infty }]}{1+{\left[2\mathrm{\pi f}{\tau }_{\mathrm{sw}}(T,0)b(N,T)\right]}^2},$$B=\frac{{\sigma }_i}{2\pi {ϵ}_{0}f},$τ_sw(T，0)是盐水的弛豫时间，ε_sw0(T，0)是盐水的静态介电常数，ε_sw∞是盐水的介电常数高频限值，a(N)和b(T，N)是与含盐量有关的参量，σ_i是与含水量和含盐量有关的电导率。

还包括以下步骤：

根据试验区含水含盐土壤介电模型采用以下公式进行试验区域土壤粒度分布的反演：

$(m_v-V_{\mathrm{bw}}){ϵ}_{\mathrm{fw}}^{{\alpha }^{\prime }}+V_{\mathrm{bw}}{ϵ}_{\mathrm{bw}}^{{\alpha }^{\prime }}=m_v^{{\beta }^{\prime }}{ϵ}_{\mathrm{fw}}^{{\alpha }^{\prime }}\Rightarrow V_{\mathrm{bw}}=(m_v-m_v^{{\beta }^{\prime }})/[1-{\left(\frac{{ϵ}_{\mathrm{bw}}}{{ϵ}_{\mathrm{fw}}}\right)}^{{\alpha }^{\prime }}]\Rightarrow \mathrm{Clay}$

$\left(\begin{array}{c}\mathrm{Clay}={\int }_{-\infty }^2\mathrm{aexp}(-0.5\frac{x-x_0}{b})\mathrm{dx}\\ 1={\int }_{-\infty }^{+\infty }\mathrm{aexp}(-0.5\frac{x-x_0}{b})\mathrm{dx}\\ a=f\left(b\right)\end{array}\right)$

其中，V_bw是束缚水体积，与粘土(clay)含量具有强相关关系，ε_bw和ε_fw分别为束缚水和自由水的复介电常数实部，为35和80。

步骤F中，复介电常数的实部和虚部分别存储，结构一致。

步骤I进一步包括步骤：

将含水量图和复介电常数虚部图进行对比分析，提取出含水量差异较大的不少于三组(ε”，M_v)；

计算ε”/M_v，以M_v由小到大顺序，排序各个比值；

比值逐渐增大的区域为盐土，比值逐渐减小的该区域为碱土。

步骤C中，复介电常数测量频率与合成孔径雷达系统工作频率一致。

采用了本发明的技术方案，提出了一套针对于盐渍区土壤的可行性含水量和含盐量反演技术方案，模型对地面辅助数据的依赖程度减小，并可以很好地刻画含水含盐土壤的介电行为。利用实验室自行配置的含水含盐土壤样本的测试数据，对模型精度进行验证，实部精度可以达到89.5％，高频和低频条件下虚部精度分别为93.6％和91.7％。由于地物的介电性质是影响雷达后向散射能力的主要影响因素之一，因此本发明为利用雷达遥感技术开展土壤盐碱化监测提供了技术基础和理论依据，而且本技术采用全极化SAR数据进行分析，不涉及多源数据的空间尺度和波谱尺度的统一问题，数据获取周期短，不受天气状况影响，便于实施快速监测。另外，本技术方案对先验数据依赖性小，中间参量不采取经验公式，使其适用于大范围实施，在空间尺度上与遥感技术吻合。同时，本发明还提出了反演土壤粒度分布信息和区分土壤盐碱化性质的方法流程。

附图说明

图1是本发明具体实施方式中土壤盐碱化的雷达遥感监测流程图。

具体实施方式

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。

本发明技术方案的主要思想是首先基于全极化SAR数据进行同极化相关分析和极化度的计算，并根据极化参数进行裸地(包括低植被覆盖区)与植被覆盖区的划分，生成掩膜图像；其次根据裸地(包括低植被覆盖区)掩膜图划定试验区，利用若干样点测量数据进行含水含盐土壤介电模型的建立；再利用SAR全极化数据反演试验区复介电常数图及相应的距离-多普勒(Range-Doppler，RD)地理编码、数据信息分离与转化等操作；然后基于建立好的介电模型及复介电常数输入数据，利用遗传算法(Genetic Algorithm，GA)进行多参数非线性反演，得到区域含水量和含盐量，综合地理坐标信息和植被覆盖掩膜图，生成最终的含水量效果图和含盐量效果图；最后利用区域复介电常数数据和含水量数据，进行盐土和碱土的区分。

图1是本发明具体实施方式中土壤盐碱化的雷达遥感监测流程图。如图1所示，该流程包括以下步骤：

步骤101、由于植被覆盖对监测土壤盐碱化是一种干扰，需要利用极化参数分析地物散射机制，划分区域植被覆盖度，为后续技术的开展提供先验知识，因此首先基于全极化合成孔径雷达数据进行同极化相关分析，获得相关系数a，相位差phi和极化度p。

步骤102、这些极化参数可以表征地物的散射机制，因此可以根据其值的分布进行裸地(包括低植被覆盖区)与植被覆盖区的划分，生成掩膜图像。如果相关系数a大于0.7、相位差phi接近于0并且极化度p大于0.7，划分为裸地，相反划分为植被覆盖区域。当然这些取值均为经验值，在实际操作过程中需要根据情况调整。裸地(包括低植被覆盖区)掩膜图为划分试验区提供参考，植被覆盖区掩膜图为最后生成含水量和含盐量效果图提供无效数据标注依据。

步骤103、根据裸地(包括低植被覆盖区)掩膜图划定试验区，在试验区中均匀采集若干试验点土壤样本(采集的试验点数量越多，建立的模型越准确，一般不少于10个)，测定其含水量、含盐量和复介电常数。其中复介电常数测量频率要与SAR系统工作频率一致。

步骤104、利用实测的试验点土壤样本数据对含水含盐土壤介电模型的拟合参数进行特定范围内的全局最优拟合，进而建立适用于试验区含水含盐土壤介电模型，根据SAR系统工作频率，选择利用虚部模型中的高频或低频公式进行训练：

利用试验点土壤样本数据训练含水含盐土壤介电模型实部，公式为：

${ϵ}^{\prime }={[1+\frac{{\rho }_b}{{\rho }_s}({ϵ}_s^{\alpha }-1)+m_v^{{\beta }^{\prime }}{ϵ}_{\mathrm{sw}}^{{\prime \alpha }^{\prime }}-m_v]}^{1/\alpha };$

利用试验点土壤样本数据训练含水含盐土壤介电模型虚部，其中，低频公式为：

${ϵ}_{\mathrm{low}}^{\prime \prime }=m_v^{{\beta }^{\prime \prime }/\alpha }[B^{{\alpha }^{\prime \prime }/\alpha }+B^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha B}}A+O\left(A^2\right)]=m_v^{{\beta }^{\prime \prime }/\alpha }{B}^{{\alpha }^{\prime \prime }/\alpha }+m_v^{{\beta }^{\prime \prime }/\alpha }{B}^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha B}}A,$

高频公式为：

${ϵ}_{\mathrm{high}}^{\prime \prime }=m_v^{{\beta }^{\prime \prime }/\alpha }[A^{{\alpha }^{\prime \prime }/\alpha }+A^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha A}}B+O\left(B^2\right)]=m_v^{{\beta }^{\prime \prime }/\alpha }{A}^{{\alpha }^{\prime \prime }/\alpha }+m_v^{{\beta }^{\prime \prime }/\alpha }{A}^{{\alpha }^{\prime \prime }/\alpha }\frac{{\alpha }^{\prime \prime }}{\mathrm{\alpha A}}B;$

其中，α′，α″，β′和β″为具有实际物理意义的拟合参量，

$A=\frac{2\mathrm{\pi f}{\tau }_{\mathrm{sw}}(T,0)b(N,T)[{ϵ}_{\mathrm{sw}0}(T,0)a\left(N\right)-{ϵ}_{\mathrm{sw}\infty }]}{1+{\left[2\mathrm{\pi f}{\tau }_{\mathrm{sw}}(T,0)b(N,T)\right]}^2},$$B=\frac{{\sigma }_i}{2\pi {ϵ}_{0}f},$τ_sw(T，0)是盐水的弛豫时间，ε_sw0(T，0)是盐水的静态介电常数，ε_sw∞是盐水的介电常数高频限值，a(N)和b(T，N)是与含盐量有关的参量，σ_i是与含水量和含盐量有关的电导率。

步骤105、根据试验区含水含盐土壤介电模型采用以下公式进行试验区域土壤粒度分布的反演：

$(m_v-V_{\mathrm{bw}}){ϵ}_{\mathrm{fw}}^{{\alpha }^{\prime }}+V_{\mathrm{bw}}{ϵ}_{\mathrm{bw}}^{{\alpha }^{\prime }}=m_v^{{\beta }^{\prime }}{ϵ}_{\mathrm{fw}}^{{\alpha }^{\prime }}\Rightarrow V_{\mathrm{bw}}=(m_v-m_v^{{\beta }^{\prime }})/[1-{\left(\frac{{ϵ}_{\mathrm{bw}}}{{ϵ}_{\mathrm{fw}}}\right)}^{{\alpha }^{\prime }}]\Rightarrow \mathrm{Clay}$

$\left(\begin{array}{c}\mathrm{Clay}={\int }_{-\infty }^2\mathrm{aexp}(-0.5\frac{x-x_0}{b})\mathrm{dx}\\ 1={\int }_{-\infty }^{+\infty }\mathrm{aexp}(-0.5\frac{x-x_0}{b})\mathrm{dx}\\ a=f\left(b\right)\end{array}\right)$

其中，V_bw是束缚水体积，与粘土(clay)含量具有强相关关系，ε_bw和ε_fw分别为束缚水和自由水的复介电常数实部，为35和80。一般来说，自然界中的土壤粒度分布都具有一定的区域一致性，这使得高斯分布参数a和b在试验区域内存在某种函数关系。另外，考虑到粘土含量定义为粒径小于2微米的土粒比例，因此综合以上分析可以建立上述方程组，解出相应参数a、b和x₀即可确定与实际土壤粒度分布相似的高斯分布，进而得到粒度分布方面的信息。

步骤106、针对试验区，基于全极化合成孔径雷达数据反演出区域复介电常数，并基于距离-多普勒算法进行地理编码，生成图件。

步骤107、将复介电常数和地理坐标信息分离，分别存储，复介电常数转化为文本形式，复介电常数的实部和虚部分别存储，结构一致，并保存映射关系，以便将反演结果重新转化为图像形式。

步骤108、根据试验区含水含盐土壤介电模型和复介电常数，利用遗传算法进行多参数非线性反演，辅助通用全局优化法，使结果达到全局最优，获得试验区区域含水量和含盐量，并利用映射关系进行反变换，形成区域含水量图和含盐量图。

步骤109、基于植被覆盖区域掩膜图，对区域含水量图和含盐量图进行无效数据的标注。

步骤110、利用地理坐标信息对区域含水量图和含盐量图进行地理编码，最终生成区域含水量效果图和含盐量效果图。

步骤111、利用区域复介电常数和含水量，进行盐土和碱土的区分。

根据实验室盐碱化土壤样本测量数据分析，发现盐土和碱土的ε”随M_v的变化规律均符合抛物线分布，但盐土的分布曲线斜率逐渐增大，而碱土正好相反。

因此首先将含水量图和复介电常数虚部图进行对比分析，提取出含水量差异较大的不少于三组(ε”，M_v)；然后计算ε”/M_v，以M_v由小到大顺序，排序各个比值；最后判断比值逐渐增大的区域为盐土，比值逐渐减小的该区域为碱土。

采用了上述技术方案，不需要土壤的粒度数据作为先验基础数据，并可以通过模型中具有实际物理意义的参数进行试验区域土壤粒度分布信息的反演；中间辅助参量不采用经验公式进行表达，而是令各参量在特定范围内自动拟合，使得各参量具有的实际的物理意义，并可以由此开展扩展性应用(如盐土和碱土的区分)；基于含水含盐土壤介电性质的详细分析，对土壤复介电常数虚部的影响因素进行了充分考虑，并对其在高频和低频下的不同规律进行了定量化的描述，使其适用于含水含盐土壤，进而建立土壤含水量和含盐量反演的技术流程；采用SAR全极化分析技术，数据获取周期短，并可以对地表粗糙度和地表介电性质进行较好的分离。另外，极化参数可以先期辅助进行地物的简单区分，分析结果作为后续反演流程的先验知识，进一步提高参数反演精度。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。