多频双极化林地积雪被动微波混合像元分解方法

摘要

本发明公开了一种多频双极化林地积雪被动微波混合像元分解方法，属于遥感图像处理的技术领域，针对现有技术利用被动微波数据对林地积雪进行雪深反演时误差较大的问题。本发明首先对林地下垫面土地分类数据的获取和重新分类，然后根据分类结果基础上建立多频双极化林地积雪被动微波混合像元分解模型，最后基于动态窗口数据选择策略的欠定性方程组求解，得到分解后各下垫面对应的组分亮温数据和误差数据，过程中本发明考虑到微波像元的空间相关性，提出了8邻域窗口数据输入和4邻域窗口数据输入两种输入数据的方案，并根据以上两种输入方案并对求解结果的优选提出了四种解决方案。

说明书

技术领域

本发明属于遥感图像处理的技术领域，具体涉及一种基于多频双极化林地积雪被动微波混合像元分解方法。

背景技术

积雪是地表最活跃的自然因素之一，地球上陆地有3/4的淡水资源以冰雪形式存在，地表表面1/3的地区存在季节性积雪，欧亚大陆和北美洲地区在冬季至少有80％的地表被积雪所覆盖。积雪是决定辐射平衡的重要因素，不仅具有强烈的气候效应，而且对能量及水分循环过程具有极其重要的影响。由于被动微波遥感具有很高的时间分辨率，能够迅速覆盖全球，因此，它在监测全球和大陆尺度的积雪时空变化中，作用尤为突出。它不仅能够全天候地观测积雪，也能够穿透大部分积雪层从而探测到雪深和雪水当量的信息。其中，积雪深度是反映积雪总量的主要因子之一，也是雪灾，融雪性洪水预警、监测和评估的重要因子。

国内外积雪研究表明下垫面性质对雪深、雪水当量的反演结果影响很大，需要对影响积雪被动微波亮温值的主要下垫面类型进行分析，在此基础上进行被动微波混合像元分解。Foster等(参见J.L.Foster,D.K.Hall,andA.T.C.Chang,(1984).“Anoverviewofpassivemicrowavesnowresearchandresults,”Rev.Geophys.SpacePhys.22(2),195-208)引入森林覆盖度参数,一定程度上提高了森林地区雪水当量的反演精度。Tait(参见A.B.Tait(1998).“EstimationofsnowwaterequivalentusingpassivemicrowaveradiationData,”RemoteSensingofEnvironment,64(3):286-291)利用SSM/I传感器的被动微波亮温数据,结合美国和俄罗斯地面雪深观测数据,根据下垫面性质(是否有湿雪、是否有深霜层的发育、是否有复杂的山区地形和是否有森林存在),将下垫面分成16类。在每一类中,利用19GHz与36GHz、19GHz与85GHz两个频率差和多种极化方式组合进行回归反演,其结果表明当下垫面无湿雪和深霜层的存在,以及无森林覆盖，平坦地区反演误差最小,相关系数最高达0.75；而在湿雪、深霜层发育的无森林覆盖区复杂地形下垫面反演效果最差,相关系数仅0.22,这也说明下垫面性质对雪水当量的反演结果影响很大。Derkesen等(参见C.Derksen,A.Walker,andB.Goodison(2005).“Evaluationofpassivemicrowavesnowwaterequivalentretrievalsacrosstheborealforest/tundraofwesternCanada,”RemoteSensingofEnvironment,96(3-4),315-327)研究加拿大西部森林区的积雪量时,发展了对不同地表覆盖类型敏感的雪深反演算法,主要考虑裸地、针叶林、落叶林和稀疏森林等四种不同下垫面,其权重即为各类型在像元内的覆盖度。针对反演效果不佳的加拿大苔原区,Derkesen(参见C.Derksen,andA.Walker(2003).“CombiningSMMRandSSM/IdatafortimeseriesanalysisofcentralNorthAmericansnowwaterJournalofHydrometeorology,4(2),304–316)采用2002/03-2006/07的冬季数据重新建立了只依靠37V极化亮温的雪深反演方法。近年来,国内学者对适用于国内的积雪雪深反演算法也进行了探索和发展。Che等人(参见T.Che,X.Li,R.Jin,R.Armstrong,andT.J.Zhang(2008).“Snowdepthderivedfrompassivemicrowaveremote-sensingdatainChina,”AnnalsofGlaciology,49(1),145-154)以Chang算法为基础,通过气象站观测数据,分别修正了中国区域的SMMR(1980年和1981年气象站点观测雪深数据)、SMM/I(2003年气象站雪深观测数据)系数,得到SMMR和SMM/I的雪深反演回归算法,反演雪深的均方根误差分别为6.22cm和5.99cm。Jiang等(参见L.M.Jiang,P.Wang,L.X.Zhang,H.Yang,andJ.Y.Yang(2014).“ImprovementofsnowdepthretrievalforFY3B-MWRIinChina,”ScienceChinaEarthSciences,44(3),531-547)利用2002-2009年的全国气象站点的地面积雪参数观测资料和相应时间、空间的AMSREL2A数据，结合中国区域的下垫面特征，把中国区域划分成森林、农田、草地和裸地四种下垫面类型，先反演这四种不同下垫面类型的雪深，然后结合权重计算混合像元内的雪深，应用到FY3B/MWRI业务化的雪深产品。Gu等(参见L.J.Gu,R.Z.Ren,K.ZhaoandX.F.Li(2014).“SnowdepthandsnowcoverretrievalfromFengYun3BmicrowaveradiationimagerybasedonasnowpassivemicrowaveunmixingmethodinNortheastChina,”J.Appl.RemoteSens.8(1),084682)基于1km空间分辨率的中国土地利用数据，结合被动微波天线增益函数，提出积雪被动微波混合像元分解模型，将中国东北地区的主要下垫面类型分为草地、农田、裸地、森林和水体五种下垫面类型，将其应用于风云三号B星搭载的微波成像仪数据。由于星载被动微波数据的空间分辨率较低(≥10km)，混合像元普遍存在于被动微波遥感数据中,它的存在是雪深后期反演精度难以达到使用要求的主要原因。在绝大部分观测区域，被动微波天线足印内均由不同的地表介质组成，微波辐射计获取的下垫面亮温是代表了一定尺度的混合像元综合亮温，而不同类型下垫面的下垫面微波辐射特性具有较大的差异，其中森林是像元内的异质性中影响最大的。

发明内容

为了解决林地像元内的异质性影响较大，而导致其在使用雪深反演方法时误差较大的问题，本发明根据林地被动微波观测像元内下垫面数据，结合被动微波天线增益函数，提供了一种多频双极化林地积雪被动微波混合像元分解方法。

本发明所采用的技术方案具体如下：

一种多频双极化林地积雪被动微波混合像元分解方法，该方法的应用条件是冬季林地观测地区的被动微波遥感数据和分辨率光谱遥感数据中的下垫面土地分类数据，方法包括如下过程：1)林地下垫面土地分类数据的获取和重新分类，2)建立多频双极化林地积雪被动微波混合像元分解模型，3)基于动态窗口数据选择策略的欠定性方程组求解。

1)林地下垫面土地分类数据的获取和重新分类

获取分辨率光谱遥感数据中的土地覆盖分类产品，将分辨率光谱遥感数据中的下垫面土地分类数据进行重新分类，具体如下：

a)将原分类数据中定义分类属于针叶植被类型的数据统一重新定义为针叶林类型B；

b)将原分类数据中定义为常绿阔叶植被、落叶阔叶植被和一年生阔叶植被类型的数据统一重新定义为阔叶林类型N；

c)其他分类数据S₁,S₂,…S_n保持原类型不变；

2)建立多频双极化林地积雪被动微波混合像元分解模型

从网站下载多频双极化被动微波遥感数据，该数据最好已经过标定、大气校正、地理校正和标准化预处理。

首先，被动微波遥感数据即卫星上辐射计天线测得地面的混合亮温值，而卫星上辐射计天线测得地面的混合亮温值是地表实际亮温和天线增益函数的积分：

$>T_F^M(x_0,y_0)=\int \int \bar{G}(x,y)T_F^C(x,y)dxdy---\left(1\right)$>

公式(1)中，是辐射计天线测量的混合亮温值，是在地面位置(x,y)实际亮温值，下脚标F是被动微波观测频段，是天线增益函数，积分覆盖整个足印的椭圆区域，天线增益函数的公式定义为：

$>\bar{G}(x,y)=\frac{ln2}{\pi }·\exp [-ln2·(\frac{x^2}{a^2}+\frac{y^2}{b^2})]---\left(2\right)$>

a和b表示被动微波数据中不同频段对应的足印。

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

$>T_F^C(x,y)=\underset{i=1}{\overset{n}{\Sigma }}{L}^i(x,y)T_i---\left(3\right)$>

其中T_i代表不同下垫面组分亮温值，Lⁱ代表第i个下垫面分类数据。

本发明中所使用的星载被动微波的频段为K频(18.7GHz)和Ka频(36.5GHz)，K频和Ka频是常用判识积雪的有效频段，每个频段均包含水平极化和垂直极化两种方式。

本发明针对林地的积雪厚度测量，因此只考虑森林主要的下垫面组成，如：阔叶林、针叶林和可能包含的其他下垫面组成，则地面位置(x,y)处的林地积雪被动微波在K频和Ka频的辐射值可以表示为：

$>\left(\begin{array}{c}{T}_{KH}^C(x,y)\\ T_{KV}^C(x,y)\\ T_{KaH}^C(x,y)\\ T_{KaV}^C(x,y)\end{array}\right)=L^B(x,y)\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+L^N(x,y)\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+L^{S_1}(x,y)\left(\begin{array}{c}{T}_{KH}^{S_1}\\ T_{KV}^{S_1}\\ T_{KaH}^{S_1}\\ T_{KaV}^{S_1}\end{array}\right)+L^{S_2}(x,y)\left(\begin{array}{c}{T}_{KH}^{S_2}\\ T_{KV}^{S_2}\\ T_{KaH}^{S_2}\\ T_{KaV}^{S_2}\end{array}\right)+...+L^{S_n}(x,y)\left(\begin{array}{c}{T}_{KH}^{S_n}\\ T_{KV}^{S_n}\\ T_{KaH}^{S_n}\\ T_{KaV}^{S_n}\end{array}\right)+ϵ---\left(4\right)$>

式(4)中T代表亮温值；L代表下垫面分类数据，若(x,y)位置存在该下垫面，则L(x,y)＝1；否则，L(x,y)＝0。上脚标C、B和S₁,S₂,…S_n分别代表该位置的下垫面土地类型、阔叶林下垫面土地类型、针叶林下垫面土地类型和其他可能存在的下垫面土地类型；下脚标K和Ka分别代表被动微波K和Ka频段，其中下脚标H和V分别代表水平极化和垂直极化；ε代表残差量；L^B(x,y)、L^N(x,y)和分别代表地面位置(x,y)处阔叶林、针叶林和其它下垫面的分布，且有：

$>L^B(x,y)+L^N(x,y)+L^{S_1}(x,y)+...+L^{S_n}(x,y)=1---\left(5\right)$>

星载被动微波辐射计在不同频段(F)和极化方式(V、H)观测到的混合亮温值亮温值可以看成是被动微波天线和地面实际下垫面类型亮温的卷积，将式(4)代入式(1)中则有：

$>\left(\begin{array}{c}{T}_{KH}^M(x_0,y_0)\\ T_{KV}^M(x_0,y_0)\\ T_{KaH}^M(x_0,y_0)\\ T_{KaV}^M(x_0,y_0)\end{array}\right)=\int \int \bar{G}(x,y)\left(L^B\right(x,y)\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+L^N(x,y)\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+L^{S_1}(x,y)\left(\begin{array}{c}{T}_{KH}^{S_1}\\ T_{KV}^{S_1}\\ T_{KaH}^{S_1}\\ T_{KaV}^{S_1}\end{array}\right)+...+L^{S_n}(x,y\left)\left(\begin{array}{c}{T}_{KH}^{S_n}\\ T_{KV}^{S_n}\\ T_{KaH}^{S_n}\\ T_{KaV}^{S_n}\end{array}\right)\right)dxdy+ϵ---\left(6\right)$>

令α是天线增益函数和下垫面分类数据L的卷积，则有：

$>{\alpha }^C(x_0,y_0)=\int \int \bar{G}(x,y)L^C(x,y)dxdy,$>C＝B,N,S(7)

利用式(7)将式(6)变形为：

$>\left(\begin{array}{c}{T}_{KH}^M(x_0,y_0)\\ T_{KV}^M(x_0,y_0)\\ T_{KaH}^M(x_0,y_0)\\ T_{KaV}^M(x_0,y_0)\end{array}\right)={\alpha }^B(x_0,y_0)\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+{\alpha }^N(x_0,y_0)\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+{\alpha }^{S_1}(x_0,y_0)\left(\begin{array}{c}{T}_{KH}^{S_1}\\ T_{KV}^{S_1}\\ T_{KaH}^{S_1}\\ T_{KaV}^{S_1}\end{array}\right)+{\alpha }^{S_2}(x_0,y_0)\left(\begin{array}{c}{T}_{KH}^{S_2}\\ T_{KV}^{S_2}\\ T_{KaH}^{S_2}\\ T_{KaV}^{S_2}\end{array}\right)+...+{\alpha }^{S_n}(x_0,y_0)\left(\begin{array}{c}{T}_{KH}^{S_n}\\ T_{KV}^{S_n}\\ T_{KaH}^{S_n}\\ T_{KaV}^{S_n}\end{array}\right)+ϵ---\left(8\right)$>

式(8)即多频双极化林地积雪被动微波混合像元分解模型。

3)基于动态窗口数据选择策略的欠定性方程组求解

根据式(8)中计算得到α值，属于线性欠定行方程组求解问题，理论上式(8)应该有无数个解。因此，应采用以下多点观测约束求解方法，如观测点个数为p时，式(8)为：

$>\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right)=\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\end{array}\right]+\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^N& ...& {\alpha }_p^N\end{array}\right]+\left(\begin{array}{c}{T}_{KH}^{S_1}\\ T_{KV}^{S_1}\\ T_{KaH}^{S_1}\\ T_{KaV}^{S_1}\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^{S_1}& ...& {\alpha }_p^{S_1}\end{array}\right]+...+\left(\begin{array}{c}{T}_{KH}^{S_n}\\ T_{KV}^{S_n}\\ T_{KaH}^{S_n}\\ T_{KaV}^{S_n}\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^{S_n}& ...& {\alpha }_p^{S_n}\end{array}\right]+E---\left(9\right)$>

式(9)中E代表残差矩阵。

通常情况下由于微波像元的空间相关性，观测点周围内所有的同类下垫面有相同的亮值，因此这里可以假设在观测点周围内所有的同类下垫面有相同的亮值，该条件可以作为公式(8)的约束，将式(8)写成：

$>\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right)=\left(\begin{array}{ccccc}{T}_{KH}^B& T_{KH}^N& T_{KH}^{S_1}& ...& T_{KH}^{S_n}\\ T_{KV}^B& T_{KV}^N& T_{KV}^{S_1}& ...& T_{KV}^{S_n}\\ T_{KaH}^B& T_{KaH}^N& T_{KaH}^{S_1}& ...& T_{KaH}^{S_n}\\ T_{KaV}^B& T_{KaV}^N& T_{KaV}^{S_1}& ...& T_{KaV}^{S_n}\end{array}\right)\left(\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\\ {\alpha }_1^N& ...& {\alpha }_p^N\\ {\alpha }_1^S& ...& {\alpha }_p^S\\ .& & .\\ .& & .\\ .& & .\\ {\alpha }_1^{S_n}& ...& {\alpha }_p^{S_n}\end{array}\right)+E---\left(10\right)$>

在这种情况下通常考虑通过最小二乘法求解：

$>T^{*}=\underset{T>0}{\mathrm{argmin}}\left|\right|Y-TA\left|\right|---\left(11\right)$>

其中，

$>Y=\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right),---\left(12\right)$>

$>T=\left(\begin{array}{ccccc}{T}_{KH}^B& T_{KH}^N& T_{KH}^{S_1}& ...& T_{KH}^{S_n}\\ T_{KV}^B& T_{KV}^N& T_{KV}^{S_1}& ...& T_{KV}^{S_n}\\ T_{KaH}^B& T_{KaH}^N& T_{KaH}^{S_1}& ...& T_{KaH}^{S_n}\\ T_{KaV}^B& T_{KaV}^N& T_{KaV}^{S_1}& ...& T_{KaV}^{S_n}\end{array}\right)---\left(13\right)$>

和

$>A=\left(\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\\ {\alpha }_1^N& ...& {\alpha }_p^N\\ {\alpha }_1^{S_1}& ...& {\alpha }_p^{S_1}\\ .& & .\\ .& & .\\ .& & .\\ {\alpha }_1^{S_n}& ...& {\alpha }_p^{S_n}\end{array}\right).---\left(14\right)$>

此外，对于某频段同一位置求解得到的被动微波组分亮温，还应该满足以下条件：

$>T_{KV}^C>T_{KH}^C---\left(15\right)$>

$>T_{KaV}^C>T_{KaH}^C---\left(16\right)$>

式(11)可以通过非负最小二乘法求解。对于解混问题，将其转换为普通的最小二乘法其满足亮温解正值的条件。

T＝(A^tA)^-1A^tY(17)

分解误差为：

E＝T-(A^tA)^-1A^tY(18)

对于公式(10)中p值选取，采用动态窗口数据选择策略：

首先，考虑到微波像元的空间相关性，以待分解的位置(x₀,y₀)处的被测的F频段的被动微波混合像元为中心,选取其周围8邻域窗口的被动微波混合像元和周围4邻域窗口的被动微波混合像元作为方程组输入数据进行求解。

对应两种输入数据的方案为：

(A)8邻域窗口数据输入：按照图1,将含当前像元在内的9个像元的多频双极化被动微波数据(即混合亮温值)转变为一维矢量代入式(12)。如果计算得到的结果满足公式(15)和公式(16),获得合理求解，得到分解后各下垫面对应的组分亮温数据T₁，得到误差数据E₁；

(B)4邻域窗口数据输入：按照图2,将含当前像元在内的5个多频双极化被动微波数据(即混合亮温值)转变为一维矢量代入公式(12)。如果计算得到的结果满足公式(15)和公式(16),获得合理求解，得到分解后各下垫面对应的组分亮温数据T₂，得到误差数据E₂。

但是根据上述方法所的组分亮温数据并不一定是最优解，因此根据(A)、(B)方案的求解结果，考虑可能出现的以下四种情况，决定最后输出的组分亮温：

情况一：选择(A)方案，若求解不合理，则采用(B)方案；若(B)方案求解仍不合理，说明输入的混合像元不能分解，保留原像元的亮温值；

情况二：选择(A)方案，若求解不合理，则采用(B)方案；若(B)方案求解合理，输出组分亮温值T₂和误差E₂，最终组分亮温值是T₂、误差数据是E₂；

情况三：选择(A)方案，若求解合理，输出组分亮温值T₁和误差矩阵E₁；继续采用(B)方案，若(B)方案求解不合理，则最终组分亮温值是T₁、误差是E₁；

情况四：选择(A)方案，若求解合理，输出组分亮温值T₁和误差矩阵E₁；继续采(B)方案，若(B)方案求解合理，输出组分亮温值T₂和误差矩阵E₂。分别对误差数据E₁和E₂的矩阵求和结果记为E₁′和E₂′，比较E₁′和E₂′的大小。若E₁′＞E₂′，则最终组分亮温值是T₂、误差数据是E₂；若E₁′＜E₂′，则最终组分亮温值是T₁、误差数据是E₁。

本发明的有益效果：

在森林积雪遥感方面，微波遥感能够发挥其穿透林木冠层的能力，有效地探测到森林冠层以下的积雪信息，这一点是光学遥感所无法企及的。利用被动微波辐射计的K频(18.7GHz)和Ka频(36.5GHz)微波辐射亮温数据进行差值计算，能够有效地获取林下积雪信息，确认了微波遥感对于林下积雪探测的有效性。但是，由于被动微波数据空间分辨率较低，微波像元内不均质性都会使得雪深的反演复杂化，特别是植被的存在会降低微波对积雪的敏感性。本发明重点针对冬季观测地区为林地类型的区域进行雪深反演研究，提出多频双极化林地积雪被动微波混合像元分解方法，按照微波像元内所包含的主要类型分为针叶林、阔叶林和其它下垫面，分别获得各下垫面上的被动微波组分亮温，将其用于针叶林、阔叶林下垫面的雪深反演，能够及时地处理出全球尺度下更为精确的不同林地类型的积雪参数数据，解决目前林地积雪参数反演精度较低的问题。

附图说明

图1是本发明8邻域窗口输入数据位置选择。

图2是本发明4邻域窗口输入数据位置选择。

图3是本发明实施例1林地观测地区被动微波辐射计K频水平极化数据。

图4是本发明实施例1林地观测地区被动微波辐射计Ka频水平极化数据。

图5是本发明实施例1林地观测地区被动微波辐射计K频水平极化组分亮温数据。

图6是本发明实施例1林地观测地区被动微波辐射计Ka频水平极化组分亮温数据。

图7是本发明实施例1原始被动微波数据和被动微波组分亮温数据进行雪深反演结果对比。

图8是本发明方法步骤流程图。

具体实施方式

下面以具体实施例的形式对本发明技术方案进一步说明。

实施例1

本发明利用FY-3BMWRI被动微波遥感数据和MODIS土地覆盖分类数据产品，结合多频双极化林地积雪被动微波混合像元分解方法，实现了2012年1月中国黑龙江省伊春市林地积雪被动微波混合像元的有效分解。

具体步骤如下：

(一)林地分类数据的获取和重新分类

本发明使用的土地分类数据是MODIS中分辨率光谱遥感数据的土地覆盖类型产品。在MODIS官方网站下载中分辨率MODIS光谱遥感数据的土地覆盖类型产品。MODIS陆地研究小组开发了MODISAqua和Terra卫星数据合成的年度土地覆盖分类产品MCD12Q1(空间分辨率500m)，选择MCD12Q1分类产品中的NPP(净初级生产力)土地覆盖分类数据集，时间是2012年。

NPP(净初级生产力)土地覆盖分类数据集将下垫面共分9类，包括：水，常绿针叶植被，常绿阔叶植被，落叶针叶植被，落叶阔叶植被，一年生阔叶植被，一年生草本植被，无植被的土地和城市。

根据观测地区获得其数据主要包括：常绿针叶植被，常绿阔叶植被，落叶针叶植被，落叶阔叶植被，一年生阔叶植被，一年生草本植被和城市，因此将分类数据产品重新分类具体如下：

(a)将原分类数据中定义为常绿针叶植被和落叶针叶植被的类型统一重新定义为针叶林类型N；

(b)将原分类数据中定义为常绿阔叶植被、落叶阔叶植被和一年生阔叶植被的类型统一重新定义为阔叶林类型B；

(c)一年生草本植被类型定义为Q；

(d)城市类型定义为U。

采用上述分类方法，可得到研究地区阔叶林类型、阔叶林类型、一年生草本植被和城市类型四种下垫面分类数据L^N、L^B、L^Q和L^U。

(二)建立多频双极化林地积雪被动微波混合像元分解模型

风云三号B星上装载的微波成像仪(MWRI)为我国第一个星载微波遥感仪器,扫描方式为圆锥扫描，其设计频率为10.65～150GHz，其中150GHz为试验通道。每个频率都有垂直和水平两种不同极化模式，这些频率的遥感成像能提供全天候、全天时地表温度、土壤水分、洪涝干旱、积雪深度、台风结构、大气含水量等丰富的信息。从网站下载MWRI被动微波遥感数据，其空间分辨率为10km。该产品已经过标定、大气校正、地理校正和标准化预处理。

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

$>T_F^M(x_0,y_0)=\int \int \bar{G}(x,y)T_F^C(x,y)dxdy---\left(1\right)$>

公式(1)中，是辐射计天线测量的混合亮温值，是在地面位置(x,y)实际亮温值，下脚标F是被动微波观测频段，是天线增益函数，积分覆盖整个足印的椭圆区域，天线增益函数的公式定义为：

$>\bar{G}(x,y)=\frac{ln2}{\pi }·\exp [-ln2·(\frac{x^2}{a^2}+\frac{y^2}{b^2})]---\left(2\right)$>

a和b表示被动微波数据中频段对应的足印，定义为a＝5km和b＝5km。

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

$>T_F^C(x,y)=\underset{i=1}{\overset{4}{\Sigma }}{L}^i(x,y)T_i---\left(3\right)$>

其中T_i代表不同下垫面(组分)亮温，Lⁱ代表第i个下垫面分类数据。

星载被动微波K频(18.7GHz)和Ka频(36.5GHz)是常用判识积雪的有效频段，考虑到观测地区的主要由阔叶林、针叶林、一年生草本植被和城市下垫面组成，地面位置(x,y)处的林地积雪被动微波在K频和Ka频的辐射值可以表示为：

$>\left(\begin{array}{c}{T}_{KH}^C(x,y)\\ T_{KV}^C(x,y)\\ T_{KaH}^C(x,y)\\ T_{KaV}^C(x,y)\end{array}\right)=L^B(x,y)\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+L^N(x,y)\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+L^Q(x,y)\left(\begin{array}{c}{T}_{KH}^Q\\ T_{KV}^Q\\ T_{KaH}^Q\\ T_{KaV}^Q\end{array}\right)+L^U(x,y)\left(\begin{array}{c}{T}_{KH}^U\\ T_{KV}^U\\ T_{KaH}^U\\ T_{KaV}^U\end{array}\right)+ϵ---\left(4^{,}\right)$>

其中T代表亮温值；L代表下垫面分类数据，若(x,y)位置存在该下垫面，则L(x,y)＝1；否则，L(x,y)＝0。上脚标C、B、N、Q和U分别代表观测位置下垫面类型、阔叶林、针叶林、一年生草本植被和城市下垫面。下脚标K和Ka分别代表被动微波K和Ka频段，其中下脚标H和V分别代表水平极化和垂直极化。ε代表残差量。L^B、L^N、L^Q和L^U分别代表地面位置(x,y)处阔叶林、针叶林、一年生草本植被和城市下垫面的分布，且有：

L^B(x,y)+L^N(x,y)+L^Q(x,y)+L^U(x,y)＝1(5’)

星载被动微波辐射计在不同频段(F)和极化方式(V、H)观测到的的林地积雪亮温值可以看成是被动微波天线和地面实际下垫面类型亮温的卷积，将(4)代入(1)有：

$>\left(\begin{array}{c}{T}_{KH}^M(x_0,y_0)\\ T_{KV}^M(x_0,y_0)\\ T_{KaH}^M(x_0,y_0)\\ T_{KaV}^M(x_0,y_0)\end{array}\right)=\int \int \bar{G}(x,y)\left(L^B\right(x,y)\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+L^N(x,y)\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+L^S(x,y\left)\left(\begin{array}{c}{T}_{KH}^S\\ T_{KV}^S\\ T_{KaH}^S\\ T_{KaV}^S\end{array}\right)\right)dxdy+ϵ---\left(6^{,}\right)$>

令α是天线增益函数和下垫面分类数据L的卷积，有：

$>{\alpha }^C(x_0,y_0)=\int \int \bar{G}(x,y)L^C(x,y)dxdy,$>C＝B,N,Q,U(7)

则(6)式变为：

$>\left(\begin{array}{c}{T}_{KH}^M\\ T_{KV}^M\\ T_{KaH}^M\\ T_{KaV}^M\end{array}\right)={\alpha }^B\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)+{\alpha }^N\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)+{\alpha }^Q\left(\begin{array}{c}{T}_{KH}^Q\\ T_{KV}^Q\\ T_{KaH}^Q\\ T_{KaV}^Q\end{array}\right)+{\alpha }^U\left(\begin{array}{c}{T}_{KH}^U\\ T_{KV}^U\\ T_{KaH}^U\\ T_{KaV}^U\end{array}\right)+ϵ---\left(8^{,}\right)$>

公式(8)中包括4个方程，已知4个观测到的林地混合亮温值和4个不同下垫面分类计算得到的α值，需要得到16个未知量,即4个阔叶林亮温，4个针叶林亮温，4个一年生草本植被亮温和4个城市亮温，这属于线性欠定行方程组求解问题。

(三)基于动态窗口数据选择策略的欠定性方程组求解

理论上，公式(7)应该有无数个解。因此，考虑采用以下多点观测约束求解方法:

$>\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right)=\left(\begin{array}{c}{T}_{KH}^B\\ T_{KV}^B\\ T_{KaH}^B\\ T_{KaV}^B\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\end{array}\right]+\left(\begin{array}{c}{T}_{KH}^N\\ T_{KV}^N\\ T_{KaH}^N\\ T_{KaV}^N\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^N& ...& {\alpha }_p^N\end{array}\right]+\left(\begin{array}{c}{T}_{KH}^Q\\ T_{KV}^Q\\ T_{KaH}^Q\\ T_{KaV}^Q\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^Q& ...& {\alpha }_p^Q\end{array}\right]+\left(\begin{array}{c}{T}_{KH}^U\\ T_{KV}^U\\ T_{KaH}^U\\ T_{KaV}^U\end{array}\right)\left[\begin{array}{ccc}{\alpha }_1^U& ...& {\alpha }_p^U\end{array}\right]+E---\left(9^{,}\right)$>

其中p是观测点的个数，E代表残差矩阵。这里假设在观测点周围内所有的同类下垫面有相同的亮值，该条件可以作为公式(8)的约束。作为超定线性系统的公式(8)可以写成：

$>\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right)=\left(\begin{array}{cccc}{T}_{KH}^B& T_{KH}^N& T_{KH}^Q& T_{KH}^U\\ T_{KV}^B& T_{KV}^N& T_{KV}^Q& T_{KH}^U\\ T_{KaH}^B& T_{KaH}^N& T_{KaH}^Q& T_{KH}^U\\ T_{KaV}^B& T_{KaV}^N& T_{KaV}^Q& T_{KH}^U\end{array}\right)\left(\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\\ {\alpha }_1^N& ...& {\alpha }_p^N\\ {\alpha }_1^Q& ...& {\alpha }_p^Q\\ {\alpha }_1^U& ...& {\alpha }_p^U\end{array}\right)+E---\left({10}^{,}\right)$>

考虑最小二乘法：

$>T^{*}=\underset{T>0}{\mathrm{argmin}}\left|\right|Y-TA\left|\right|---\left(11\right)$>

其中

$>Y=\left(\begin{array}{ccc}{T}_{KH}^M\left(1\right)& ...& T_{KH}^M\left(p\right)\\ T_{KV}^M\left(1\right)& ...& T_{KV}^M\left(p\right)\\ T_{KaH}^M\left(1\right)& ...& T_{KaH}^M\left(p\right)\\ T_{KaV}^M\left(1\right)& ...& T_{KaV}^M\left(p\right)\end{array}\right),---\left(12\right)$>

$>T=\left(\begin{array}{cccc}{T}_{KH}^B& T_{KH}^N& T_{KH}^Q& T_{KH}^U\\ T_{KV}^B& T_{KV}^N& T_{KV}^Q& T_{KH}^U\\ T_{KaH}^B& T_{KaH}^N& T_{KaH}^Q& T_{KH}^U\\ T_{KaV}^B& T_{KaV}^N& T_{KaV}^Q& T_{KH}^U\end{array}\right)---\left({13}^{,}\right)$>

和

$>A=\left(\begin{array}{ccc}{\alpha }_1^B& ...& {\alpha }_p^B\\ {\alpha }_1^N& ...& {\alpha }_p^N\\ {\alpha }_1^Q& ...& {\alpha }_p^Q\\ {\alpha }_1^U& ...& {\alpha }_p^U\end{array}\right).---\left({14}^{,}\right)$>

此外，对于某频段同一位置求解得到的被动微波组分亮温，还应该满足以下条件：

$>T_{KV}^C>T_{KH}^C---\left(15\right)$>

$>T_{KaV}^C>T_{KaH}^C---\left(16\right)$>

理论上，公式(11)可以通过非负最小二乘法求解。对于解混问题，将其转换为普通的最小二乘法其满足亮温解正值的条件。

T＝(A^tA)^-1A^tY(17)

分解误差为：

E＝T-(A^tA)^-1A^tY(18)

对于公式(10)中p值选取，采用动态窗口数据选择策略：

首先，考虑到微波像元的空间相关性，以待分解的位置(x₀,y₀)处的F频段的被动微波混合像元T_F^M为中心，选取其周围8邻域窗口的被动微波混合像元(如图1所示)和周围4邻域窗口的被动微波混合像元(如图2所示)作为方程组输入数据进行求解，F分别取K和Ka频段微波数据。

对应两种输入数据的方案为：

(a)8邻域窗口数据输入：按照图1,将含当前像元在内的9像元的多频被动微波数据转变为一维矢量代入公式(12)。如果计算得到的结果满足公式(15)和公式(16),获得合理求解，得到分解后各下垫面对应的组分亮温数据T₁，得到误差数据E₁；

(b)4邻域窗口数据输入：按照图2,将含当前像元在内的5个多频被动微波数据转变为一维矢量代入公式(12)。如果计算得到的结果满足公式(15)和公式(16),获得合理求解，得到分解后各下垫面对应的组分亮温数据T₂，得到误差数据E₂；

设计动态窗口数据选择策略，求解欠定性方程组。根据(a)、(b)方案的求解结果，考虑以下四种情况，决定最后输出的组分亮温：

情况一：选择(a)方案，若求解不合理，则采用(b)方案；若(b)方案求解仍不合理，说明输入的混合像元不能分解，保留原像元的亮温值；

情况二：选择(a)方案，若求解不合理，则采用(b)方案；若(b)方案求解合理，输出组分亮温值T₂和误差数据E₂，最终组分亮温值是T₂、误差数据是E₂；

情况三：选择(a)方案，若求解合理，输出组分亮温值T₁和误差数据E₁；继续采用(b)方案，若(b)方案求解不合理，则最终组分亮温值是T₁、误差数据是E₁；

情况四：选择(a)方案，若求解合理，输出组分亮温值T₁和误差数据E₁；继续采用(b)方案，若(b)方案求解合理，输出组分亮温值T₂和误差数据E₂。分别计算各下垫面对应的组分的误差E₁总和记为E₁′，误差E₂的总和记为E₂′，比较E₁′和E₂′的大小。若E₁′＞E₂′，则最终组分亮温值是T₂、误差数据是E₂；若E₁′＜E₂′，则最终组分亮温值是T₁、误差数据是E₁。

效果验证：

根据Chang的微波雪深反演方法，假设积雪为均一、单层的干雪，雪密度为0.3g/cm³，雪粒径为0.3mm，雪深反演公式如下：

SD＝C×(T_KH-T_KaH)(19)

其中SD为积雪深度(cm)；C为回归系数，通常取1.59；T_KH和T_KaH别为被动微波K和Ka频段的水平极化数据。由中国黑龙江省伊春市林地观测地区的被动微波辐射计K和Ka频段水平极化数据(如图3和图4所示)，利用本发明提出的多频双极化林地积雪被动微波混合像元分解方法，可以得到K和Ka频段分解后的水平极化组分亮温数据(如图5和图6所示)。

采用原始被动微波K和Ka频段水平极化数据和对应的混合像元分解后的组分亮温数据，在测试点位置所在像元并根据其实际下垫面土地分类所属组分亮温数据进行Chang雪深反演的结果(如图7所示)。由图7可见，Chang雪深反演方法中，采用分解后的组分亮温数据进行雪深反演得到的结果，相对于原始数据雪深反演的结果，更接近于实测林地的雪深数据。采用本发明提出的多频双极化林地积雪被动微波混合像元分解方法，可以获得针叶林、阔叶林下垫面的积雪深度，有效地提高后期雪深反演的精度。