CN103808736B - 基于被动微波混合像元分解技术的盐碱地特性探测方法 - Google Patents

Abstract

本发明的基于被动微波混合像元分解技术的盐碱地特性探测方法，属于遥感图像处理的技术领域。通过光谱遥感数据获得观测地区地物分类数据，建立针对盐碱地类型的被动微波混合像元分解模型；采用带约束条件的非负最小二乘法迭代运算求解方程组，实现盐碱地被动微波混合像元分解；进一步建立分解后的被动微波遥感数据与盐碱地含盐量的关系模型，实现盐碱地特性的有效探测。本发明为深入研究苏打盐碱地水盐特性与波谱信息的相关性，及苏打盐碱土的导电特性与微波辐射/散射信息的相关性，掌握苏打盐碱地的基本形态、时空动态变化过程、演变规律和地理分布特性，及苏打盐碱地的改良与合理利用，提供理论方法和技术手段。<pb pnum="1" />

Description

基于被动微波混合像元分解技术的盐碱地特性探测方法

技术领域

本发明属于遥感图像处理的技术领域，利用光学和微波遥感在盐碱地识别上的优势，提出一种基于被动微波混合像元分解技术的盐碱地特性探测方法，为掌握苏打盐碱地的基本形态、时空动态变化过程、演变规律和地理分布特性，及盐碱地的改良与合理利用，提供理论方法和技术手段。

背景技术

土壤盐碱化是干旱、半干旱区主要的土地退化问题，是当今土地荒漠化的主要问题之一。土地的盐碱化作为一种环境灾害，导致土地的退化，从而削弱和破坏了土地的生产力，使农业区粮食产量下降，严重威胁着生态及国民经济的可持续发展。传统的野外定位观测费时费力，而且观测点少，无法进行大面积动态监测。遥感技术能大面积地获取观测地区的有用信息，为监测地面的信息变化提供了可靠的依据，被广泛地应用于盐碱地的动态监测和资源调查中。

由于澳大利亚、美国、埃及等许多国家存在着不同程度的土地盐碱化，而且面积都很大，所以土地盐碱化研究很早就成为国外的热门话题。主要集中在光学遥感盐碱地特性和微波遥感盐碱地特性研究。

（1）光学遥感盐碱地特性的研究进展

20世纪70年代，国外利用卫星遥感技术对土地盐碱化进行监测研究。Rao和Dwivedi等在利用遥感技术研究盐碱土时，使用OIF(Optimum Index Factor)技术，分析了TM数据的波段组合模式与含盐碱化信息量的关系，并指出提取盐碱化土壤信息的精度与遥感数据信息量的大小不成正比关系（参见Rao B RM,Dwivedi R S,etal.Mapping the Magnitude of Sodicity in Part of the Indo-GangeticPlain of Uttar Pradesh,Northern India Using Landsat-TM Data[J].InternationalJournal of Remote Sensing,1991,12(3):419-425）。Rao对盐碱化土壤光谱特征分析发现，与未盐碱化土相比，盐碱化土壤在可见光和近红外波段有很强的反射率，并且土壤盐碱化程度越高，反射光谱越强（参见Rao B R M，Sankar T R DwivediR,et al.Spectral Behaviour of Salt-affected Soils.International Journal of RemoteSensing1995,16(12):2125-2136）。2003年，Bui E.N通过对澳大利亚东北部的植被种类，以及植被的分布与集群进行了深入地研究，确定了土地盐碱化的程度（参见Bui E N,Henderson B L.Vegetation indicators of salinity in northern Queensland[J].Austral Ecology,2003,28:539–552）。国内开展土壤盐渍化卫星遥感监测研究始于上世纪80年代，比国外大约晚10年，但是盐渍化方面的研究不断地得到新的成果。我国学者对土壤盐渍化的研究工作集中在盐渍化土壤的植被及景观生态、土壤盐分积累及运移模型、生态地质环境及水化学环境等方面。利用遥感影像进行目视判读是进行盐渍土定性、定量和动态分析的重要手段，数字图像处理技术在早期的盐渍土监测研究方面也发挥了一定的作用（参见扶卿华等.土壤盐分含量的遥感反演研究.农业工程学报,2007,(01)；吴景坤,章兆兴,王爱军.库尔勒盐渍土的遥感图像处理.遥感信息,1987,(1):26）。

（2）微波遥感盐碱地特性的研究进展

从色谱学角度而言，在实验室内利用高光谱信息能够有效地识别盐碱土的含盐类型。在自然条件下，由于水和其它环境条件的作用降低了含盐光谱信息的灵敏度，导致利用高光谱信息难以实现盐碱土特性的定量反演。在微波波段，含水含盐土壤的变化会影响土壤的导电性，进而影响土壤介电常数，而介电常数的变化直接影响其后向散射系数或微波辐射亮温值（参见Thomas J.Jackson P E.NellO.Salinity Effects on the Microwave Emission of Siols[J].IEEE Transaction onGeoscience and Remote Sensing.1987,28(2):214-220；邵芸,吕远,董庆,等.含水含盐土壤的微波介电特性分析研究[J],遥感学报,2002,6(6):416-423），这为通过微波遥感监测土壤含盐量提供了可能。利用微波遥感手段对盐碱化土壤含盐量的研究比较少，而且大多是初步的定性的探讨。Jackson and Qneiu(参见Thomas J.JacksonP E.Nell O.Salinity Effects on the Microwave Emission of Siols[J].IEEETransaction on Geoscience and Remote Sensing.1987,28(2):214-220.)从他们的野外实验中发现在一定的土壤含水量之下，盐度的增加会造成辐射的降低。

综上所述，国内外学者在土壤盐分的遥感监测上己做了很多研究，取得了许多研究成果。利用微波遥感手段对盐碱化土壤含盐量的研究大多是初步的探讨，因此，在这一领域的研究主要还是在单一光谱遥感数据源层面上，没有将光学与微波遥感相结合，共同实现观测地区盐碱特性探测。在遥感探测盐碱地特性方面所取得的技术成果还无法满足盐碱地改良的应用需求。

发明内容

本发明要解决的技术问题是，提供一种基于被动微波混合像元分解技术的盐碱地特性探测方法，将光学和被动微波遥感技术相结合，共同实现观测地区盐碱地特性的有效探测。

通过光谱遥感数据获得观测地区地物分类数据，建立针对盐碱地类型的被动微波混合像元分解模型；采用带约束条件的非负最小二乘法迭代运算求解方程组，实现盐碱地被动微波混合像元分解；进一步建立分解后的被动微波遥感数据与盐碱地含盐量的关系模型，实现盐碱地特性的有效探测。

为解决本发明要解决的技术问题，给出技术方案如下：

一种基于被动微波混合像元分解技术的盐碱地特性探测方法，该方法的应用条件是被动微波遥感数据和光谱遥感数据，有如下过程：1)实现观测地区的地物分类，2)建立针对盐碱地类型的被动微波混合像元分解模型，3)被动微波混合像元分解模型求解，4）建立分解后被动微波遥感数据与盐碱地含盐量的关系模型；

所述的实现观测地区的地物分类，是对光谱遥感数据采用植被指数划分观测地区地物类型，确定观测地区水体、植被、盐碱地三类地物分类数据的结果；具体的地物分类方法是，从官方网站下载7、8和9月份分辨率光谱遥感影像，计算植被指数NDVI₇、NDVI₈和NDVI₉，其中，NDVI=（b2-b1）/（b2+b1）,式中:b2是光谱遥感数据第二个波段，b1是光谱遥感数据第一个波段；由7、8和9月连续3个月NDVI指数直方图统计结果得到阈值T₁、T₂、T₃、T₄、T₅、T₆，则满足T₁≤NDVI_7-9<T₂位置的数据定义为水体，同时满足T₂≤NDVI_8-9<T₃和T₄<NDVI₇≤T₅位置的数据定义为植被，同时满足T₂≤NDVI_8-9<T₃和T₆≤NDVI₇≤T₄位置的数据定义为盐碱地；

所述的建立针对盐碱地类型的被动微波混合像元分解模型，是从网站下载被动微波遥感数据，根据被动微波遥感数据的空间分辨率，根据观测地区的地理息，利用ArcGIS软件加载该地区shapefile矢量文件，得到该地区的被动微波遥感亮温数据，并通过ArcGIS软件将该数据转换为Grid网格数据；具体过程是：第一步是实现被动微波数据和地物分类数据配准，定位与一个被动微波混合像元匹配的地物分类信息，结合研究地区地物分类的结果，得到被动微波混合像元中该类地物的比例L；第二步被动微波混合像元模型应满足如下数学表达公式：

T_b(x,y)=T_v(x,y)L_v(x,y)+T_w(x,y)L_w(x,y)+T_s(x,y)L_s(x,y)   （3）

其中T_b代表观测地区被动微波混合像元亮温值，(x,y)代表被动微波混合像元位置；T_v代表地物分类为植被的亮温值，L_v代表与该被动微波混合像元空间位置匹配的分类数据中植被像元的比例；T_w代表地物分类为水体的亮温值，L_w代表与该被动微波混合像元空间位置匹配的分类数据中水体像元的比例；T_s代表地物分类为盐碱地的亮温值，L_s代表与该被动微波混合像元空间位置匹配的分类数据中盐碱地像元的比例；且被动微波混合像元中，所有地物分类类型出现的比例L应满足：

L_v(x,y)+L_w(x,y)+L_s(x,y)=1    （4）

第三步选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵；按公式（5）对m×n范围内的被动微波混合像元进行分解。通过构建方程组和非负最小二乘法迭代运算，求解m×n范围内被动微波混合像元分解后的各类地物的组分亮温T_c；

T_B=P_c.T_c+E   (5)

式中：c=w,v,s是地物分类的种类,即水体、植被和盐碱地，c的值为3；T_B是一个(m×n)×1的矢量，是m×n个被动微波混合像元亮温值；T_c是一个c×1的矩阵，是m×n窗口中每个被动微波混合像元对应的三类地物的组分亮温，这里m×n窗口中同类地物的组分亮温完全相同；P_c是一个(m×n)×c的矩阵，是m×n窗口中每个被动微波混合像元对应的地物分布比例；E是一个(m×n)×1的矢量，是m×n个残差数据；

所述的针对盐碱地被动微波混合像元分解模型求解，是将观测地区被动微波混合像元数据，分解成水体、植被和盐碱地三种类型的被动微波组分亮温数据；具体过程是采用迭代自组织数据分析算法（ISODATA算法）对观测地区连续一段时间的被动微波混合像元亮温进行统计分类；迭代自组织数据分析算法的分类的中心值作为水体、植被和盐碱地三类地物组分亮温初值的参考，定义各类地物分类的亮温初值X_w、X_v、X_s和各类地物分类的亮温变化阈值Y_c，确定某一地物分类的组分亮温T_c的选取范围是[X_c-Y_c,X_c+Y_c]，其中c=w,v,s代表地物分类的种类,即水体、植被和盐碱地；

所述的建立分解后被动微波遥感数据与盐碱地含盐量的关系模型，是利用C波段和X波段的被动微波混合像元分解后的亮温数据的差异性，采用差分的方法，分析其亮温差值与盐碱地含盐量间的关系；具体过程是先计算C波段和X波段的被动微波混合像元分解后的亮温数据，进一步计算亮温差值，然后通过二次多项式拟合，建立亮温差值和地面采样数据含盐量的关系模型，如下：

L=A×(T_X-T_C)²+B×(T_X-T_C)+C   （8）

其中T_C和T_X分别代表C波段和X波段的被动微波混合像元分解后的亮温数据，L代表地面采样数据的含盐量，A、B和C代表二次多项式拟合系数；通过（8）式，获得关系模型中的回归系数A、B和C，进一步可以得到观测地区盐碱地的含盐量，实现盐碱地特性的有效探测。

在2)建立针对盐碱地类型的被动微波混合像元分解模型中，所述的从网站下载被动微波遥感数据，是经过标定、大气校正、地理校正和标准化预处理的数据。

在2)建立针对盐碱地类型的被动微波混合像元分解模型中，所述的被动微波混合像元中该类地物的比例L，计算公式如（2）式：

$L=\frac{N}{{(R_{\mathrm{low}}/R_{\mathrm{high}})}^2}---\left(2\right)$

其中，L代表被动微波混合像元中该类地物的比例，N代表被动微波混合像元中某类地物的像元数量，R_high代表光谱地物分类数据的空间分辨率，R_low代表被动微波混合像元的空间分辨率。

在被动微波混合像元分解模型求解中，还设计一个目标函数R作为判断可能解优劣的标准，目标函数R定义为：

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_b(x,y)-T_v(x,y)L_v(x,y)-T_w(x,y)L_w(x,y)-T_s(x,y)L_s(x,y))}^2}<\mathrm{\&xi;}$

ξ为事先给定的阈值，它的取值大小取决于对解的精度要求。

为解决本发明要解决的技术问题，更具体的技术方案叙述如下。

一种基于被动微波混合像元分解技术的盐碱地特性探测方法，该方法的应用条件是被动微波遥感数据和光谱遥感数据，方法包括如下过程：1)实现观测地区的地物分类，2)建立针对盐碱地类型的被动微波混合像元分解模型，3)被动微波混合像元分解模型求解，4）建立分解后被动微波遥感数据与盐碱地含盐量的关系模型。具体过程如下：

（1）观测地区的地物分类

从官方网站下载中分辨率光谱遥感数据的土地覆盖类型产品，其空间分辨率为R_high，对下载的全球光谱遥感数据采用MRT软件进行图像拼接、等经纬度投影，采样方法为邻近法，椭球为WGS-84体系；根据观测地区的地理信息，在ArcGIS软件下加载该地区shapefile矢量文件，得到该地区的光谱遥感数据的土地覆盖类型数据，并通过ArcGIS软件将该数据转换为Grid网格数据。

考虑到被动微波遥感数据的空间分辨率较低，观测地区的地物数据主要分为以下几种类型：水体、盐碱地和植被。对光谱遥感数据采用植被指数划分观测地区地物类型，确定观测地区三种地物分类数据的结果。

（a）植被指数计算

由地物在各个波段的光谱反射特征可以看出，植被的反射率在近红外（b2）波段与（b1）波段的反射率差要明显大于其它地物类型。植被指数就是根据这一特征，通过比值变换，使植被信号放大，从而提取植被信息。植被指数计算有许多方法，其中NDVI最为普遍，植物叶面在可见光红光波段（Red）有很强的吸收特性，在近红外波段(NIR，Near-infrared bands)有很强的反射特性，通过这2个波段的计算可得到不同的NDV。NDVI公式表达为：

NDVI=(b2-b1)/(b2+b1)   （1）

式中:b2是光谱遥感数据第二个波段，b1是光谱遥感数据第一个波段。

(b)地物分类方法

遥感影像监测研究区盐碱地的动态变化，主要是借助于盐碱地、植物光谱的反射特征表现出来的。考虑到研究区水盐运动的规律以及农业种植制度对盐碱地和地物上覆盖植被的影响，研究选择中分辨率光谱遥感影像的采集时间为每年7、8和9月。这个时期的遥感影像更容易反映地物的细节。一般NDVI数值越大,表明植被的覆盖程度越高，植被活动越强。

令观测地区7月、8月和9月份NDVI指数为NDVI₇,NDVI₈和NDVI₉，地物分类方法如下：

步骤一：对于满足T₁≤NDVI_7-9<T₂位置的数据定义为水体，

步骤二：对于同时满足T₂≤NDVI_8-9<T₃和T₄<NDVI₇≤T₅位置的数据定义为植被，

步骤三：对于同时满足T₂≤NDVI_8-9<T₃和T₆≤NDVI₇≤T₄位置的数据定义为盐碱地；

其中的NDVI_7-9分别代表NDVI₇,NDVI₈和NDVI₉，是7月、8月和9月份NDVI指数值；其中的NDVI_8-9代表NDVI₈和NDVI₉，是8月和9月份NDVI指数值；

其中T₁～T₆是根据观测地区7、8和9月连续3个月NDVI指数直方图统计结果得到的阈值。采用上述分类方法，可得到研究地区地物分类的结果，该分类数据主要包括水体、植被和盐碱地三类。

（2）建立针对盐碱地类型的被动微波混合像元分解模型

从网站下载被动微波遥感数据，其空间分辨率为R_low，选择被动微波遥感数据的高级产品，该产品已经过标定、大气校正、地理校正和标准化等预处理。根据被动微波遥感数据的空间分辨率，利用ENVI软件对其实现等经纬度投影；根据观测地区的地理信息，利用ArcGIS软件加载该地区shapefile矢量文件，得到该地区的被动微波遥感亮温数据，并通过ArcGIS软件将该数据转换为Grid网格数据。

在ArcGIS软件下，实现被动微波数据和地物分类数据配准，可以定位与一个被动微波混合像元匹配的地物分类信息。结合研究地区地物分类的结果，可以得到被动微波混合像元中该类地物的比例。具体计算公式如（2）：

$L=\frac{N}{{(R_{\mathrm{low}}/R_{\mathrm{high}})}^2}---\left(2\right)$

其中，L代表被动微波混合像元中该类地物的比例，N代表被动微波混合像元中某类地物的像元数量，R_high代表光谱地物分类数据的空间分辨率，R_low代表被动微波混合像元的空间分辨率。

根据观测地区地物分类的结果，针对盐碱地地区的被动微波混合像元模型应满足如下数学表达公式：

T_b(x,y)=T_v(x,y)L_v(x,y)+T_w(x,y)L_w(x,y)+T_s(x,y)L_s(x,y)   （3）

其中T_b代表观测地区被动微波混合像元亮温值，(x,y)代表被动微波混合像元位置；T_v代表地物分类为植被的亮温值，L_v代表与该被动微波混合像元空间位置匹配的分类数据中植被像元的比例；T_w代表地物分类为水体的亮温值，L_w代表与该被动微波混合像元空间位置匹配的分类数据中水体像元的比例；T_s代表地物分类为盐碱地的亮温值，L_s代表与该被动微波混合像元空间位置匹配的分类数据中盐碱地像元的比例。

被动微波混合像元中，所有地物分类类型出现的比例应满足：

L_v(x,y)+L_w(x,y)+L_s(x,y)=1   （4）

选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵，并令m×n窗口中同类地物的组分亮温完全相同。已知被动微波混合像元和地物分布比例矩阵，可以按公式（5）对m×n范围内的被动微波混合像元进行分解。通过构建方程组和非负最小二乘法迭代运算，求解m×n范围内被动微波混合像元分解后的各类地物的组分亮温T_c。

T_B=P_c.T_c+E   (5)

式中：

c=w,v,s是地物分类的种类,即水体、植被和盐碱地，c的值为3；

T_B是一个(m×n)×1的矢量，是m×n个被动微波混合像元亮温值；

T_c是一个c×1的矩阵，是m×n窗口中每个被动微波混合像元对应的三类地物的组分亮温，这里m×n窗口中同类地物的组分亮温完全相同。

P_c是一个(m×n)×c的矩阵，是m×n窗口中每个被动微波混合像元对应的地物分布比例；

E是一个(m×n)×1的矢量，是m×n个残差数据。

上述的第二步骤，即（2）建立针对盐碱地类型的被动微波混合像元分解模型，是本发明的核心内容，其目的是为了求得公式（3）中的Tv，Tw和Ts,即三种地物的组分亮温；其中（3）式中已知的是Tb（卫星数据得到），Lv、Ls、Lw三个量由公式（2）计算得出，因此一个公式中有3个未知量待解，单用公式（3）求解不到未知量，因此引入公式（5）构成分解模型方程组。公式（5）可以看成是由m×n个公式（3）构成，即求解公式（5）的同时，公式（3）中的未知量都会得到，最终可以得到m×n窗口中每个被动微波混合像元对应的三类地物的组分亮温。由于解可能有很多个，求解方法有很多，所以引入下面的第三步骤，即，（3）被动微波混合像元分解模型求解，利用公式（6）和（7），通过亮温初值选取，限定公式（5）方程组中解的范围，都是为了最终得到公式（5）中m×n窗口中每个被动微波混合像元对应的三类地物的组分亮温。

（3）被动微波混合像元分解模型求解

采用ISODATA算法（迭代自组织数据分析算法）对观测地区连续一段时间的被动微波混合像元亮温进行统计分类。ISODATA算法的分类的中心值可以作为水体、盐碱地和植被三类地物组分亮温初值的参考。定义各类地物分类的亮温初值X_w、X_v、X_s，同时定义各类地物分类的亮温变化阈值Y_c，确定某一地物分类的组分亮温T_c的选取范围是[X_c-Y_c,X_c+Y_c]，其中c=w,v,s是地物分类的种类,即水体、植被和盐碱地。

根据相关研究，三类地物的被动微波组分亮温基本满足下式：

0<T_W<T_S<T_V   （6）

此外，还需设计一个目标函数作为判断可能解优劣的标准，这里选取R作为目标函数，ξ为事先给定的阈值，它的取值大小取决于对解的精度要求。目标函数R定义为：

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_b(x,y)-T_v(x,y)L_v(x,y)-T_w(x,y)L_w(x,y)-T_s(x,y)L_s(x,y))}^2}<\mathrm{\&xi;}---\left(7\right)$

根据目标函数调整搜索窗口尺寸m×n，合理选择各类组分的亮温初值X_w、X_v和X_s，代入方程（5）中，在（6）式、（7）式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解方程组。由于初值的选取范围考虑到不同地物分类的亮温变化空间，且存在多个条件限制方程组的求解范围，因此采用本发明提出的方法可以得到较好的最优解。

（4）建立分解后被动微波遥感数据与盐碱地含盐量的关系模型

微波遥感可以提供多频率、多角度的观测数据，还可以提供不同极化方式的数据，多极化特点是微波遥感所特有的。不同频率和极化的微波辐射对地表参数的响应不一样，可以利用这些差异对盐碱地特性进行检测。本发明利用C波段和X波段的被动微波混合像元分解后的亮温数据的差异性，采用差分的方法，分析其亮温差值与盐碱地含盐量间的关系。

具体实施方式是先计算C波段和X波段的被动微波混合像元分解后的亮温数据，进一步计算亮温差值，然后通过二次多项式拟合，建立亮温差值和地面采样数据含盐量的关系模型，如下：

L=A×(T_X-T_C)²+B×(T_X-T_C)+C   （8）

其中T_C和T_X分别代表C波段和X波段的被动微波混合像元分解后的亮温数据，L代表地面采样数据的含盐量，A、B和C代表二次多项式拟合系数。通过（8）式，获得关系模型中的回归系数A、B和C，进一步可以得到观测地区盐碱地的含盐量，实现盐碱地特性的有效探测。

本发明的有益效果：

本发明利用光学和微波遥感在盐碱地识别上的优势，紧密结合多时相的光学与被动微波遥感观测数据，提出一种基于被动微波混合像元分解技术的盐碱地特性探测方法。通过建立分解后被动微波遥感数据与盐碱地含盐量的关系模型，实现观测地区盐碱地特性的有效探测，为深入研究苏打盐碱地水盐特性与波谱信息的相关性，及苏打盐碱土的导电特性与微波辐射/散射信息的相关性，掌握苏打盐碱地的基本形态、时空动态变化过程、演变规律和地理分布特性，及苏打盐碱地的改良与合理利用，提供理论方法和技术手段。

附图说明

图1是本发明实施例1的吉林省西部地区的地物分类的结果。

图2是本发明实施例1的MWRI10.65GHz水平极化被动微波遥感数据。

图3是本发明实施例1的MWRI18.7GHz水平极化被动微波遥感数据。

图4是本发明实施例1的MWRI10.65GHz水平极化被动微波混合像元分解数据。

图5是本发明实施例1的MWRI18.7GHz水平极化被动微波混合像元分解数据。

图6是本发明实施例1的2011年8月吉林省西部地区盐碱地含盐量特性分布图。

具体实施方式

实施例1：

研究区位于吉林省西部，属于松嫩平原的西南部。该区所处的经纬度范围：121°38′-126°11′E，43°59′-46°18′N，区域包括含大安市、镇赉县、通榆县等八个区域，面积约4.7万km²。虽然在整个区域内都分布有盐碱化土地,但各地的面积、发育程度都不尽相同。本发明利用MWRI被动微波遥感数据和MODIS的8日合成地表反射率数据产品，结合提出的盐碱地被动微波混合像元分解方法，实现了2011年8月中国吉林省西部地区盐碱地特性探测。

具体包括以下步骤：

（1）观测地区的地物分类：

在MODIS官方网站下载中分辨率MODIS光谱遥感数据的8日合成地表反射率数据MOD09A1，其空间分辨率为500m，对下载的MOD09A1数据采用MRT软件进行图像拼接、等经纬度投影，采样方法为邻近法，椭球为WGS-84体系；根据观测地区的地理信息，在ArcGIS软件下加载该地区shapefile矢量文件，得到该地区的MODIS的8日地表反射率数据，并通过ArcGIS软件将该数据转换为Grid网格数据,空间分辨率为0.0005°×0.0005°（500m×500m）。

观测地区的地物数据主要分为以下几种地物分类类型：水体、盐碱地和植被。对MOD09A1数据采用植被指数划分观测地区地物类型，确定观测地区三种地物分类数据的结果。

（a）植被指数计算

地物在各个波段的光谱反射特征可以看出植被的反射率在近红外（b2）波段与（b1）波段的反射率差要明显大于其它地物类型。植被指数就是根据这一特征,通过比值变换,使植被信号放大,从而提取植被信息。植被指数计算有许多方法，其中NDVI.最为普遍，植物叶面在可见光红光波段（Red）有很强的吸收特性，在近红外波段(NIR，Near-infrared bands)有很强的反射特性，通过这2个波段的计算可得到不同的NDV。NDVI公式表达为：

NDVI=(b2-b1)/(b2+b1)   （1）

式中:b2是MODIS的8日地表反射率数据MOD09A1的第二个波段，b1是MODIS的8日地表反射率数据MOD09A1的第一个波段。

(b)地物分类方法

遥感影像监测研究区盐碱地的动态变化，主要是借助于盐碱地、植物光谱的反射特征表现出来的。为了充分考虑研究区水盐运动的规律以及农业种植制度对盐碱地和地物上覆盖植被的影响，研究选择中分辨率MODIS光谱遥感影像的采集时间为2011年7、8和9月。这个时期的遥感影像更容易反映地物的细节。一般NDVI数值越大,表明植被的覆盖程度越高，植被活动越强。

令吉林省西部地区7月、8月和9月份NDVI指数为NDVI₇,NDVI₈和NDVI₉，结合吉林省西部地区7、8和9月连续3个月NDVI直方图统计结果得到分类阈值，地物分类方法如下：

步骤一：对于满足-1≤NDVI_7-9<0位置的数据定义为水体；

步骤二：对于同时满足0≤NDVI_8-9<1和0.4<NDVI₇≤1位置的数据定义为植被；

步骤三：对于同时满足0≤NDVI_8-9<1和0≤NDVI₇≤0.4位置的数据定义为盐碱地；

其中NDVI_7-9分别代表NDVI₇,NDVI₈和NDVI₉，是7月、8月和9月份NDVI指数值；NDVI_8-9代表NDVI₈和NDVI₉，是8月和9月份NDVI指数值。

采用上述分类方法，可得到研究地区地物分类的结果。如图1所示，该地物分类数据主要包括水体、植被和盐碱地三类，地物分类数据的空间分辨率为0.005°×0.005°（500m×500m）。

（2）建立针对盐碱地类型的被动微波混合像元分解模型

风云三号B星上装载的微波成像仪(MWRI）为我国第一个星载微波遥感仪器,扫描方式为圆锥扫描，其设计频率为10.65-150GHz,其中150GHz为试验通道。每个频率都有垂直和水平两种不同极化模式，这些频率的遥感成像能提供全天候、全天时地表温度、土壤水分、洪涝干旱、积雪深度、台风结构、大气含水量等丰富的信息。

从网站下载MWRI被动微波遥感数据，其空间分辨率为10km，选择被动微波遥感数据的高级产品，数据采集时间是2011年8月3日，其C和X频段分别选取其10.65GHz水平极化被动微波遥感数据和18.7GHz水平极化被动微波遥感数据(图2和图3所示)。该产品已经过标定、大气校正、地理校正和标准化等预处理。根据被动微波遥感数据的空间分辨率，利用ENVI软件对其实现等经纬度投影；根据吉林省西部地区的地理信息，利用ArcGIS软件加载该地区shapefile矢量文件，得到该地区的MWRI被动微波遥感亮温数据。MWRI被动微波数据的空间分辨率为0.1°×0.1°（10km×10km）。

在ArcGIS软件下，实现MWRI被动微波数据和地物分类数据的配准，可以定位与一个MWRI被动微波混合像元匹配的地物分类信息，即每个MWRI数据对应20×20个地物分类数据。结合研究地区地物分类的结果，可以得到MWRI被动微波混合像元中该类地物的比例。具体计算公式如（2）：

$L=\frac{N}{{(R_{\mathrm{low}}/R_{\mathrm{high}})}^2}---\left(2\right)$

其中，L代表MWRI被动微波混合像元中该类地物的比例，N代表MWRI被动微波混合像元中某类地物的像元数量，R_high代表地物分类数据的空间分辨率为500m;R_low代表MWRI被动微波混合像元的空间分辨率为10km。

根据观测地区地物分类的结果，针对盐碱地地区的被动微波混合像元模型应满足如下数学表达公式：

T_b(x,y)=T_v(x,y)L_v(x,y)+T_w(x,y)L_w(x,y)+T_s(x,y)L_s(x,y)   （3）

其中T_b代表观测地区被动微波混合像元亮温值，(x,y)代表被动微波混合像元位置；T_v代表地物分类为植被的亮温值，L_v代表与该被动微波混合像元空间位置匹配的分类数据中植被像元的比例；T_w代表地物分类为水体的亮温值，L_w代表与该被动微波混合像元空间位置匹配的分类数据中水体像元的比例；T_s代表地物分类为盐碱地的亮温值，L_s代表与该被动微波混合像元空间位置匹配的分类数据中盐碱地像元的比例。

MWRI被动微波混合像元中，所有地物分类类型出现的比例应满足：

L_v(x,y)+L_w(x,y)+L_s(x,y)=1   （4）

选取2×2范围的被动微波混合像元构成一个搜索窗口，记录2×2窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵，并令2×2窗口中同类地物的组分亮温完全相同。已知被动微波混合像元和地物分布比例矩阵，可以按公式（5）对2×2范围内的被动微波混合像元进行分解，通过构建方程组和非负最小二乘法迭代运算，求解2×2范围内被动微波混合像元分解后的各类地物的组分亮温T_c。

T_B=P_c.T_c+E   (5)

式中：

c=v,w,s是地物分类的种类,即水体、植被和盐碱地，c的值为3；

T_B是一个(2×2)×1的矢量，是2×2个被动微波混合像元亮温值；

T_c是一个3×1的矩阵，是2×2窗口中每个被动微波混合像元对应的三类地物的组分亮温，这里2×2窗口中同类地物的组分亮温完全相同；

P_c是一个(2×2)×3的矩阵，是2×2窗口中每个被动微波混合像元对应的地物分布比例；

E是一个(2×2)×1的矢量，是2×2个残差数据。

（3）被动微波混合像元分解模型求解

采用ISODATA算法（迭代自组织数据分析算法）对观测地区连续一段时间的被动微波混合像元亮温进行统计分类。ISODATA算法的分类的中心值可以作为水体、盐碱地和植被三类地物亮温初值的参考。

对于10.65GHz水平极化被动微波数据，三类地物的亮温初始值为X_w=237k、X_v=268k和X_s=254k。同时定义各类地物分类的亮温变化阈值Y_c，确定某一地物分类的组分亮温T_c的选取范围是[X_c-Y_c,X_c+Y_c]，其中c=v,w,s是地物分类的种类,即水体、植被和盐碱地。定义水体的亮温变化阈值Y_w为5k，其组分亮温T_w的选取范围是237±5k;植被的亮温变化阈值Y_v为7k，其组分亮温T_v的选取范围是268±7k；盐碱地的亮温变化阈值Y_S为6k，其组分亮温T_s的选取范围是254±6k。

对于18.7GHz GHz水平极化被动微波数据，三类地物的亮温初始值为X_w=239k、X_v=274k和X_s=262k。同时定义各类地物分类的亮温变化阈值Y_c，确定某一地物分类的组分亮温T_c的选取范围是[X_c-Y_c,X_c+Y_c]，其中c=v,w,s是地物分类的种类,即水体、植被和盐碱地。定义水体的亮温变化阈值Y_w为6k，其组分亮温T_w的选取范围是239±6k;植被的亮温变化阈值Y_v为6k，其组分亮温T_v的选取范围是274±6k；盐碱地的亮温变化阈值Y_S为5k，其组分亮温T_s的选取范围是262±5k。

根据相关研究，三类地物的被动微波组分亮温基本满足下式：

0<T_W<T_S<T_V   （6）

此外，还需设计一个目标函数作为判断可能解优劣的标准，这里选取R作为目标函数，ξ为事先给定的阈值为0.2K，它的取值大小取决于对解的精度要求。目标函数R定义为：

$R=\sqrt{\underset{x=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{2}{\mathrm{\&Sigma;}}}{(T_b(x,y)-T_v(x,y)L_v(x,y)-T_w(x,y)L_w(x,y)-T_s(x,y)L_s(x,y))}^2}<\mathrm{\&xi;}---\left(7\right)$

在（6）式、（7）式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得（5）式最优解，应用matlab软件中的fsolve函数求解方程组。由于初值的选取范围考虑到不同地物分类的亮温变化空间，且存在多个条件限制方程组的求解范围，因此采用本发明提出的方法可以得到较好的最优解。图4和图5分别是10.65GHz和18.7GHz水平极化被动微波混合像元分解数据。

（4）建立分解后被动微波遥感数据与盐碱地含盐量的关系模型

微波遥感可以提供多频率、多角度的观测数据，还可以提供不同极化方式的数据，多极化特点是微波遥感所特有的。不同频率和极化的微波辐射对地表参数的响应不一样，可以利用这些差异对盐碱地特性进行检测。本发明利用C波段（10.65GHz水平极化）和X波段（18.7GHz水平极化）的被动微波混合像元分解后的亮温数据的差异性，采用差分的方法，分析其亮温差值与盐碱地含盐量间的关系。

具体实施方式是先计算10.65GHz和18.7GHz水平极化的被动微波混合像元分解后的亮温数据，进一步计算亮温差值，然后通过二次多项式拟合，建立亮温差值和地面采样数据含盐量的关系模型，如下：

L=A×(T_X-T_C)²+B×(T_X-T_C)+C   （8）

其中T_C和T_X分别代表10.65GHz和18.7GHz水平极化的被动微波混合像元分解后的亮温数据，L代表地面实测采样数据的含盐量，如表1所示；A、B和C代表二次多项式拟合系数。通过（8）式，获得关系模型中的拟合系数，其中A的值为0.0015，B的值为-0.0590，C的值为1.5699。进一步将拟合系数代入公式（8），可以得到观测地区盐碱地的含盐量（图6所示），实现盐碱地特性有效探测。

表1地面实测采样数据的含盐量

标号	经度	纬度	含盐量（g/kg）
				1	123.741550	45.630400	8.4
2	123.691967	45.637667	18.96552
				3	123.710983	45.456917	8.282759
4	123.596017	45.456133	16.72414
				5	123.323433	45.391433	10.22759
6	123.076633	45.368950	10.5
				7	122.845417	45.331817	10.24138
8	122.876500	45.189667	9.651724
				9	123.047183	44.909200	11.69655
10	123.219467	44.778033	13.15862
				11	123.309500	44.772500	12.8
12	123.458817	44.785817	8.341379
				13	123.724850	44.886383	9.3
14	123.814167	44.938550	12.01379
				15	123.973833	44.972200	9.87931
16	124.319467	45.026267	9.506897
				17	123.981767	45.599517	10.14483
18	123.919300	45.508883	12.44483
				19	124.181867	45.408083	13.68966
20	123.511817	46.145333	8.182759

实验结果：

通过选择具有典型盐碱地特征的中国吉林省西部地区作为研究区域，采用专利提出的一种基于被动微波混合像元分解技术的盐碱地特性探测方法，将光学和被动微波遥感技术相结合，共同实现观测地区盐碱地特性的有效探测。通过MODIS光谱遥感数据获得观测地区地物分类数据，可以实现粗略定位盐碱地分布的地理位置。建立针对盐碱地类型的MWRI被动微波混合像元分解模型，采用带约束条件的非负最小二乘法迭代运算求解方程组，实现盐碱地被动微波混合像元分解。结合地面实测盐碱地含盐量数据，进一步建立分解后的被动微波遥感数据与盐碱地含盐量的关系模型，实现盐碱地特性的有效探测。实验结果表明进入夏季7月和8月，夏季雨水多而集中，大量可溶性盐随水渗到下层或流走，进入“脱盐”季节，不同地区的盐碱地的含盐量比较接近。提出的方法能够快速、较准确地给出盐碱化程度的判定信息，为被动微波遥感数据在盐碱地研究领域的广泛深入应用提供理论与数据支持，对合理开发利用盐碱地资源,保护生态环境，具有重要的意义。

Claims (4)

1.一种基于被动微波混合像元分解技术的盐碱地特性探测方法，该方法的应用条件是被动微波遥感数据和光谱遥感数据，有如下过程：1)实现观测地区的地物分类，2)建立针对盐碱地类型的被动微波混合像元分解模型，3)被动微波混合像元分解模型求解，4)建立分解后被动微波遥感数据与盐碱地含盐量的关系模型；

①所述的实现观测地区的地物分类，是对光谱遥感数据采用植被指数划分观测地区地物类型，确定观测地区水体、植被、盐碱地三类地物分类数据的结果；具体的地物分类方法是，从官方网站下载7、8和9月份分辨率光谱遥感影像，计算植被指数NDVI₇、NDVI₈和NDVI₉，其中，NDVI＝(b2-b1)/(b2+b1),式中:b2是光谱遥感数据第二个波段，b1是光谱遥感数据第一个波段；由7、8和9月连续3个月NDVI指数直方图统计结果得到阈值T₁、T₂、T₃、T₄、T₅、T₆，则满足T₁≤NDVI_7-9<T₂位置的数据定义为水体，同时满足T₂≤NDVI_8-9<T₃和T₄<NDVI₇≤T₅位置的数据定义为植被，同时满足T₂≤NDVI_8-9<T₃和T₆≤NDVI₇≤T₄位置的数据定义为盐碱地；

②所述的建立针对盐碱地类型的被动微波混合像元分解模型，是从网站下载被动微波遥感数据，根据被动微波遥感数据的空间分辨率，根据观测地区的地理信息，利用ArcGIS软件加载该地区shapefile矢量文件，得到该地区的被动微波遥感亮温数据，并通过ArcGIS软件将该数据转换为Grid网格数据；具体过程是：第一步是实现被动微波数据和地物分类数据配准，定位与一个被动微波混合像元匹配的地物分类信息，结合研究地区地物分类的结果，得到被动微波混合像元中该类地物的比例L；第二步被动微波混合像元模型应满足如下数学表达公式：

T_b(x,y)＝T_v(x,y)L_v(x,y)+T_w(x,y)L_w(x,y)+T_s(x,y)L_s(x,y)   (3)

其中T_b代表观测地区被动微波混合像元亮温值，(x,y)代表被动微波混合像元位置；T_v代表地物分类为植被的亮温值，L_v代表与该被动微波混合像元空间位置匹配的分类数据中植被像元的比例；T_w代表地物分类为水体的亮温值，L_w代表与该被动微波混合像元空间位置匹配的分类数据中水体像元的比例；T_s代表地物分类为盐碱地的亮温值，L_s代表与该被动微波混合像元空间位置匹配的分类数据中盐碱地像元的比例；且被动微波混合像元中，所有地物分类类型出现的比例L应满足：

L_v(x,y)+L_w(x,y)+L_s(x,y)＝1   (4)

第三步选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵；按公式(5)对m×n范围内的被动微波混合像元进行分解；通过构建方程组和非负最小二乘法迭代运算，求解m×n范围内被动微波混合像元分解后的各类地物的组分亮温T_c；

T_B＝P_c.T_c+E   (5)

式中：c＝w,v,s是地物分类的种类,即水体、植被和盐碱地，c的值为3；T_B是一个(m×n)×1的矢量，是m×n个被动微波混合像元亮温值；T_c是一个c×1的矩阵，是m×n窗口中每个被动微波混合像元对应的三类地物的组分亮温，这里m×n窗口中同类地物的组分亮温完全相同；P_c是一个(m×n)×c的矩阵，是m×n窗口中每个被动微波混合像元对应的地物分布比例；E是一个(m×n)×1的矢量，是m×n个残差数据；

③所述的针对盐碱地被动微波混合像元分解模型求解，是将观测地区被动微波混合像元数据，分解成水体、植被和盐碱地三种类型的被动微波组分亮温数据；具体过程是采用迭代自组织数据分析算法(ISODATA算法)对观测地区连续一段时间的被动微波混合像元亮温进行统计分类；迭代自组织数据分析算法的分类的中心值作为水体、植被和盐碱地三类地物组分亮温初值的参考，定义各类地物分类的亮温初值X_w、X_v、X_s和各类地物分类的亮温变化阈值Y_c，确定某一地物分类的组分亮温T_c的选取范围是[X_c-Y_c,X_c+Y_c]，其中c＝w,v,s代表地物分类的种类,即水体、植被和盐碱地；

④所述的建立分解后被动微波遥感数据与盐碱地含盐量的关系模型，是利用C波段和X波段的被动微波混合像元分解后的亮温数据的差异性，采用差分的方法，分析其亮温差值与盐碱地含盐量间的关系；具体过程是先计算C波段和X波段的被动微波混合像元分解后的亮温数据，进一步计算亮温差值，然后通过二次多项式拟合，建立亮温差值和地面采样数据含盐量的关系模型，如下：

H＝A×(T_X-T_C)²+B×(T_X-T_C)+C   (8)

其中T_C和T_X分别代表C波段和X波段的被动微波混合像元分解后的亮温数据，H代表地面采样数据的含盐量，A、B和C代表二次多项式拟合系数；通过(8)式，获得关系模型中的回归系数A、B和C，进一步可以得到观测地区盐碱地的含盐量，实现盐碱地特性的有效探测。

2.根据权利要求1所述的基于被动微波混合像元分解技术的盐碱地特性探测方法，其特征是，所述的从网站下载被动微波遥感数据，是经过标定、大气校正、地理校正和标准化预处理的数据。

3.根据权利要求1所述的基于被动微波混合像元分解技术的盐碱地特性探测方法，其特征是，所述的被动微波混合像元中该类地物的比例L，计算公式如(2)式：

$L=\frac{N}{{(R_{low}/R_{high})}^2}---\left(2\right)$

其中，L代表被动微波混合像元中该类地物的比例，N代表被动微波混合像元中某类地物的像元数量，R_high代表光谱地物分类数据的空间分辨率，R_low代表被动微波混合像元的空间分辨率。

4.根据权利要求1所述的基于被动微波混合像元分解技术的盐碱地特性探测方法，其特征是，在被动微波混合像元分解模型求解中，还设计一个目标函数R作为判断可能解优劣的标准，目标函数R定义为：

$R=\sqrt{\underset{x=1}{\overset{m}{\&Sigma;}}\underset{y=1}{\overset{n}{\&Sigma;}}{\left(T_b\left(x,y\right)-T_v\left(x,y\right)L_v\left(x,y\right)-T_w\left(x,y\right)L_w\left(x,y\right)-T_s\left(x,y\right)L_s\left(x,y\right)\right)}^2}<\mathrm{\&xi;}$

ξ为事先给定的阈值，它的取值大小取决于对解的精度要求。