CN102608592B - 基于五种地物分类信息的积雪被动微波混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，提供了一种基于五种地物分类信息的积雪被动微波混合像元分解方法。本发明根据影响积雪辐射特性的典型五种下垫面类型，通过获得观测地区中较高空间分辨率的地物分类数据，利用积雪在不同下垫面微波辐射之间差异特性，选择微波天线增益函数和采样率，建立积雪微波混合像元分解模型，采用具有约束条件的最小二乘法迭代计算求解欠定性方程组，实现积雪被动微波混合像元分解。本发明可以有效地解决积雪被动微波混合像元问题，改进积雪参数反演精度，在气候和水文研究以及积雪灾害评估等领域中具有重要应用价值。

Description

基于五种地物分类信息的积雪被动微波混合像元分解方法

技术领域

本发明属于遥感图像处理的技术领域，具体涉及一种基于五种地物分类信息的积雪被动微波混合像元分解方法，可解决由于积雪微波混合像元导致积雪参数反演精度较低的问题。

背景技术

积雪的变化与特征是气候研究、天气预报和水资源管理的重要参数(参见下面列出的参考文献1)。积雪遥感中最为常用的波段是可见光、近红外和微波，其中，可见光和近红外主要用于提取积雪覆盖范围，它们最大的弱点是不能用于反演雪深和雪水当量。微波在积雪遥感中处于不可缺少的位置，它不仅能够全天观测积雪，也能够穿透大部分积雪层，从而探测到雪深和雪水当量信息。然而，合成孔径雷达(SAR)卫星数据受到重复访问周期较长、山区地形影响以及费用较高等多方面制约，很难应用于研究。被动微波遥感具有很高的时间分辨率，能够迅速覆盖全球而且数据资料免费，因此，它在监测全球和大陆尺度的积雪时空变化中作用尤为突出(参见下面列出的参考文献2、3)。自20世纪70年代以来，几代微波辐射计已经获取了全球积雪深度的较为可信的资料，各种成熟的积雪遥感数据产品也已经被应用于气候和水文研究以及灾害评估等领域(参见下面列出的参考文献4)。

被动微波遥感反演积雪深度的基础是雪粒子在毫米波段对来自下垫面微波信号的散射大于在厘米波的散射，且对雪盖深度敏感，因此可以利用多频段被动微波遥感信息探测雪盖深度(参见下面列出的参考文献5)。被动微波遥感是提供全球尺度雪水当量研究的唯一有效途径，但是由于被动微波遥感数据的空间分辨率较低(≥10km)，其天线波束覆盖区大多为复合像元，使得积雪深度的反演精度还不能满足实际需求(参见下面列出的参考文献6)。

目前，国内外提出的被动微波混合像元分解方法多数是针对沿海区域微波混合像元的分解(参见下面列出的参考文献7、8、9)，对于足印内只含有海洋和陆地两种地表类型的情况比较适用，基本实现了陆地和水体的微波混合像元分解。然而，当天线足印内显著包含两种以上地物时，采用目前提出的被动微波混合像元分解模型均会产生较大的误差。

被动微波辐射计获取的积雪亮温是代表了一定尺度的混合像元综合亮温，而不同下垫面的积雪微波辐射特性具有较大的差异，因此导致了被动微波遥感反演积雪深度的精度难以提高。为了改进积雪参数反演精度，需要针对积雪下垫面的不同分类建立积雪微波混合像元分解模型，从而解决积雪微波混合像元问题，为更多的研究和应用领域服务。

与本发明相关的现有技术有如下参考文献：

1.Rees W G，2006.Remote Sensing of Snow and Ice[M].Boca Raton，FL：CRCPress，Taylor&Francis，1-285.

2.Foster，J.L.，Hall，D.K.，Kelly，R.E.J.，2009.Seasonal snow extent and snow mass inSouth America using SMMR and SSM/I passive microwave data(1979-2006)[J]Remote Sensing of Environment，113(2)：291-305.

3.Evora Noёl Dacruz，Tapsoba Dominique，De Seve Danielle，2008.Combiningartificial neural network models，geostatistics，and passive microwave data forsnow water equivalent retrieval and mapping[J].IEEE Transactions on Geoscienceand Remote Sensing，46(7)：1925-1939.

4.李新，车涛.积雪被动微波遥感研究进展[J].冰川冻土.2007，29(3)：487-496.

5.Savoie Matthew H.，Armstrong，Richard L.，Brodzik，Mary J.，2009.Atmosphericcorrections for I mproved satellite passive microwave snow cover retrievals overthe Tibet Plateau[J]Remote Sensing of Environment，113(12)：2661-2669.

6.刘宝康，冯蜀青，杜玉娥.积雪被动微波遥感研究进展与前景展望[J].草业科学.2009，26(11)：37-43.

7.Tim Bellerby，Malcolm Taberner，Andrea Wilmshurst，1998.Retrieval of Land andSea Brightness Temperatures from Mixed Coastal Pixels in Passive MicrowaveData[J].IEEE Trans.Geosci.RemoteSensing，36(6)：1844-1851.

8.MaaβNina，Kaleschke Lars，2010.Improving passive microwave sea iceconcentration algorithms for coastal areas：applications to the Baltic Sea[J].TellusSeries A Dynamic Meteorology and Oceanography，62(4)：393-410.

9.Gu Lingjia，Zhao Kai，Zhang Shuang，Zheng Xingming，2011.An AMSR-E DataUnmixing Method for Monitoring Flood and Waterlogging Disaster[J].ChineseGeographical Science.2011，21(6)：666-675。

发明内容

本发明要解决的技术问题在于提供了一种基于五种地物分类信息的积雪被动微波混合像元分解方法，可以有效地解决积雪被动微波混合像元问题，改进积雪参数反演精度，在气候和水文研究以及灾害评估等领域中具有重要应用价值。

本发明根据影响积雪辐射特性的典型五种下垫面类型，通过获得观测地区较高空间分辨率的地物分类数据，利用积雪在不同下垫面微波辐射之间差异特性，选择微波天线增益函数和采样率，建立积雪微波混合像元分解模型，采用具有约束条件的最小二乘法迭代计算求解欠定性方程组，实现积雪被动微波混合像元分解。本方法可以实现全球范围内积雪微波混合像元分解，通过计算来获得五种下垫面的积雪辐射亮温，提高后期积雪参数反演精度。

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

一种基于五种地物分类信息的积雪被动微波混合像元分解方法，该方法的应用条件是冬季积雪被动微波遥感数据，其特征在于，方法包括如下过程：1)确定观测地区冬季积雪情况下的五种的地物分类数据，2)确定微波天线增益函数和采样率，3)建立被动微波积雪混合像元模型，4)积雪被动微波混合像元分解模型求解；

所述的确定观测地区冬季积雪情况下的五种的地物分类数据，是将已有的观测地区的地物分类数据产品重新分类为水体、草地、林地、农田和裸土五种下垫面信息，确定观测地区冬季积雪情况下的五种地物分类数据的结果L；

所述的确定微波天线增益函数和采样率，是根据被动微波辐射计的特性，选择被动微波数据中不同频率对应的椭圆足印a和b，设计微波天线增益函数水平采样率Δx和垂直采样率Δy，结合地物分类数据的结果L，利用公式(12)求得系数α：

$\mathrm{\α}=\underset{x=-a}{\overset{a}{\mathrm{\Σ}}}\left(\underset{y=-\mathrm{bf}\left(x\right)}{\overset{\mathrm{bf}\left(x\right)}{\mathrm{\Σ}}}\stackrel{\‾}{G}(x,y)L(x,y)\mathrm{\Δy}\right)\mathrm{\Δx}---\left(12\right)$

α是天线增益函数

和对应足印内地物分类数据的结果L的卷积，(x，y)是被动微波数据的位置；

所述的建立被动微波积雪混合像元模型，是将地物分类数据和微波天线增益函数及采样率的计算结果代入积雪被动微波混合像元模型，确定积雪被动微波混合像元分解模型求解方法；具体的积雪被动微波混合像元分解模型为：

$T_B=\∫\∫\stackrel{\‾}{G}(x,y)\·\{L_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+L_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+L_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+L_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+L_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)\}\mathrm{dxdy}---\left(16\right)$

$={\mathrm{\α}}_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+{\mathrm{\α}}_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+{\mathrm{\α}}_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+{\mathrm{\α}}_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+{\mathrm{\α}}_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)$

其中：

${\mathrm{\α}}_c(x,y)=\∫\∫\stackrel{\‾}{G}(x,y)L_c(x,y)\mathrm{dxdy}$ c＝land，water，forrest，grass，crop  (17)

公式(16)中，T_B代表积雪混合像元亮温值，(x，y)代表混合像元空间位置，即，被动微波数据的位置；T_land、T_water、T_forest、T_grass和T_crop分别代表下垫面为裸土、水体、林地、草地和农田的亮温值，L_land、L_water、L_forest、L_grass和L_crop分别代表与混合像元空间位置(x，y)匹配的分类数据中裸土、水体、林地、草地和农田像元的比例；公式(17)的结果由公式(12)计算得到；

所述的积雪被动微波混合像元分解模型求解，是最终将积雪被动微波混合像元数据分解成五种类型的被动微波组分亮温数据；具体过程是利用建立的积雪被动微波混合像元分解模型，结合被动微波遥感数据，建立联立方程组，对m×n范围内的被动微波混合像元进行分解，通过构建方程组和约束性最小二乘法求解m×n范围内被动微波混合像元分解后的五种地物的组分亮温T_c。

上述的积雪被动微波混合像元分解模型求解，其过程可以进一步说明如下：

所述的积雪被动微波混合像元分解模型求解方法，

首先选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物分类的出现比例，构成地物分类分布比例矩阵，

T_BN＝P_CT_C+E                                     (18)

式中：c是地物分类的种类，c的值为5；T_BN是一个(m×n)×1的矢量，是m×n个被动微波混合像元亮温值；T_c是一个(m×n×c)×1的矩阵，是m×n窗口中每个被动微波混合像元对应的五种地物分类的组分亮温；P_c是一个(m×n)×(m×n×c)的矩阵，是m×n窗口中每个被动微波混合像元对应的五种地物分类的分布与天线方向图的卷积；E是一个(m×n)×1的矢量，是m×n个残差数据；m选取1～10，n选取1～10；

其次求解m×n范围内被动微波混合像元分解后的五种地物的组分亮温T_c；

m×n方程组中有(m×n)×c个未知量，该方程组属于欠定方程组，有无穷多组解；对每个下垫面亮温值设定其对应的取值区间，在该区间内，任意一点都对应着一个可能的解，从而构成完整的解空间；

通过地基和星载微波辐射计对不同下垫面类型进行观测，得到不同下垫面的亮温变化区间，采用k-means聚类算法对观测地区连续一段时间的被动微波混合像元亮温进行统计分类；k-means聚类后的聚类中心值作为不同下垫面亮温初值的参考；将地基微波数据建立的数据库中的数据、被动微波遥感数据统计得到的各类下垫面的亮温相结合，定义各类下垫面的亮温初值S_grass、S_water、S_forest、S_crop、S_land，同时定义各类下垫面的亮温变化阈值L_c，确定某一下垫面的初值亮温T_c的选取范围是[S_c-L_c，S_c+L_c]，其中c代表grass，water，land，forest，crop；

不同下垫面的微波亮温变化基本满足下式：

0＜T_water＜T_land＜T_crop＜T_forest                        (19)

草地亮温T_grass变化与草的高矮和土壤湿度有关，需要结合具体观测地区进行分析；

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

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_B(x,y)-\underset{c=1}{\overset{5}{\mathrm{\&Sigma;}}}{P}_c(x,y)\&CenterDot;T_c(x,y))}^2}<\mathrm{\&xi;}---\left(20\right)$

根据目标函数调整搜索窗口尺寸m×n，合理选择各类组分的亮温初值S_grass、S_water、S_forest、S_crop、S_land，代入方程(18)中，在(19)式、(20)式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解欠定性方程组。

本发明更细致具体的技术方案叙述为以下4个步骤：

(1)地物分类数据的获取

积雪微波混合像元分解方法，需要已知观测地区的下垫面地物分类信息。积雪微波亮温值受其下垫面影响较大，通过地基实验发现林地上积雪亮温受树干的影响较大，草地上积雪亮温主要受草的散射影响，农田上积雪亮温受农田几何形状的影响较大，冰上的积雪亮温随雪深增加亮温值增加，裸土上的积雪亮温随雪深增加亮温值降低，影响积雪亮温的其他下垫面类型基本都包含在这五种典型下垫面类型中。因此，积雪下垫面分类可以定义为水体(冰)、草地、林地、农田和裸土，对同一类型的地物不需要更加细致的划分。

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

目前已有分类数据产品的土地覆盖类型多于研究所需的五种类型，需要将地物分类数据产品按照这五种类型重新分类。重新分类的总体原则如下：

(a)原分类数据中定义为农田、阔叶作物或谷类作物的类型统一重新定义为农田类型；

(b)原分类数据中定义为阔叶林、针叶林、混合林或灌丛的类型统一重新定义为林地类型；

(c)原分类数据中定义为草原、草地、植被类型统一重新定义为草地类型；

(d)原分类数据中定义为水体类型仍旧定义为水体类型；

(e)原分类数据中定义为城建用地、荒漠或荒地的类型统一重新定义为裸土类型；

采用上述分类方法，可得到研究地区地物分类的结果，该分类数据主要包括水体(包括冰)、草地、林地、农田和裸土五类。

(2)微波天线增益函数和采样率的选择

卫星上辐射计天线测得的亮温是地表实际亮温和天线增益函数的积分：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)T_b(x,y)\mathrm{dxdy}---\left(1\right)$

公式(1)中，T_B是辐射计天线测量的亮温值，T_b是在位置(x，y)处的实际亮温值，

是天线增益函数(允许仪器在测试期间扫描移动)，积分覆盖整个足印的椭圆区域。

对于地面真实亮温T_b，可以认为是其足印内n种下垫面亮温值的线性叠加：

$T_b(x,y)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)---\left(2\right)$

其中T_i代表不同下垫面(组分)亮温，L_i代表微波混合像元中不同下垫面出现的比例，且满足：

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{L}_i=1---\left(3\right)$

将(2)代入(1)式有：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)\mathrm{dxdy}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{T}_i---\left(4\right)$

其中，α是天线增益函数

和对应足印内下垫面分布L的卷积，相应的积分覆盖整个椭圆足印：

$\mathrm{\&alpha;}(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{dxdy}---\left(5\right)$

实际的天线增益方向图可以近似为高斯分布函数。因此，天线增益函数的公式定义为：

$\stackrel{\&OverBar;}{G}(x,y)=g_0\exp [-g_1(\frac{x^2}{a^2}+\frac{y^2}{b^2})]=g_0\exp [-g_1{r}^2]---\left(6\right)$

其中g₀和g₁是描述高斯分布的常数，x和y是当前位置到测量中心点的水平和垂直距离。a和b表示-3dB椭圆的半轴，可以根据被动微波数据中不同频率对应的足印来确定。椭圆半径r定义为：

$r=\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}}---\left(7\right)$

-3dB点定义为从中心点到天线增益降低到最大值一半的位置，因此有：

G(-3dB)＝0.5g₀→e^-g1＝0.5→g₁＝ln2                      (8)

定义整个天线增益为1，g₀可以通过下式计算：

$1={\&Integral;}_0^{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}{\mathrm{rg}}_0\mathrm{drd\&phi;}---\left(9\right)$

$g_0=\frac{\ln 2}{\mathrm{\&pi;}}---\left(10\right)$

位于-3dB的足印椭圆内点(x，y)满足：

$\frac{x^2}{{\left(a\right)}^2}+\frac{y^2}{{\left(b\right)}^2}\&le;1---\left(11\right)$

积分可以近似为：

$\mathrm{\&alpha;}=\underset{x=-a}{\overset{a}{\mathrm{\&Sigma;}}}\left(\underset{y=-\mathrm{bf}\left(x\right)}{\overset{\mathrm{bf}\left(x\right)}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{\&Delta;y}\right)\mathrm{\&Delta;x}---\left(12\right)$

其中：

$f\left(x\right)=\sqrt{1-\frac{x^2}{{\left(a\right)}^2}}---\left(13\right)$

式(12)中水平采样率Δx和垂直采样率Δy的选择由地物分类数据的空间分辨率确定。

(3)建立被动微波积雪混合像元模型

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

考虑到影响积雪亮温的五种主要下垫面类型，地表积雪混合像元模型应满足如下数学表达公式：

T_b(x，y)＝T_land(x，y)L_land(x，y)+T_water(x，y)L_water(x，y)+T_forest(x，y)L_forest(x，y)+T_grass(x，y)L_grass(x，y)+T_crop(x，y)L_crop(x，y)                   (14)

其中T_b代表积雪混合像元亮温值，(x，y)代表混合像元位置；T_land代表下垫面为裸土的亮温值，L_land代表与该混合像元空间位置匹配的分类数据中裸土像元的比例；T_water代表下垫面为水体(冰)的亮温值，L_water代表与该混合像元空间位置匹配的分类数据中水体(冰)像元的比例；T_forest代表下垫面为林地的亮温值，L_forest代表与该混合像元空间位置匹配的分类数据中林地像元的比例；T_grass代表下垫面为草地的亮温值，L_grass代表与该混合像元空间位置匹配的分类数据中草地像元的比例；T_crop代表下垫面为农田的亮温值，L_crop代表与该混合像元空间位置匹配的分类数据中农田像元的比例。

微波混合像元中，所有下垫面类型出现的比例应满足：

L_land(x，y)+L_water(x，y)+L_forest(x，y)+L_grass(x，y)+L_crop(x，y)＝1    (15)

卫星辐射计获得的积雪被动微波混合像元T_B可以通过天线增益函数和地表积雪混合像元模型的卷积获得，将(14)代入(1)得到积雪被动微波混合像元分解模型：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\&CenterDot;\{L_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+L_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+L_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+L_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+L_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)\}\mathrm{dxdy}---\left(16\right)$

$={\mathrm{\&alpha;}}_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+{\mathrm{\&alpha;}}_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)$

其中：

${\mathrm{\&alpha;}}_c(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L_c(x,y)\mathrm{dxdy}$ c＝land，water，forest，grass，crop  (17)

可以根据公式(12)，计算α_c的值。

其中的c＝land，water，forest，grass，crop，即，c可以分别表示裸土，水体，林地，草地，农田。

选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵。已知被动微波混合像元和地物分布比例矩阵，可以按公式(18)对m×n范围内的被动微波混合像元进行分解，通过构建方程组和约束性最小二乘法求解m×n范围内被动微波混合像元分解后的各类地物的组分亮温T。

T_B＝P_c.T_c+E                            (18)

式中：

c是地物分类的种类，对于本研究设计的积雪微波混合模型c的值为5；

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

T_c是一个(m×n×c)×1的矩阵，是m×n窗口中每个被动微波混合像元对应的不同地物的组分亮温；

P_c是一个(m×n)×(m×n×c)的矩阵，是m×n窗口中每个被动微波混合像元对应的不同地物的分布与天线方向图的卷积；

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

若m×n窗口内所有被动微波混合像元的下垫面类型完全相同，则认为m×n个被动微波混合像元为纯像元，不需要分解。由(18)式可见，该混合像元分解模型中，m×n方程组中有(m×n)×c个未知量，对于这种方程个数小于未知量个数的方程组，属于欠定方程组。欠定方程组有无穷多组解，如何求得最优解是该算法需要解决的问题。

(4)欠定性方程组的求解

研究提出在整个解空间中搜索最优解的计算方法，对每个组分亮温值设定其对应的取值区间。在该区间内，任意一点都对应着一个可能的解，这样可以构成完整的解空间。其中，各下垫面初始亮温值的选取对于能否求得最优解尤为重要。

由于不同下垫面亮温在某一段时间内基本在一定区间内变化，可以通过地基和星载微波辐射计对不同下垫面类型进行观测，得到不同下垫面的亮温变化区间。考虑冬季气候在一段时间的稳定性，首先采用k-means聚类算法对观测地区连续一段时间的被动微波混合像元亮温进行统计分类。k-means聚类后的聚类中心值可以作为不同下垫面亮温初值的参考。将地基微波数据建立的数据库中的数据、被动微波遥感数据统计得到的各类下垫面的亮温相结合，定义各类下垫面的亮温初值S_grass、S_water、S_forest、S_crop、S_land，同时定义各类下垫面的亮温变化阈值L_c，确定某一下垫面的初值亮温T_c的选取范围是[S_c-L_c，S_c+L_c]，其中c代表grass，water，land，forest，crop。

根据相关研究，不同下垫面的微波亮温变化基本满足下式：

0＜T_water＜T_land＜T_crop＜T_forest                   (19)

其中，草的亮温T_grass变化比较复杂，与草的高矮和土壤湿度有关，需要结合具体观测地区进行分析。

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

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_B(x,y)-\underset{c=1}{\overset{5}{\mathrm{\&Sigma;}}}{P}_c(x,y)\&CenterDot;T_c(x,y))}^2}<\mathrm{\&xi;}---\left(20\right)$

根据目标函数调整搜索窗口尺寸m×n，合理选择各类组分的亮温初值S_gtrass、S_water、S_forest、S_crop、S_land，代入方程(18)中，在(19)式、(20)式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解欠定性方程组。由于初值的选取范围考虑到不同下垫面的亮温变化空间，且存在多个条件限制方程组的求解范围，因此采用本发明提出的方法可以得到较好的最优解。

本发明的有益效果：

根据影响积雪辐射特性的典型五种下垫面类型，对已有的光谱遥感土地覆盖分类产品进行重新分类，利用积雪在不同下垫面微波辐射之间差异特性，建立积雪被动微波混合像元分解模型，通过混合像元分解，获得基于典型下垫面类型的更为准确的积雪被动微波亮温数据。传统方法直接使用积雪微波混合像元数据进行积雪参数(如雪深、雪水当量)的反演，混合像元的存在使积雪参数反演的结果存在很大的误差。通过采用本发明提出的方法，可以使原始像元根据其足印内不同下垫面类型分解成多个微波像元。分解后的微波像元亮温值根据其下垫面类型不同而有各自的微波辐射特征，能够及时地处理出大尺度范围下更为精确的积雪参数数据，在气候和水文研究以及积雪灾害评估等领域中具有重要应用价值。

附图说明

图1是本发明实施例1黑龙江省兴凯湖地区五种地物分类数据。

图2是本发明实施例1采样率为1km的AMSR-E天线增益函数。

图3是本发明实施例1AMSR-E 18.7GHz水平极化被动微波积雪混合像元数据。

图4是本发明实施例1AMSR-E 36.5GHz水平极化被动微波积雪混合像元数据。

图5是本发明实施例1AMSR-E 18.7GHz水平极化被动微波积雪混合像元分解数据。

图6是本发明实施例1AMSR-E 36.5GHz水平极化被动微波积雪混合像元分解数据。

图7是本发明实施例1被动微波积雪混合像元数据获得兴凯湖附近积雪覆盖区域(白色为雪)。

图8是本发明实施例1被动微波积雪混合像元分解数据获得兴凯湖附近积雪覆盖区域(白色为雪)。

图9是本发明实施例2吉林省地区五种地物分类数据。

图10是本发明实施例2采样率为5km的AMSR-E天线增益函数。

图11是本发明实施例2AMSR-E 18.7GHz水平极化被动微波积雪混合像元数据。

图12是本发明实施例2AMSR-E 36.5GHz水平极化被动微波积雪混合像元数据。

图13是本发明实施例2AMSR-E 18.7GHz水平极化被动微波积雪混合像元分解数据。

图14是本发明实施例2AMSR-E 36.5GHz水平极化被动微波积雪混合像元分解数据。

图15是本发明实施例2被动微波积雪混合像元数据获得吉林省积雪覆盖区域(白色为雪)。

图16是本发明实施例2被动微波积雪混合像元分解数据获得吉林省积雪覆盖区域(白色为雪)。

具体实施方式

实施例1：

本发明利用AMSR-E被动微波遥感数据和MODIS土地覆盖分类UMD数据产品，结合提出的积雪被动微波遥感混合像元分解方法，实现了2010年11月26日中国黑龙江省兴凯湖地区积雪被动微波混合像元的有效分解。

具体包括以下步骤：

一、地物分类数据的获取：

积雪微波混合像元分解方法，需要已知观测地区的下垫面地物分类信息。积雪微波亮温值受其下垫面影响较大，通过地基实验发现林地上积雪亮温受树干的影响较大，草地上积雪亮温主要受草的散射影响，农田上积雪亮温受农田几何形状的影响较大，冰上的积雪亮温随雪深增加亮温值增加，裸土上的积雪亮温随雪深增加亮温值降低，影响积雪亮温的其他下垫面类型基本都包含在这五种典型下垫面类型中。因此，积雪下垫面分类可以定义为水体(冰)、草地、林地、农田和裸土，对同一类型的地物不需要更加细致的划分。

在MODIS官方网站下载中分辨率MODIS光谱遥感数据的土地覆盖类型产品。MODIS陆地研究小组(MODIS Land Team)开发了MODIS Aqua和Terra卫星数据合成的年度土地覆盖分类产品MCD12Q1(空间分辨率500m)，包含IGBP(国际地圈生物圈计划)、UMD(马里兰大学修订版IGBP)、LA I/FPAR(叶面积指数/光合有效辐射分量)、NPP(净初级生产力)和PFT(植物功能类)5种土地覆盖分类数据集，信息提取主要技术是监督决策树分类。研究中选择UMD年度地物分类数据集，时间是2009年。

对下载的UMD地物分类数据集使用MRT软件进行图像拼接、等经纬度投影，采样方法为邻近法，椭球为WGS-84体系；根据兴凯湖地区的地理信息，在ArcGIS软件下加载该地区shapefile矢量文件，得到该地区的光谱遥感数据的土地覆盖类型数据，并通过ArcGIS软件将该数据转换为Grid网格，空间分辨率为0001°×0001°(1km×1km)。

UMD地物分类数据产品中土地覆盖类型多于研究所需的五种类型，需要将分类数据产品按照这五种类型重新分类。重新分类的总体原则如下：

(a)原分类数据中定义为农田、阔叶作物或谷类作物的类型统一重新定义为农田类型，因此，UMD中代码为12农田的数据重新定义为代码1农田。

(b)原分类数据中定义为阔叶林、针叶林、混合林或灌丛的类型统一重新定义为林地类型，因此，UMD中代码分别为1常绿针叶林、2常绿阔叶林、3落叶针叶林、4落叶阔叶林、5混生林、6郁闭灌丛和7稀疏灌丛的数据重新定义为代码2林地。

(c)原分类数据中定义为草原、草地、植被类型统一重新定义为草地类型，因此，UMD中代码分别为8多树的草原、9稀树草原和10草原的数据重新定义为代码3草地。

(d)原分类数据中定义为水体类型仍旧定义为水体类型，因此，UMD中代码为0水体的数据重新定义为代码4水体；

(e)原分类数据中定义为城建用地、荒漠或荒地的类型统一重新定义为裸土类型，因此，UMD中代码分别为13城建用地和16荒漠或荒地的数据重新定义为代码5裸土。

采用上述分类方法，可得到图1所示的兴凯湖地区的土地覆盖分类结果，数据主要包括水体(冰)、草地、林地、农田和裸土五类。

二、微波天线增益函数和采样率的选择

卫星上辐射计天线测得的亮温是地表实际亮温和天线增益函数的积分：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)T_b(x,y)\mathrm{dxdy}---\left(1\right)$

公式(1)中，T_B是辐射计天线测量的亮温值，T_b是在位置(x，y)处的实际亮温值，

是天线增益函数(允许仪器在测试期间扫描移动)，积分覆盖整个足印的椭圆区域。

对于地面真实亮温T_b，可以认为是其足印内五种下垫面亮温值的线性叠加：

$T_b(x,y)=\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)---\left(2\right)$

其中T_i代表不同下垫面(组分)亮温，L_i代表微波混合像元中五种不同下垫面出现的比例，且满足：

$\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i=1---\left(3\right)$

将(2)代入(1)式有：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)\mathrm{dxdy}=\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{T}_i---\left(4\right)$

其中，α是天线增益函数

和对应足印内下垫面分布L的卷积，相应的积分覆盖整个椭圆足印：

$\mathrm{\&alpha;}(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{dxdy}---\left(5\right)$

实际的天线增益方向图可以近似为高斯分布函数。因此，天线增益函数的公式定义为：

$\stackrel{\&OverBar;}{G}(x,y)=g_0\exp [-g_1(\frac{x^2}{a^2}+\frac{y^2}{b^2})]=g_0\exp [-g_1{r}^2]---\left(6\right)$

其中g₀和g₁是描述高斯分布的常数，x和y是当前位置到测量中心点的水平和垂直距离。a和b表示-3dB椭圆的半轴，可以根据被动微波数据中不同频率对应的足印来确定。椭圆半径r定义为：

$r=\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}}---\left(7\right)$

-3dB点定义为从中心点到天线增益降低到最大值一半的位置，因此有：

G(一3dB)＝0.5-g₀→e^-g1＝0.5→g₁＝ln 2              (8)

定义整个天线增益为1，g₀可以通过下式计算：

$1={\&Integral;}_0^{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}{\mathrm{rg}}_0\mathrm{drd\&phi;}---\left(9\right)$

$g_0=\frac{\ln 2}{\mathrm{\&pi;}}---\left(10\right)$

位于-3dB的足印椭圆内点(x，y)满足：

$\frac{x^2}{{\left(a\right)}^2}+\frac{y^2}{{\left(b\right)}^2}\&le;1---\left(11\right)$

积分可以近似为：

$\mathrm{\&alpha;}=\underset{x=-a}{\overset{a}{\mathrm{\&Sigma;}}}\left(\underset{y=-\mathrm{bf}\left(x\right)}{\overset{\mathrm{bf}\left(x\right)}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{\&Delta;y}\right)\mathrm{\&Delta;x}---\left(12\right)$

其中：

$f\left(x\right)=\sqrt{1-\frac{x^2}{{\left(a\right)}^2}}---\left(13\right)$

式(12)中水平采样率Δx和垂直采样率Δy的选择由地物分类数据的空间分辨率确定，研究选择采样率为1km。同时，根据AMSR-E被动微波辐射计的参数指标，(12)式中a和b选择为25km。根据所选采样率，获得图2所示的AMSR-E被动微波天线增益函数。

三、建立被动微波积雪混合像元模型

从网站上下载冬季期间AMSR-E被动微波遥感数据，该数据已经过标定、大气校正、地理校正和标准化等预处理过程；数据采集时间是2010年11月26日，为全球降轨数据。考虑到AMSR-E的综合空间分辨率为25km，利用ENVI软件将其对于积雪辐射特性明显的两个频段：图3所示的18.7GHZ水平极化数据和图4所示的36.5GHZ水平极化数据，重新投影到全球0.25°×0.25°等经纬度网格；根据黑龙江省兴凯湖地区的shapefile矢量文件，利用ArcGIS软件获取AMSR-E在兴凯湖地区的微波亮温值。

AMSR-E被动微波数据的空间分辨率为0.25°×0.25°(25km×25km)，地物分类数据的空间分辨率为0.001°×0.001°(1km×1km)。在ArcGIS软件下，实现被动微波数据和地物分类数据的配准，即每个AMSR-E数据对应25×25个地物分类数据。

考虑到影响积雪亮温的五种主要下垫面类型，地表积雪混合像元模型应满足如下数学表达公式：

T_b(x，y)＝T_land(x，y)L_land(x，y)+T_water(x，y)L_water(x，y)+T_forest(x，y)L_forest(x，y)+T_grass(x，y)L_grass(x，y)+T_crop(x，y)L_crop(x，y)                (14)

其中T_b代表积雪混合像元亮温值，(x，y)代表混合像元位置；T_land代表下垫面为裸土的亮温值，L_land代表与该混合像元空间位置匹配的分类数据中裸土像元的比例；T_water代表下垫面为水体(冰)的亮温值，L_water代表与该混合像元空间位置匹配的分类数据中水体(冰)像元的比例；T_forest代表下垫面为林地的亮温值，L_forest代表与该混合像元空间位置匹配的分类数据中林地像元的比例；T_grass代表下垫面为草地的亮温值，L_grass代表与该混合像元空间位置匹配的分类数据中草地像元的比例；T_crop代表下垫面为农田的亮温值，L_crop代表与该混合像元空间位置匹配的分类数据中农田像元的比例。

微波混合像元中，所有下垫面类型出现的比例应满足：

L_land(x，y)+L_water(x，y)+L_forest(x，y)+L_grass(x，y)+L_crop(x，y)＝1    (15)

卫星辐射计获得的积雪被动微波混合像元T_B可以通过天线增益函数和地表积雪混合像元模型的卷积获得，将(14)代入(1)得到积雪被动微波混合像元分解模型：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\&CenterDot;\{L_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+L_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+L_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+L_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+L_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)\}\mathrm{dxdy}$ (16)

$={\mathrm{\&alpha;}}_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+{\mathrm{\&alpha;}}_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)$

其中：

${\mathrm{\&alpha;}}_c(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L_c(x,y)\mathrm{dxdy}$ c＝land，water，forest，grass，crop  (17)

可以根据公式(12)，计算α_c的值。

选取2×2范围的被动微波混合像元构成一个搜索窗口，记录2×2窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵。已知被动微波混合像元和地物分布比例矩阵，可以按公式(18)对2×2范围内的被动微波混合像元进行分解，通过构建方程组和约束性最小二乘法求解2×2范围内被动微波混合像元分解后的各类地物的组分亮温T_c。

T_B＝P_c.T_c+E                       (18)

式中：

c是五种地物覆盖类型，即c＝land，water，forest，grass，crop；

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

T_c是一个(2×2×5)×1的矩阵，是2×2窗口中每个被动微波混合像元对应的五种地物的组分亮温；

P_c是一个(2×2)×(2×2×5)的矩阵，是2×2窗口中每个被动微波混合像元对应的五种地物的分布与天线方向图的卷积结果；

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

若2×2窗口内所有被动微波混合像元的下垫面类型完全相同，则认为2×2个被动微波混合像元为纯像元，不需要分解。由(18)式可见，该混合像元分解模型中，2×2方程组中有(2×2)×5个未知量，对于这种方程个数小于未知量个数的方程组，属于欠定方程组。欠定方程组有无穷多组解，如何求得最优解是该算法需要解决的问题。

四、欠定性方程组的求解

研究提出在整个解空间中搜索最优解的计算方法，对每个下垫面亮温值设定其对应的取值区间。在该区间内，任意一点都对应着一个可能的解，这样可以构成完整的解空间。其中，各下垫面初始亮温值的选取对于能否求得最优解尤为重要。

由于不同下垫面亮温在某一段时间内基本在一定区间内变化，可以通过地基和星载微波辐射计对不同下垫面类型进行观测，得到不同下垫面的亮温变化区间。考虑冬季气候在一段时间的稳定性，首先采用k-means聚类算法对观测地区连续一段时间的被动微波混合像元亮温进行统计分类。k-means聚类后的聚类中心值可以作为不同下垫面亮温初值的参考。k-means聚类后的聚类中心值代表五种下垫面亮温初值，即S_grass、S_water、S_forest、S_crop、和S_land，经计算S_grass＝205k、S_water＝160k、S_forest＝247k、S_crop＝229k、和S_land＝180k。将地基微波数据建立的数据库中的数据、被动微波遥感数据统计得到的各类下垫面的亮温相结合，考虑昼夜温差对于五种下垫面微波遥感数据的影响，定义草地的亮温变化阈值为10k，其初值亮温S_grass的选取范围是205±10k；水体的亮温变化阈值为3k，其初值亮温S_water的选取范围是160±3k；林地的亮温变化阈值为6k，其初值亮温S_forest的选取范围是247±6k；农田的亮温变化阈值为8k，其初值亮温S_crop的选取范围是229±8k；裸土的亮温变化阈值为9k，其初值亮温S_land的选取范围是180±9k。

观测地区中五种下垫面的微波亮温变化基本满足下式：

0＜T_water＜T_land＜T_grass＜T_crop＜T_forest                  (19)

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

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_B(x,y)-\underset{c=1}{\overset{5}{\mathrm{\&Sigma;}}}{P}_c(x,y)\&CenterDot;T_c(x,y))}^2}<\mathrm{\&xi;}---\left(20\right)$

根据目标函数调整搜索窗口尺寸2×2，合理选择各类组分的亮温初值S_grass、S_water、S_forest、S_crop、S_land，代入方程(18)中，在(19)式、(20)式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解欠定性方程组。由于初值的选取范围考虑到不同下垫面的亮温变化空间，且存在多个条件限制方程组的求解范围，因此采用本研究提出的方法可以得到较好的最优解。图5是积雪被动微波混合像元分解后获得18.7GHZ水平极化数据，图6是积雪被动微波混合像元分解后36.5GHZ水平极化数据。

实验结果：

为了评价积雪被动微波混合像元对于积雪参数反演的影响，本实施例使用积雪被动微波混合像元数据和积雪被动微波混合像元模型分解后的数据来实现兴凯湖附近积雪覆盖区域的判定。

根据Chang的微波雪深反演方法，积雪越厚，散射效率越高，表现在18.7GHZ和36.5GHZ两个频段上的亮度温度的差值越大，可以利用这两个频段差值大于零来判定观测地区是否有雪覆盖。

图7是利用AMSR-E两个频段的积雪被动微波混合像元数据计算得到的兴凯湖地区雪盖区域，图8是利用被动微波混合像元分解数据计算得到的兴凯湖地区雪盖区域。由图7和图8对比可见，由于被动微波数据的空间分辨率较低，原始数据中存在大量混合像元，因此利用该数据计算得到的雪盖区域存在不准确性，而通过采用提出的混合像元分解方法，考虑到下垫面地物对于积雪辐射特性的影响，对微波混合像元数据进行分解，使用分解后的数据进行雪盖判定的结果更加符合观测地区真实积雪覆盖情况。

实施例2：

本发明利用AMSR-E被动微波遥感数据和MODIS土地覆盖分类LAI/FPAR数据产品，结合提出的积雪被动微波遥感混合像元分解方法，实现了2010年11月13日中国吉林省地区积雪被动微波混合像元的有效分解。

具体包括以下步骤：

一、地物分类数据的获取：

积雪微波混合像元分解方法，需要已知观测地区的下垫面地物分类信息。积雪微波亮温值受其下垫面影响较大，通过地基实验发现林地上积雪亮温受树干的影响较大，草地上积雪亮温主要受草的散射影响，农田上积雪亮温受农田几何形状的影响较大，冰上的积雪亮温随雪深增加亮温值增加，裸土上的积雪亮温随雪深增加亮温值降低，影响积雪亮温的其他下垫面类型基本都包含在这五种典型下垫面类型中。因此，积雪下垫面分类可以定义为水体(冰)、草地、林地、农田和裸土，对同一类型的地物不需要更加细致的划分。

在MODIS官方网站下载中分辨率MODIS光谱遥感数据的土地覆盖类型产品。MODIS陆地研究小组(MODIS Land Team)开发了MODIS Aqua和Terra卫星数据合成的年度土地覆盖分类产品MCD12Q1(空间分辨率500m)，包含IGBP(国际地圈生物圈计划)、UMD(马里兰大学修订版IGBP)、LAI/FPAR(叶面积指数/光合有效辐射分量)、NPP(净初级生产力)和PFT(植物功能类)5种土地覆盖分类数据集，信息提取主要技术是监督决策树分类。研究选择LAI/FPAR年度地物分类数据集，时间是2009年。

对下载的LAI/FPAR地物分类数据集使用MRT软件进行图像拼接、等经纬度投影，采样方法为邻近法，椭球为WGS-84体系；根据吉林省的地理信息，在ArcGIS软件下加载该地区shapefile矢量文件，得到该地区的光谱遥感数据的土地覆盖类型数据，并通过ArcGIS软件将该数据转换为Grid网格，空间分辨率为0.005°×0.005°(5km×5km)。

LAI/FPAR地物分类数据产品中土地覆盖类型多于研究所需的五种类型，需要将分类数据产品按照这五种类型重新分类。重新分类的总体原则如下：

(a)原分类数据中定义为农田、阔叶作物或谷类作物的类型统一重新定义为农田类型，因此，LAI/FPAR中代码为3阔叶作物的数据重新定义为代码1农田。

(b)原分类数据中定义为阔叶林、针叶林、混合林或灌丛的类型统一重新定义为林地类型，因此，LAI/FPAR中代码分别为5阔叶林和6针叶林的数据重新定义为代码2林地。

(c)原分类数据中定义为草原、草地、植被类型统一重新定义为草地类型，因此，LAI/FPAR中代码分别为1草地/谷物作物、2灌丛和4稀树草原的数据重新定义为代码3草地。

(d)原分类数据中定义为水体类型仍旧定义为水体类型，因此，LAI/FPAR中代码为0水体的数据重新定义为代码4水体；

(e)原分类数据中定义为城建用地、荒漠或荒地的类型统一重新定义为裸土类型，因此，LAI/FPAR中代码分别为8城镇和7无植被陆地的数据重新定义为代码5裸土。

采用上述分类方法，可得到图9所示的吉林省的土地覆盖分类结果，数据主要包括水体(冰)、草地、林地、农田和裸土五类。

二、微波天线增益函数和采样率的选择

卫星上辐射计天线测得的亮温是地表实际亮温和天线增益函数的积分：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)T_b(x,y)\mathrm{dxdy}---\left(1\right)$

公式(1)中，T_B是辐射计天线测量的亮温值，T_b是在位置(x，y)处的实际亮温值，

是天线增益函数(允许仪器在测试期间扫描移动)，积分覆盖整个足印的椭圆区域。

对于地面真实亮温T_b，可以认为是其足印内五种下垫面亮温值的线性叠加：

$T_b(x,y)=\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)---\left(2\right)$

其中T_i代表不同下垫面(组分)亮温，L_i代表微波混合像元中五种不同下垫面出现的比例，且满足：

$\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i=1---\left(3\right)$

将(2)代入(1)式有：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{L}_i(x,y)T_i(x,y)\mathrm{dxdy}=\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{T}_i---\left(4\right)$

其中，α是天线增益函数

和对应足印内下垫面分布L的卷积，相应的积分覆盖整个椭圆足印：

$\mathrm{\&alpha;}(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{dxdy}---\left(5\right)$

实际的天线增益方向图可以近似为高斯分布函数。因此，天线增益函数的公式定义为：

$\stackrel{\&OverBar;}{G}(x,y)=g_0\exp [-g_1(\frac{x^2}{a^2}+\frac{y^2}{b^2})]=g_0\exp [-g_1{r}^2]---\left(6\right)$

其中g₀和g₁是描述高斯分布的常数，x和y是当前位置到测量中心点的水平和垂直距离。a和b表示-3 dB椭圆的半轴，可以根据被动微波数据中不同频率对应的足印来确定。椭圆半径r定义为：

$r=\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}}---\left(7\right)$

-3 dB点定义为从中心点到天线增益降低到最大值一半的位置，因此有：

G(-3dB)＝0.5g₀→e^-g1＝0.5→g₁＝ln2               (8)

定义整个天线增益为1，g₀可以通过下式计算：

$1={\&Integral;}_0^{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}{\mathrm{rg}}_0\mathrm{drd\&phi;}---\left(9\right)$

$g_0=\frac{\ln 2}{\mathrm{\&pi;}}---\left(10\right)$

位于-3 dB的足印椭圆内点(x，y)满足：

$\frac{x^2}{{\left(a\right)}^2}+\frac{y^2}{{\left(b\right)}^2}\&le;1---\left(11\right)$

积分可以近似为：

$\mathrm{\&alpha;}=\underset{x=-a}{\overset{a}{\mathrm{\&Sigma;}}}\left(\underset{y=-\mathrm{bf}\left(x\right)}{\overset{\mathrm{bf}\left(x\right)}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{\&Delta;y}\right)\mathrm{\&Delta;x}---\left(12\right)$

其中：

$f\left(x\right)=\sqrt{1-\frac{x^2}{{\left(a\right)}^2}}---\left(13\right)$

式(12)中水平采样率Δx和垂直采样率Δy的选择由地物分类数据的空间分辨率确定，研究选择采样率为5 km。同时，根据AMSR-E被动微波辐射计的参数指标，(12)式中a和b选择为25 km。根据所选采样率，获得图10所示的AMSR-E被动微波天线增益函数方向图。

三、建立被动微波积雪混合像元模型

从网站上下载冬季期间AMSR-E被动微波遥感数据，该数据已经过标定、大气校正、地理校正和标准化等预处理过程；数据采集时间是2010年11月13日，为全球降轨数据。考虑到AMSR-E的综合空间分辨率为25 km，利用ENVI软件将其对于积雪辐射特性明显的两个频段：图11所示的18.7GHZ水平极化数据和图12所示的36.5GHZ水平极化数据，重新投影到全球0.25°×0.25°等经纬度网格；根据吉林省的shapefile矢量文件，利用ArcGIS软件获取AMSR-E在吉林省地区的微波亮温值。

AMSR-E被动微波数据的空间分辨率为0.25°×0.25°(25km×25km)，地物分类数据的空间分辨率为0.05°×0.05°(5km×5km)。在ArcGIS软件下，实现被动微波数据和地物分类数据的配准，即每个AMSR-E数据对应5×5个地物分类数据。

考虑到影响积雪亮温的五种主要下垫面类型，地表积雪混合像元模型应满足如下数学表达公式：

T_b(x，y)＝T_land(x，y)L_land(x，y)+T_water(x，y)L_water(x，y)+T_forest(x，y)L_forest(x，y)+T_grass(x，y)L_grass(x，y)+T_crop(x，y)L_crop(x，y)                  (14)

其中T_b代表积雪混合像元亮温值，(x，y)代表混合像元位置；T_land代表下垫面为裸土的亮温值，L_land代表与该混合像元空间位置匹配的分类数据中裸土像元的比例；T_water代表下垫面为水体(冰)的亮温值，L_water代表与该混合像元空间位置匹配的分类数据中水体(冰)像元的比例；T_forest代表下垫面为林地的亮温值，L_forest代表与该混合像元空间位置匹配的分类数据中林地像元的比例；T_grass代表下垫面为草地的亮温值，L_grass代表与该混合像元空间位置匹配的分类数据中草地像元的比例；T_crop代表下垫面为农田的亮温值，L_crop代表与该混合像元空间位置匹配的分类数据中农田像元的比例。

微波混合像元中，所有下垫面类型出现的比例应满足：

L_land(x，y)+L_water(x，y)+L_forest(x，y)+L_grass(x，y)+L_crop(x，y)＝1    (15)

卫星辐射计获得的积雪被动微波混合像元T_B可以通过天线增益函数和地表积雪混合像元模型的卷积获得，将(14)代入(1)得到积雪被动微波混合像元分解模型：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\&CenterDot;\{L_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+L_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+L_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+L_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+L_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)\}\mathrm{dxdy}---\left(16\right)$

$={\mathrm{\&alpha;}}_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+{\mathrm{\&alpha;}}_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)$

其中：

${\mathrm{\&alpha;}}_c(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L_c(x,y)\mathrm{dxdy}$ c＝land，water，forest，grass，crop  (17)

可以根据公式(12)，计算α_c的值。

选取4×4范围的被动微波混合像元构成一个搜索窗口，记录4×4窗口中每个被动微波像元中不同地物的出现比例，构成地物分布比例矩阵。已知被动微波混合像元和地物分布比例矩阵，可以按公式(18)对4×4范围内的被动微波混合像元进行分解，通过构建方程组和约束性最小二乘法求解4×4范围内被动微波混合像元分解后的各类地物的组分亮温T_c。

T_B＝P_c.T_c+E                      (18)

式中：

c是五种地物覆盖类型，即c＝land，water，forest，grass，crop；

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

T_c是一个(4×4×5)×1的矩阵，是4×4窗口中每个被动微波混合像元对应的

五种地物的组分亮温；

P_c是一个(4×4)×(4×4×5)的矩阵，是4×4窗口中每个被动微波混合像元对

应的五种地物的分布与天线方向图的卷积结果；

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

若4×4窗口内所有被动微波混合像元的下垫面类型完全相同，则认为4×4个被动微波混合像元为纯像元，不需要分解。由(18)式可见，该混合像元分解模型中，4×4方程组中有(4×4)×5个未知量，对于这种方程个数小于未知量个数的方程组，属于欠定方程组。欠定方程组有无穷多组解，如何求得最优解是该算法需要解决的问题。

四、欠定性方程组的求解

研究提出在整个解空间中搜索最优解的计算方法，对每个下垫面亮温值设定其对应的取值区间。在该区间内，任意一点都对应着一个可能的解，这样可以构成完整的解空间。其中，各下垫面初始亮温值的选取对于能否求得最优解尤为重要。

由于不同下垫面亮温在某一段时间内基本在一定区间内变化，可以通过地基和星载微波辐射计对不同下垫面类型进行观测，得到不同下垫面的亮温变化区间。考虑冬季气候在一段时间的稳定性，首先采用k-means聚类算法对观测地区连续一段时间的被动微波混合像元亮温进行统计分类。k-means聚类后的聚类中心值可以作为不同下垫面亮温初值的参考。k-means聚类后的聚类中心值代表五种下垫面亮温初值，即S_grass、S_water、S_forest、S_crop、和S_land，经计算S_grass＝229k、S_water＝180k、S_forest＝270k、S_crop＝252k、和S_land＝210k。将地基微波数据建立的数据库中的数据、被动微波遥感数据统计得到的各类下垫面的亮温相结合，考虑昼夜温差对于五种下垫面微波遥感数据的影响，定义草地的亮温变化阈值为8k，其初值亮温S_grass的选取范围是229±8k；水体的亮温变化阈值为5k，其初值亮温S_water的选取范围是180±5k；林地的亮温变化阈值为6k，其初值亮温S_forest的选取范围是270±6k；农田的亮温变化阈值为10k，其初值亮温S_crop的选取范围是252±10k；裸土的亮温变化阈值为8k，其初值亮温S_land的选取范围是210±8k。

观测地区中五种下垫面的微波亮温变化基本满足下式：

0＜T_water＜T_land＜T_grass＜T_crop＜T_forest                  (19)

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

$R=\sqrt{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_B(x,y)-\underset{c=1}{\overset{5}{\mathrm{\&Sigma;}}}{P}_c(x,y)\&CenterDot;T_c(x,y))}^2}<\mathrm{\&xi;}---\left(20\right)$

根据目标函数调整搜索窗口尺寸4×4，合理选择各类组分的亮温初值S_grass、S_water、S_forest、S_crop、S_land，代入方程(18)中，在(19)式、(20)式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解欠定性方程组。由于初值的选取范围考虑到不同下垫面的亮温变化空间，且存在多个条件限制方程组的求解范围，因此采用本研究提出的方法可以得到较好的最优解。积雪被动微波混合像元分解后，获得图13所示的吉林省地区被动微波18.7GHZ水平极化数据和图14所示的36.5GHZ水平极化数据。

实验结果：

为了评价积雪被动微波混合像元对于积雪参数反演的影响，本实施例使用积雪被动微波混合像元数据和积雪被动微波混合像元模型分解后的数据来实现吉林省地区积雪覆盖区域的判定。

根据Chang的微波雪深反演方法，积雪辐射特性在18.7GHZ和36.5GHZ两个频段上比较显著，且随频段的增加辐射亮温值降低，因此可以利用这两个频段差值大于零来判定观测地区是否有雪覆盖。

图15是利用AMSR-E两个频段的积雪被动微波混合像元数据计算得到的吉林省地区雪盖区域，图16是利用被动微波混合像元分解数据计算得到的吉林省地区雪盖区域。由图15和图16对比可见，由于被动微波数据的空间分辨率较低，原始数据中存在大量混合像元，因此利用该数据计算得到的雪盖区域存在不准确性和分辨率较低的特点，采用提出的混合像元分解方法，可以利用下垫面地物的辐射特性对积雪混合像元数据进行分解，使用分解后的数据进行的雪盖判定结果更加符合吉林省地区真实积雪覆盖情况。

Claims (3)

1.一种基于五种地物分类信息的积雪被动微波混合像元分解方法，该方法的应用条件是冬季积雪被动微波遥感数据，其特征在于，方法包括如下过程：1)确定观测地区冬季积雪情况下的五种的地物分类数据，2)确定微波天线增益函数和采样率，3)建立积雪被动微波混合像元分解模型，4)积雪被动微波混合像元分解模型求解；

所述的确定观测地区冬季积雪情况下的五种的地物分类数据，是将已有的观测地区的地物分类数据产品重新分类为水体、草地、林地、农田和裸土五种下垫面信息，确定观测地区冬季积雪情况下的五种地物分类数据的结果L；

所述的确定微波天线增益函数和采样率，是根据被动微波辐射计的特性，选择被动微波数据中不同频率对应的椭圆足印a和b，设计微波天线增益函数

水平采样率Δx和垂直采样率Δy，结合地物分类数据的结果L，利用公式（12）求得系数α：

$\mathrm{\&alpha;}=\underset{x=-a}{\overset{a}{\mathrm{\&Sigma;}}}\left(\underset{y=-\mathrm{bf}\left(x\right)}{\overset{\mathrm{bf}\left(x\right)}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{G}(x,y)L(x,y)\mathrm{\&Delta;y}\right)\mathrm{\&Delta;x}---\left(12\right)$

其中 $f\left(x\right)=\sqrt{1-\frac{x^2}{{\left(a\right)}^2}};$

α是天线增益函数

和对应足印内地物分类数据的结果L的卷积，(x,y)是被动微波数据的位置；

所述的建立积雪被动微波混合像元分解模型，是将地物分类数据和微波天线增益函数及采样率的计算结果代入积雪被动微波混合像元分解模型，确定积雪被动微波混合像元分解模型求解方法；具体的积雪被动微波混合像元分解模型为：

$T_B=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)\&CenterDot;\{L_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+L_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+L_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+L_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+L_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)\}\mathrm{dxdy}---\left(\text{16}\right)$

$={\mathrm{\&alpha;}}_{\mathrm{land}}(x,y)T_{\mathrm{land}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{water}}(x,y)T_{\mathrm{water}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{forest}}(x,y)T_{\mathrm{forest}}(x,y)$

$+{\mathrm{\&alpha;}}_{\mathrm{grass}}(x,y)T_{\mathrm{grass}}(x,y)+{\mathrm{\&alpha;}}_{\mathrm{crop}}(x,y)T_{\mathrm{crop}}(x,y)$

其中：

${\mathrm{\&alpha;}}_c(x,y)=\&Integral;\&Integral;\stackrel{\&OverBar;}{G}(x,y)L_c(x,y)\mathrm{dxdy},c=\mathrm{land},\mathrm{water},\mathrm{forest},\mathrm{grass},\mathrm{crop}---\left(17\right)$

公式（16）中，T_B代表积雪混合像元亮温值，(x,y)代表混合像元空间位置，即，被动微波数据的位置；T_land、T_water、T_forest、T_grass和T_crop分别代表下垫面为裸土、水体、林地、草地和农田的亮温值，L_land、L_water、L_forest、L_grass和L_crop分别代表与混合像元空间位置(x,y)匹配的分类数据中裸土、水体、林地、草地和农田像元的比例；公式（17）的结果由公式（12）计算得到；

所述的积雪被动微波混合像元分解模型求解，是最终将积雪被动微波混合像元数据分解成五种类型的被动微波组分亮温数据；具体过程是利用建立的积雪被动微波混合像元分解模型，结合被动微波遥感数据，建立联立方程组，对m×n范围内的被动微波混合像元进行分解，通过构建方程组和约束性最小二乘法求解m×n范围内被动微波混合像元分解后的五种地物的组分亮温T_c。

2.根据权利要求1所述的一种基于五种地物分类信息的积雪被动微波混合像元分解方法，其特征在于，所述的地物分类数据产品重新分类，总体原则是：

（a）原分类数据中定义为农田、阔叶作物或谷类作物的类型统一重新定义为农田类型；

（b）原分类数据中定义为阔叶林、针叶林、混合林或灌丛的类型统一重新定义为林地类型；

（c）原分类数据中定义为草原、草地、植被类型统一重新定义为草地类型；

（d）原分类数据中定义为水体类型仍旧定义为水体类型；

（e）原分类数据中定义为城建用地、荒漠或荒地的类型统一重新定义为裸土类型。

3.根据权利要求1或2所述的一种基于五种地物分类信息的积雪被动微波混合像元分解方法，其特征在于，所述的积雪被动微波混合像元分解模型求解方法，

首先选取m×n范围的被动微波混合像元构成一个搜索窗口，记录m×n窗口中每个被动微波像元中不同地物分类的出现比例，构成地物分类分布比例矩阵，

T_BN=P_CT_C+E      （18）

式中：c是地物分类的种类，c的值为5；T_BN是一个(m×n)×1的矢量，是m×n个被动微波混合像元亮温值；T_c是一个(m×n×c)×1的矩阵，是m×n窗口中每个被动微波混合像元对应的五种地物分类的组分亮温；P_c是一个(m×n)×(m×n×c)的矩阵，是m×n窗口中每个被动微波混合像元对应的五种地物分类的分布与天线方向图的卷积；E是一个(m×n)×1的矢量，是m×n个残差数据；m选取1～10，n选取1～10；

其次求解m×n范围内被动微波混合像元分解后的五种地物的组分亮温T_c；

m×n方程组中有(m×n)×c个未知量，该方程组属于欠定方程组，有无穷多组解；对每个下垫面亮温值设定其对应的取值区间，在该区间内，任意一点都对应着一个可能的解，从而构成完整的解空间；

通过地基和星载微波辐射计对不同下垫面类型进行观测，得到不同下垫面的亮温变化区间，采用k-means聚类算法对观测地区连续一段时间的被动微波混合像元亮温进行统计分类；k-means聚类后的聚类中心值作为不同下垫面亮温初值的参考；将地基微波数据建立的数据库中的数据、被动微波遥感数据统计得到的各类下垫面的亮温相结合，定义各类下垫面的亮温初值S_grass、S_water、S_forest、S_crop、S_land，同时定义各类下垫面的亮温变化阈值L_c，确定某一下垫面的初值亮温T_c的选取范围是[S_c-L_c,S_c+L_c]，其中c代表grass,water,land,forest,crop；

不同下垫面的微波亮温变化满足下式：

0<T_water<T_land<T_crop<T_forest      （19）

草地亮温T_grass变化与草的高矮和土壤湿度有关，需要结合具体观测地区进行分析；

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

$R=\sqrt{{\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}(T_B(x,y)-\underset{c=1}{\overset{5}{\mathrm{\&Sigma;}}}{P}_c(x,y)\&CenterDot;T_c(x,y))}^2}<\mathrm{\&xi;}---\left(20\right)$

根据目标函数调整搜索窗口尺寸m×n，合理选择各类组分的亮温初值S_grass、S_water、S_forest、S_crop、S_land，代入方程（18）中，在（19）式、（20）式共同约束下，采用带约束条件的非负最小二乘法迭代运算求得最优解，应用matlab软件中的fsolve函数求解欠定性方程组。

Images