农业可持续发展与农业资源的高效配置是世界各国政府关注的重点问题，也是联合国“2030年可持续发展议程”中非常重要的组成部分。由于全球经济仍未从新冠疫情带来的巨大冲击中恢复，部分国家和地区的局势动荡也为全球粮食供应带来了极大的不确定性。2024年全球农业可持续发展状况不容乐观，科学的农业规划与农业资源的有效统筹仍然是各国政府关注的重点。遥感技术因具有较大的覆盖范围和较为稳定的数据来源，能够大大提高农业资源调查与规划的时间和经济效率，已经成为我国农业发展的重要组成部分。光学遥感图像已经在植被指数监测^[1]、产量估计^[2]、含水量估计^[3]、病虫害监测^[4-5]等领域取得了很好的研究进展。但是，光学或高光谱遥感数据易受云雨天气影响，在农作物的关键生长期内很难及时获取充足、清晰的影像。相对地，合成孔径雷达(synthetic aperture radar，SAR)作为一种不受云雨天气影响的高分辨率遥感技术，可对农作物进行全天时、全天候的监测，微波的穿透性使得SAR信号可以穿过植被表面，获取来自植被冠层以及土壤的信息，为农作物实时监测提供了可靠的数据保障和技术支持^[6-8]，在无法获取稳定光学遥感数据的地区成为进行农业遥感观测的重要手段。

近年来，随着雷达技术的发展，星载雷达传感器的性能和成像模式不断发展进步，具备了获取多分辨率、多成像模式、多极化通道SAR图像的能力。相较于单极化SAR，全极化SAR包含更加全面的地表散射信息，显著提高了SAR系统进行农作物精确分类的能力^[9]。在农业应用领域，多时相SAR信息是刻画农作物的散射特性随生长而变化的关键。Li等^[10]利用支持向量机对4景RADARSAT-2全极化图像进行水稻种植区提取研究，并根据极化分解特征分析水稻生长过程中的散射机制变化。Jiao等^[11]获取了19景RADARSAT-2多时相全极化SAR图像进行加拿大东北部Ontario地区的地表覆盖分类，证明Freeman-Durden分解和H/A/α分解对作物识别的有效性。McNairn等^[12]提取多时相TerraSAR-X和RADARSAT-2图像后向散射系数，用决策树方法进行玉米和大豆的生长早期分类。Demarez等^[13]利用随机森林方法和增量分类实验开展法国西南部的灌溉作物识别研究。Wang等^[14]利用高分3号全极化SAR对比不同极化分解特征在玉米和棉花识别方面的表现。Xie等^[15]利用7景RADARSAT-2图像评估物理约束通用模型的极化分解方法在全极化SAR农作物分类中的表现，结合前向特征筛选策略和随机森林分类器进行分析，验证了该极化分解方法的有效性。然而，这些研究大多聚焦农作物的极化散射机制，直接使用机器学习方法，将多时相的极化特征视为彼此独立的通道，忽略了极化特征的时间相关性。近年来引发广泛关注的深度学习方法为极化SAR的农作物分类提供了新的思路。深度卷积神经网络或全卷积网络具备从空间上下文中学习高维分布特征的能力，但是忽视了作物在多时相极化SAR中的时间关联性；循环神经网络方法能提取特征的时序规律，但无法兼顾极化SAR多样的极化特征，也难以保持像素间的空间关系^[16-17]。虽然已有一些研究提出了新型网络结构以从多时相极化SAR的高维极化特征中学习时间和空间信息^[18-19]，但是由于深度学习方法需要较大的样本数量，仍需对这些方法进行进一步的验证。目前，能够兼顾农作物散射信息和时序变化规律的多时相极化SAR分类方法仍有待深入研究。

全极化SAR虽然能提供丰富的散射信息，但数据覆盖范围较小，有时无法满足农业监测的实际需求。Souyris等^[20]提出简缩极化SAR的概念，设计了一种发射45°方向线极化波，接收水平和垂直回波信号的特殊的双极化SAR信号组合；Stacy等^[21]提出一种双圆极化模式，该模式发射左旋或者右旋圆极化波，并接收左旋和右旋的圆极化回波；Raney^[22]提出一种混合极化模式，又被称为CTLR模式(circular transmit linear receive)，该模式发射圆极化波，接收水平和垂直回波。简缩极化SAR本质上是一种双极化系统，能在一定程度上兼顾极化信息和数据覆盖范围，为大范围农作物监测提供了新的机遇，已引起研究人员的广泛关注^[23-25]。Brisco等^[26]获取了5景全极化影像，对比简缩极化、双极化与全极化SAR的水稻种植区和湿地提取的精度，结果表明混合模式的简缩极化SAR能取得比HH/HV双极化SAR更好的分类结果。Xie等^[27]证明了用简缩极化SAR的线极化度和线极化比特征进行水稻提取的可行性。Uppala等^[28]用1景RISAT-1卫星获取的真实简缩极化SAR影像的后向散射系数和m-χ分解结果进行了玉米种植区提取。Robertson等^[29]利用加拿大RADARSAT星座任务在2020年6—7月获取的17景C波段简缩极化SAR数据，采用随机森林方法对草场、大麦、小麦等7类农作物进行分类评估，发现能取得与光学遥感图像相当的分类精度，并讨论了作物极化度随物候的变化情况。Guo等^[30]从7景RADARSAT-2模拟的简缩极化SAR数据中提取了20维极化参数进行水稻分类研究，发现引入简缩极化参数在时间尺度上的变化信息对水稻分类更加有利，机器学习分类的最优精度可以达到97 %。目前的研究表明，在农作物分类应用中，简缩极化能够取得与全极化SAR相当的结果^[31]。但是，简缩极化SAR的极化分解特征与全极化SAR存在一定差异，因而不能将其简单地视作全极化SAR的替代^[32-33]。Xie等^[34-35]提出一个基于H/α分解的多模式极化SAR分类方法，并针对中国南方的典型作物提出一个基于表面散射和散射熵的多模式极化SAR农作物分类方法。这些研究虽然初步建立了多模式极化SAR分类统一框架，但仍然局限于单时相情况。目前，能兼顾简缩极化的多时相极化SAR农作物分类方法仍有待进一步研究。

针对上述问题，本文提出一个适用于全极化和简缩极化SAR的多时相极化SAR分类方法，该方法引入张量空间的多线性主成分分析(multilinear principle component analysis，MPCA)方法^[36]，从多时相极化SAR的协方差矩阵中提取有利于分类的核心信息，最后用决策树方法对降维特征进行分类。本文方法的创新之处在于提出了适用于多模式极化SAR的多时相张量表示范式，该范式以多时相极化SAR的复数协方差矩阵作为输入，能同时适用于全极化和简缩极化SAR；在张量域进行降维，顺应了多时相重复观测极化SAR数据之间的时间顺序逻辑，其降维结果能够一定程度上保持复数极化SAR隐含的相位信息，以及农作物随着生长物候改变而产生的极化散射变化信息。此外，为降低分类结果中的斑噪影响，采取目标级分割策略，将简单线性迭代聚类(simple linear iterative cluster，SLIC)超像素分割引入多时相、多模式极化SAR，以提升分类结果中同质地块的一致性。通过目标级分类策略和张量空间的多时相极化特征降维，实现了空间信息、极化信息和时序变化信息的有机结合。

1 基于张量表示的多时相多模式极化SAR农作物分类方法

图 1展示了本文方法的主要步骤：首先，提取多时相多模式极化SAR的协方差矩阵，使用SLIC方法进行图像超像素分割；其次，提取目标样本的极化特征并进行张量化表示；再次，使用MPCA方法对张量目标样本集进行特征降维；最后，用决策树分类器对降维特征进行分类，获取农作物分类结果。

1.1 多时相极化SAR超像素分割

进行超像素分割之前，首先对N景覆盖同一区域的多时相多模式极化SAR图像进行必要的预处理，包括配准、多视和滤波等操作，并计算各时相数据的协方差矩阵。在互易假设条件下，Lexicographic基下的全极化SAR散射矢量k可表示为

$ \boldsymbol{k}=\left[\begin{array}{lll} S_{\mathrm{HH}} & \sqrt{2} S_{\mathrm{HV}} & S_{\mathrm{VV}} \end{array}\right]^{\mathrm{T}} . $

(1)

其中：S_HH、S_HV、S_VV分别表示全极化Sinclair散射矩阵HH、HV和VV极化通道分量，上标T表示转置操作。对于简缩极化SAR，其散射矢量可以表示为全极化SAR散射矢量元素的线性组合，因此，可以直接用全极化SAR模拟简缩极化SAR。π/4模式的散射矢量 k_π/4可以表示为

$ \boldsymbol{k}_{\pi / 4}=\left[\begin{array}{ll} S_{\mathrm{HH}} & S_{\mathrm{HV}} \\ S_{\mathrm{VH}} & S_{\mathrm{VV}} \end{array}\right] \cdot \frac{1}{\sqrt{2}}\left[\begin{array}{l} 1 \\ 1 \end{array}\right]=\frac{1}{\sqrt{2}}\left[\begin{array}{l} S_{\mathrm{HH}}+S_{\mathrm{HV}} \\ S_{\mathrm{VV}}+S_{\mathrm{HV}} \end{array}\right] . $

(2)

CTLR模式与双圆极化模式都发射圆极化波，本质上提供了相同的极化信息^[26]。以发射右旋圆极化波的CTLR模式为例，其散射矢量 k_CTLR可以表示为

$ \boldsymbol{k}_{\mathrm{CTLR}}=\left[\begin{array}{ll}S_{\mathrm{HH}} & S_{\mathrm{HV}} \\ S_{\mathrm{VH}} & S_{\mathrm{VV}}\end{array}\right] \cdot \frac{1}{\sqrt{2}}\left[\begin{array}{c}1 \\ -\mathrm{i}\end{array}\right]=\frac{1}{\sqrt{2}}\left[\begin{array}{c}S_{\mathrm{HH}}-\mathrm{i} S_{\mathrm{HV}} \\ S_{\mathrm{HV}}-\mathrm{i} S_{\mathrm{VV}}\end{array}\right]. $

(3)

将散射矢量与其自身复共轭相乘再进行空间平均即可得到协方差矩阵，全极化、π/4模式和右旋CTLR模式的协方差矩阵 C_FP、C_π/4、C_CTLR可分别表示为

$ \boldsymbol{C}_{\mathrm{FP}}=\boldsymbol{k} \cdot \boldsymbol{k}^{\mathrm{T}}=\left[\begin{array}{ccc} \left.\left.\langle | S_{\mathrm{HH}}\right|^{2}\right\rangle & \sqrt{2}\left\langle S_{\mathrm{HH}} S_{\mathrm{HV}}^{*}\right\rangle & \left\langle S_{\mathrm{HH}} S_{\mathrm{VV}}^{*}\right\rangle \\ \sqrt{2}\left\langle S_{\mathrm{HV}} S_{\mathrm{HH}}^{*}\right\rangle & \left.\left.2\langle | S_{\mathrm{HV}}\right|^{2}\right\rangle & \sqrt{2}\left\langle S_{\mathrm{HV}} S_{\mathrm{VV}}^{*}\right\rangle \\ \left\langle S_{\mathrm{VV}} S_{\mathrm{HH}}^{*}\right\rangle & \sqrt{2}\left\langle S_{\mathrm{VV}} S_{\mathrm{HV}}^{*}\right\rangle & \left.\left.\langle | S_{\mathrm{VV}}\right|^{2}\right\rangle \end{array}\right], $

(4)

$ \boldsymbol{C}_{\pi / 4}=\left\langle\boldsymbol{k}_{\pi / 4} \cdot \boldsymbol{k}_{\pi / 4}^{\mathrm{T}}\right\rangle=\frac{1}{2}\left[\begin{array}{cc} \left.\langle | S_{\mathrm{HH}}+\left.S_{\mathrm{HV}}\right|^{2}\right\rangle & \left\langle\left(S_{\mathrm{HH}}+S_{\mathrm{HV}}\right)\left(S_{\mathrm{VV}}+S_{\mathrm{HV}}\right)^{*}\right\rangle \\ \left\langle\left(S_{\mathrm{VV}}+S_{\mathrm{HV}}\right)\left(S_{\mathrm{HH}}+S_{\mathrm{HV}}\right)^{*}\right\rangle & \left.\langle | S_{\mathrm{VV}}+\left.S_{\mathrm{HV}}\right|^{2}\right\rangle \end{array}\right], $

(5)

$ \boldsymbol{C}_{\mathrm{CTLR}}=\left\langle\boldsymbol{k}_{\mathrm{CTLR}} \cdot \boldsymbol{k}_{\mathrm{CTLR}}^{\mathrm{T}}\right\rangle=\frac{1}{2}\left[\begin{array}{cc} \left.\langle | S_{\mathrm{HH}}-\left.\mathrm{i} S_{\mathrm{HV}}\right|^{2}\right\rangle & \left\langle\left(S_{\mathrm{HH}}-\mathrm{i} S_{\mathrm{HV}}\right)\left(S_{\mathrm{HV}}-\mathrm{i} S_{\mathrm{VV}}\right)^{*}\right\rangle \\ \left\langle\left(S_{\mathrm{HV}}-\mathrm{i} S_{\mathrm{VV}}\right)\left(S_{\mathrm{HH}}-\mathrm{i} S_{\mathrm{VV}}\right)^{*}\right\rangle & \left.\langle | S_{\mathrm{HV}}-\left.\mathrm{i} S_{\mathrm{VV}}\right|^{2}\right\rangle \end{array}\right] . $

(6)

同一农作物种植区内部的像素往往具有相似的后向散射特征，可以看作匀质区域。然而受到SAR图像斑噪和空间随机性的影响，分类结果图像容易出现“噪声”。为了消除由斑噪和时间随机性带来的分类误差，本文采用目标级分类策略，利用SLIC方法将散射特性相似的相邻像素划分为同一目标，使得目标内部分类结果保持一致，以抑制农作物分类结果中的噪声现象。原始的SLIC方法适用于光学图像，利用像素在CIE-Lab颜色空间的距离和像素点间的欧式距离实现超像素分割。Liu等^[37]和Qin等^[38]用Wishart距离代替颜色空间距离，将SLIC方法引入极化SAR图像分割。本文进一步将极化SAR的SLIC分割方法推广至多时相情况。在多时相情况下，农作物的后向散射和极化特征可能随生长过程发生变化，造成目标分割结果的不一致。为此，将多时相极化SAR数据视作一个整体，采用N个时相Wishart距离均值d_w而非单一时相的Wishart距离作为距离度量进行SLIC分割，因此，不需要对多景数据进行重复分割，而是直接获取统一的多时相目标分割边界。获取超像素目标边界后，以每个目标内所有像素的协方差矩阵均值代表目标的协方差矩阵，每个目标的协方差矩阵将仍然是复数的埃尔米特半正定矩阵。多时相极化SAR超像素分割过程的具体步骤如下：

1) 超像素中心初始化，给定初始超像素的大小S(S=2R+1，R为超像素半径参数)，将图像分割为规则的矩形区域，每个矩形区域的大小为S×S，矩形区域的中心像素即为超像素初始中心。

2) 初始超像素中心位置修正，为避免超像素中心位于边缘或噪声像素，计算初始超像素中心3×3邻域的梯度，并选取梯度最低的像素作为修正后的初始超像素中心。

3) 超像素的迭代，以每一个超像素中心为原点，以2R为搜索半径，计算搜索范围内每个像素与超像素中心之间的距离D，每个像素将被周围的多个超像素搜索，并最终被标记为距离最近的超像素。距离D计算公式^[37-38]如下

$ D=\sqrt{\left(\frac{\bar{d}_{\mathrm{w}}}{m}\right)^2+\left(\frac{d_{\mathrm{s}}}{S}\right)^2}, $

(7)

$ \begin{gathered} d_{\mathrm{w}}(i, J)=\ln \left(\frac{\operatorname{det}\left(\overline{\boldsymbol{C}}_J\right)}{\operatorname{det}\left(\boldsymbol{C}_i\right)}\right)+\operatorname{trace}\left(\frac{\boldsymbol{C}_i}{\overline{\boldsymbol{C}}_J}\right)-q, \\ \bar{d}_{\mathrm{w}}=\frac{1}{N} \sum\nolimits_{n=1}^N d_{\mathrm{w}}^n, \end{gathered} $

(8)

$ d_{\mathrm{s}}(i, J)=\sqrt{\left(r_{i}-r_{J}\right)^{2}+\left(c_{i}-c_{J}\right)^{2}}. $

(9)

式中：C_i表示像素i的协方差矩阵；C_J表示目标J的平均协方差矩阵，对于由k个像素构成的目标J，有$ \overline{\boldsymbol{C}}_J=\frac{1}{k} \sum_{i=1}^k \boldsymbol{C}_i$；D由Wishart距离和欧氏距离构成，其中，d_w(i, J)表示第i个像素和第J个目标中心之间的Wishart距离；m为Wishart距离的权重参数；det()表示矩阵的行列式；trace()表示矩阵的迹；q为维数参数，对于全极化SAR有q=3，对于简缩极化SAR有q=2；d_s(i, J)表示第i个像素和第J个目标中心之间的空间距离，其坐标(r_i, c_i)和(r_J, c_J)分别表示第i个像素和第J个目标中心在图像上的行列位置。

1.2 多时相极化SAR目标级样本的张量表示

一个张量$ \boldsymbol{X} \in \mathbb{R}^{I_{1} \times I_{2} \times \cdots \times I_{N}}$可以被视为一个N维数组，N被称为张量的阶数或者模数，它的第(i₁, i₂, …, i_n)个元素可表示为x_{(i1, i2, …, in)}。在数据处理中应用最为广泛的矢量和矩阵即可被分别视为1阶张量和2阶张量，3阶张量表示多个矩阵的集合，而4阶张量又表示多个3阶张量的组合，以此类推。极化SAR数据的协方差矩阵为埃尔米特半正定矩阵，属于黎曼流形，因此可被视为2阶张量^[39]。多时相全极化SAR中每个像素的协方差矩阵可被视为一个多层矩阵，随着重复观测图像数量的增加，协方差矩阵的层数也随之增加，因此每个像素可以表示成一个3阶张量，阶数I₁=I₂=3，I₃=N，其中前2阶表示协方差矩阵的维数，第3阶表示影像数量。类似地，一个由k个像素构成的目标也可以表示为3阶张量的形式，前2阶表示目标协方差矩阵维数，第3阶表示SAR图像数量。因此，多时相全极化SAR的目标样本集可表示为$ \boldsymbol{X}=\left\{\boldsymbol{X}_{1}, \boldsymbol{X}_{2}, \cdots, \boldsymbol{X}_{M}\right\}, M$表示目标样本总个数，$ \boldsymbol{X}_{m} \in \mathbb{R}^{3 \times 3 \times N}(m=1, 2, \cdots, M)$代表第m个目标样本，Cⁿ表示该目标样本在第n个时相的平均协方差矩阵

$ \begin{gather*} \boldsymbol{X}_{m}=\left[\begin{array}{llll} \overline{\boldsymbol{C}}^{1} & \overline{\boldsymbol{C}}^{2} & \cdots & \overline{\boldsymbol{C}}^{N} \end{array}\right]_{m} , \\ \overline{\boldsymbol{C}}^{n}=\frac{1}{k} \sum\nolimits_{i=1}^{k} \boldsymbol{C}_{i}(n=1, 2, \cdots, N) . \end{gather*} $

(10)

这种表示方式可以很容易地推广至简缩极化SAR，其每个目标样本可表示为$ \boldsymbol{X}_m \in \mathbb{R}^{2 \times 2 \times N}$。

1.3 张量空间MPCA特征降维

一个阶数为N的张量，固定其任意N-2维即可获取一个矩阵，称为张量切片，固定其任意N-1维即可获取一个矢量，称为张量纤维。按列优先方式将张量X的第m阶纤维重新排列成矩阵，即可得到张量的m阶展开矩阵$ \boldsymbol{X}_{(m)} \in \mathbb{R}^{I_m \times\left(I_1 \times \cdots \times I_{m-1} \times I_{m+1} \times \cdots \times I_N\right)} 。\boldsymbol{X}_{(m)}$与矩阵$ \boldsymbol{A} \in \mathbb{R}^{I_A \times I_m}$相乘可得一新张量B，其维数为I₁×…×I_m-1×I_A×I_m+1×…×I_N，可记作

$ \begin{equation*} \boldsymbol{B}=\boldsymbol{X} \times{ }_{m} \boldsymbol{A} \text {. } \end{equation*} $

(11)

张量B的各元素为

$ \begin{gather*} b_{\left(i_{1}, i_{2}, \cdots, i_{m-1}, i_{a}, i_{m+1}, \cdots, i_{N}\right)}=\sum\nolimits_{i_{m}=1}^{I_{m}} x_{\left(i_{1}, i_{2}, \cdots, i_{n}\right)} \cdot a_{\left(i_{a}, i_{m}\right)}, \\ i_{a}=1, 2, \cdots, I_{A} . \end{gather*} $

(12)

张量B的m阶展开矩阵可表示为

$ \begin{equation*} \boldsymbol{B}_{(m)}=\boldsymbol{A} \cdot \boldsymbol{X}_{(m)} . \end{equation*} $

(13)

任意张量$ \boldsymbol{X} \in \mathbb{R}^{I_{1} \times I_{2} \times \cdots \times I_{N}}$均可被分解成一个核心张量P与一系列矩阵U相乘的形式，且过程可逆

$ \boldsymbol{X}=\boldsymbol{P} \times{ }_{1} \boldsymbol{U}^{(1)} \times{ }_{2} \boldsymbol{U}^{(2)} \cdots \times{ }_{N} \boldsymbol{U}^{(N)}, $

(14)

$ \boldsymbol{P}=\boldsymbol{X} \times{ }_{1}\left(\boldsymbol{U}^{(1)}\right)^{\mathrm{T}} \times{ }_{2}\left(\boldsymbol{U}^{(2)}\right)^{\mathrm{T}} \cdots \times{ }_{N}\left(\boldsymbol{U}^{(N)}\right)^{\mathrm{T}}. $

(15)

其中: P称作核心张量; $ \boldsymbol{U}_{n=1, 2, \cdots, N}^{(n)}$表示N个维数为I_n×I_n正交方阵，称为分解矩阵。张量X的m阶展开可以表示为

$ \begin{gather*} \boldsymbol{X}_{(m)}=\boldsymbol{U}^{(m)} \cdot \boldsymbol{P}_{(m)} \cdot\left(\boldsymbol{U}^{(m+1)} \otimes \cdots \otimes \boldsymbol{U}^{(N)}\otimes\right. \\ \left.\boldsymbol{U}^{(1)} \otimes \cdots \otimes \boldsymbol{U}^{(m-1)}\right)^{\mathrm{T}} . \end{gather*} $

(16)

其中: P_(m)为核心张量P的m阶展开矩阵，$ \otimes$表示矩阵的Kronecker积。

MPCA算法^[36]是二维主成分分析(principle component analysis，PCA)在张量域的拓展，其目的是找到一系列投影矩阵$ \boldsymbol{U}_{n=1, 2, \cdots, N}^{(n)} \in \mathbb{R}^{I_{n} \times L_{n}}$，$ \left(L_{n}<I_{n}\right)$，将原始张量样本集$ \boldsymbol{X}=\left\{\boldsymbol{X}_{1}, \boldsymbol{X}_{2}, \cdots\right.\left.\boldsymbol{X}_{M}\right\}, \boldsymbol{X}_{m} \in \mathbb{R}^{I_{1} \times I_{2} \times \cdots \times I_{N}}$投影到新的张量空间，构造包含原始样本集核心信息的新样本集$ \boldsymbol{Y}=\left\{\boldsymbol{Y}_{1}\right.\left.\boldsymbol{Y}_{2}, \cdots, \boldsymbol{Y}_{M}\right\}, \boldsymbol{Y}_{m} \in \mathbb{R}^{L_{1} \times L_{2} \times \cdots \times L_{N}}$，以实现特征降维。其中每个样本Y_{m=1, 2, …，M}均为原始样本X_{m=1, 2, …，M}在矩阵$ \boldsymbol{U}_{n=1, 2, \cdots, N}^{(n)}$下的投影

$ \begin{equation*} \boldsymbol{Y}_{m}=\boldsymbol{X}_{m} \times{ }_{1} \boldsymbol{U}^{(1)} \times{ }_{2} \boldsymbol{U}^{(2)} \cdots \times{ }_{N} \boldsymbol{U}^{(N)} . \end{equation*} $

(17)

样本集Y={Y₁, Y₂, …, Y_M}的散度可以定义为

$ \varphi_{\boldsymbol{Y}}=\sum\nolimits_{m=1}^{M}\left\|\boldsymbol{Y}_{m}-\overline{\boldsymbol{Y}}\right\|_{F}^{2}, \overline{\boldsymbol{Y}}=\frac{1}{M} \sum\nolimits_{m=1}^{M} \boldsymbol{Y}_{m} . $

(18)

其中：Y表示张量样本集Y的均值。张量 Y_m的n阶总散度矩阵可以定义为

$ \begin{equation*} \boldsymbol{S}_{T A}^{(n)}=\sum\nolimits_{m=1}^{M}\left(\boldsymbol{Y}_{m}^{(n)}-\overline{\boldsymbol{Y}}^{(n)}\right)\left(\boldsymbol{Y}_{m}^{(n)}-\overline{\boldsymbol{Y}}^{(n)}\right)^{\mathrm{T}} . \end{equation*} $

(19)

其中：Y_m⁽ⁿ⁾表示张量样本 Y_m的n阶展开矩阵，Y⁽ⁿ⁾为所有张量样本n阶展开矩阵的平均值。MPCA算法意在从原始张量中捕捉观察到的大部分变化信息，这种变化信息往往是分类的核心，而散度可以视为这种变化信息的度量，因此MPCA的核心问题就是如何找到能够令张量Y总散度最大的投影矩阵

$ \begin{equation*} \boldsymbol{U}_{n=1, 2, \cdots, N}^{(n)}=\underset{U^{(1)}, \cdots, U^{(N)}}{\operatorname{argmax}} \varphi_{Y} . \end{equation*} $

(20)

投影矩阵 U_{n=1, 2, …, N}⁽ⁿ⁾对原始张量的投影由N个子空间的投影组成，式(20)的寻优过程难以对所有子空间同时展开。当给定投影矩阵 U⁽¹⁾, …, U^(n-1), U⁽ⁿ⁺¹⁾, …, U^(N)时，矩阵 U⁽ⁿ⁾能够令张量集Y的n阶散度最大化，因此可以采用逐阶求解的方式依次求解投影矩阵。式(18)的n阶展开可以表示为

$ \begin{align*} \boldsymbol{\varphi}_{\boldsymbol{Y}_{(n)}}= & \sum\limits_{m=1}^{M}\left\|\boldsymbol{Y}_{m(n)}-\overline{\boldsymbol{Y}}_{(n)}\right\|_{\mathrm{F}}^{2} \\ = & \sum\limits_{m=1}^{M}\left\|\boldsymbol{U}^{(n) \mathrm{T}} \cdot\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right) \cdot \boldsymbol{U}_{\varPhi^{(n)}}\right\|_{F}^{2} \\ = & \sum\limits_{m=1}^{M} \operatorname{trace}\left(\boldsymbol{U}^{(n) \mathrm{T}} \cdot\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right) \cdot \boldsymbol{U}_{\varPhi^{(i)}} \cdot \boldsymbol{U}_{\varPhi^{(i)}}^{\mathrm{T}} \cdot\right. \\ & \left.\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right)^{\mathrm{T}} \cdot \boldsymbol{U}^{(n)}\right) \\ = & \operatorname{trace}\left(\boldsymbol{U}^{(n) \mathrm{T}} \cdot \boldsymbol{\varPhi}^{(n)} \cdot \boldsymbol{U}^{(n)}\right) . \end{align*} $

(21)

其中：X_(n)表示张量样本集X的均值X的n阶展开矩阵，

$ \begin{aligned} & \boldsymbol{U}_{\varPhi^{(n)}}= \\ & \left(\boldsymbol{U}^{(n+1)} \otimes \cdots \otimes \boldsymbol{U}^{(N)} \otimes \boldsymbol{U}^{(1)} \otimes \cdots \otimes \boldsymbol{U}^{(n-1)}\right)^{\mathrm{T}}, \end{aligned} $

(22)

$ \begin{align*} & \boldsymbol{\varPhi}^{(n)}= \\ & \sum\nolimits_{m=1}^{M}\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right) \cdot \boldsymbol{U}_{\varPhi^{(n)}} \cdot \boldsymbol{U}_{\varPhi^{(n)}}^{\mathrm{T}} \cdot\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right)^{\mathrm{T}} . \end{align*} $

(23)

当U⁽ⁿ⁾由矩阵Φ⁽ⁿ⁾的前L_n个最大特征值所对应的特征向量所构成时，式(21)取得最大值。

MPCA算法通过一个用户自定义参数Q决定各阶投影矩阵U⁽ⁿ⁾的维度L_n。对于每一阶张量展开矩阵，降维后的特征维数与原始维数之间满足关系

$ L_n=\operatorname{argmax}\left(\frac{\prod_{n=1}^N L_n}{\prod_{n=1}^N I_n} \leqslant Q\right) . $

(24)

Q的取值范围为(0, 1]，可以视为度量降维张量保留原始张量核心信息量的指标。Q=1表示保留全体信息的投影，有L_n=I_n，此时进行的投影称为全投影; 随着Q值的减小，L_n的值也必须逐渐减小以满足式(24)的限制。Q的值越接近于0，表示保留信息的比率越小，L_n与I_n之间的差距越大，降维张量 Y的维数也就越低。在算法的求解过程中，L_n的初值由矩阵 Φ₀⁽ⁿ⁾的累积概率决定

$ \begin{equation*} \boldsymbol{\varPhi}_{0}^{(n)}=\sum\nolimits_{m=1}^{M}\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right) \cdot\left(\boldsymbol{X}_{m(n)}-\overline{\boldsymbol{X}}_{(n)}\right)^{\mathrm{T}} . \end{equation*} $

(25)

MPCA算法的迭代过程可以总结如下：

1) 求解原始张量样本集均值X，并将X的n阶展开式代入式(24)，计算 Φ₀⁽ⁿ⁾(n=1, 2, …, N)；

2) 将 Φ₀⁽ⁿ⁾的特征值λ_t(1≤t≤I_n)依照从大到小的顺序排列，计算前L_n个特征值对应的特征向量，组成投影矩阵的初值U₀⁽ⁿ⁾；

3) 在第k次迭代中，将U_k-1⁽ⁿ⁾代入式(18)和式(23)，计算φ^k_Y(n)和Φ_k⁽ⁿ⁾，并求解Φ_k⁽ⁿ⁾的前L_n个特征值所对应的特征向量，组成新的投影矩阵U_k⁽ⁿ⁾；

4) 当投影张量的散度收敛(φ_Y^k-φ_Y^k-1 < η)，或当程序达到最大迭代次数时，停止迭代。

1.4 降维特征矢量化及决策树分类

MPCA方法可以从原始目标样本集中获取时间维度和极化维度的综合信息，随后对降维张量进行矢量化，即可通过机器学习分类器进行目标样本集的分类。本文采用回归分类树构造决策树分类器^[40]，该方法以基尼系数作为节点分裂的判定指标。基尼系数越小表示节点纯度越大，节点分裂对目标的分类效果越好。当基尼系数等于0时，表示某节点下的所有数据均属于同一类别。每次选择分裂属性时将遍历所有属性，寻找各属性分裂后基尼系数最小的特征分割界限，再对所有属性进行对比，选择一个最优分裂属性进行节点的二分，直到树的深度或者节点纯度达到阈值，或者所有的属性都被遍历。

2 实验区与实验数据

本文利用一组4景RADARSAT-2 Fine Quad模式SAR图像进行算法验证。实验区位于天津市武清区(中心坐标约为117°02′02″E，39°28′29″N)，农作物类型包含夏播玉米、水稻、荷花、草地、大豆等。实验数据的基本信息如表 1所示。图 2展示实验区的地理位置以及SAR图像覆盖区域在Google Earth中的光学图像。图 3(a)展示获取于7月26日的全极化数据的Pauli分解假彩色合成图。图 3(b)是根据野外实验考察以及对光学数据对比获取的样本分布图，该区域的主要农作物类型为玉米和大豆，其中玉米主要分布在实验区中部和南部，大豆主要分布在实验区中部和北部。

以7月26日获取的数据作为参考图像进行多时相极化SAR图像配准，并进行5×5的均值滤波，提取各景影像的协方差矩阵，随后对4景全极化图像的协方差矩阵进行联合超像素分割。超像素分割的效果受SAR图像分辨率、斑噪水平、地物类型及其分布等多方面因素影响。令式(7)中的权重参数m=1，图 4(a)~4(d)分别展示搜索半径R取值不同时的超像素分割局部结果。当搜索半径较小时，超像素目标包含像素较少，目标对像素空间同质特性的提取能力相对较低；此外，由于目标面积较小，像素相关性较强，内部梯度变化较小，因此在地物边界区域有时会出现过分割的现象，如图 4(a)和4(b)中黄色方框所示。随着搜索半径的增大，超像素目标面积逐渐增大，当搜索半径较大时，超像素目标中容易包含地物边界区域，导致分割结果内部的地物类型分布不均，如图 4(d)中蓝色方框所示。此外，当搜索半径较大时，分割得到的目标个数较少，难以产生足够的训练样本。因此，对实验数据进行超像素分割时设置搜索半径R=7，此时超像素目标能较好地维持地物边界，超像素目标内部包含较多的同质像素，在提取像素空间同质性信息的同时，能够获取较为充足的训练样本。此时，图 3(b)所示的样本区域被分割为1 758个目标，各地物类型目标具体数量如表 2所示。当数据源或研究区发生变化时，可以结合作物类型、成像条件以及图像质量等多方面因素，调整R的取值，以获取较好的分割效果，令目标地块兼具较高的内部均质性和较为准确的边界。

3 实验结果与分析

本文基于决策树分类器开展了多组农作物分类对比实验以检验本文方法的有效性。如表 3所示，实验1~实验4采取不同的数据表示方法和降维策略。一般地，机器学习分类方法将待分类样本表示为“样本数量×特征数量”的二维矩阵，即每行记录表示一个样本的特征向量。实验1和实验2均采用这一表示方式，直接提取协方差矩阵中9个独立元素(4个时相共计36维独立特征)作为分类特征。实验1直接使用决策树分类器进行农作物分类；实验2则先采用PCA方法对36维多时相协方差矩阵进行特征降维；实验3将Liu等^[41]提出的多频极化SAR张量表示方法应用于多时相极化SAR，随后利用MPCA方法进行特征降维；实验4采用本文方法实现多时相极化SAR的特征降维和农作物分类。实验3和实验4中，MPCA方法降维参数Q取值均设置为0.95，此时实验3和实验4的目标样本特征维数从36维(3×3×4)下降至4维(2×2×1)，因此，实验2中PCA特征降维维数也设置为4。

所有实验均从目标样本集中随机选取20 % 作为训练样本，并进行10次重复实验获取平均分类精度。其中，由于水稻和水体的目标数量较少，因此随机选择其50 % 作为训练样本。最终，实验1~实验4的多时相全极化SAR农作物分类精度对比如表 4所示，其中PA表示生产者精度，UA表示用户精度，OA表示总体精度，KA表示Kappa系数。

由表 4可以看出，直接采用决策树分类器对所有特征进行农作物分类的OA略低于本文方法。其中，大豆、水稻、水体的PA和UA，以及荷花的UA均显著低于本文方法。此外，根据对重复实验分类精度的逐一对比，发现实验1中水稻、水体等小样本地物的分类精度波动较大，水稻的最优PA和UA为96.97 % 和75.76 %，最低值仅为63.64 % 和42.47 %，而本文方法中水稻的最优PA和UA分别达到100.00 % 和91.18 %，最低值分别为81.82 % 和62.26 %，优于直接采用决策树进行分类的结果。类似地，实验1中水体的最低PA和UA仅为86.11 % 和90.00 %，明显低于本文方法。实验2采用PCA方法，与本文方法具有相同的降维特征数量，但在实验2的分类结果中，大豆、水稻、荷花样本发生了一定程度的混分，导致大豆、水稻、荷花的PA、UA和OA均低于本文方法。

由于极化协方差矩阵是3×3埃尔米特半正定矩阵，属于黎曼流形，可被视为2阶张量^[39]。因此，Tao等^[42]将张量空间引入全极化SAR分类，首先提取48维极化特征构成阶数为“行×列×48”的张量，再在张量空间中进行特征降维和分类。这种方式虽然可以直观地表示极化SAR的特征，但每个像素的特征空间依然是一维矢量，对极化张量的第3阶进行独立成分分析以降低特征维数的操作与对二维图像矢量化后进行特征降维并无本质区别。另一方面，Liu等^[41]对多频率极化SAR分类时将极化SAR的每个像素视为一个元胞，将元胞中心像素及其N×N邻域 T₃矩阵的N²×9个独立实数元素表示为特征矢量。实验3采用Liu等提出的张量表示方法，与本文方法进行对比，这2种方法的差异在于本文方法在特征降维阶段不进行实、虚部拆分，而是保持协方差矩阵本身的复数结构。由表 4可知，本文方法的总体精度优于实验3，原因在于实验3的分类结果中存在较多被误分为水稻和荷花的建筑样本，大豆与荷花之间的混淆也较为严重，而本文方法大豆、水稻、荷花、建筑的PA和UA均明显较优。

顺应多时相极化SAR的数据逻辑，保持协方差矩阵内部的数据结构，是本文方法的主要特点。相较于对多时相极化协方差矩阵独立元素直接进行分类、对多时相极化协方差矩阵PCA降维特征进行分类，以及对Liu等^[41]提出的极化SAR张量MPCA降维特征进行分类的方法，本文提出的方法取得了最高的总体分类精度和Kappa系数，对水稻、水体、荷花等样本数较少的类型，识别效果更优，表明本文方法能够从多时相极化SAR中提取有利于农作物精确分类的核心信息。

本文方法直接以极化SAR协方差矩阵作为输入，因而可以很容易地推广到简缩极化SAR数据。实验5和实验6分别将其应用于CTLR模式和π/4模式的多时相协方差矩阵。为便于实验对比，采用与全极化SAR相同的目标样本集，表 5给出了全极化与简缩极化SAR的分类精度对比。

由表 5可知，CTLR模式对大部分地物能取得与全极化SAR相当的分类精度，荷花和草地的UA以及OA甚至优于全极化SAR，原因是CTLR模式的分类结果中被误分为荷花的建筑样本和被误分为草地的玉米样本均相对较少，但水稻与建筑的误分现象较为严重，造成水稻的PA和UA相对较低。π/4模式的OA则略低于全极化SAR和CTLR模式，其中，草地UA略优于全极化SAR，其他精度指标大多与全极化SAR相当，但由于部分大豆样本被误分为荷花，π/4模式大豆的PA和荷花UA略低于全极化SAR。此外，与CTLR模式类似，π/4模式的分类结果中也存在水稻与建筑之间的混淆，导致较低的水稻和建筑分类精度。

图 5展示了全极化SAR和2种模式的简缩极化SAR分类结果，图 6展示了部分农作物样本感兴趣区域(region of interest, ROI)的分类结果。

2种模式简缩极化SAR的分类结果与全极化SAR基本相符，表明多时相简缩极化SAR基本可以取得与多时相全极化SAR相当的分类精度。由于本文方法采用了目标级分类策略，减少了结果中的噪声，这可能是使得简缩极化SAR的分类精度与全极化SAR相当、甚至更优的原因。在本文的实验数据中，由于水稻和荷花的样本数量较小，因此与其他类别发生混分时其精度更易受到影响。

全极化SAR和简缩极化SAR对玉米和大豆的分类能力较为相近；而全极化SAR对水稻和荷花的分类效果相较于简缩极化SAR具有一定的优越性，与精度统计结果相符。可能的原因是：全极化SAR包含最全面的农作物极化散射信息，即使水稻和荷花的样本数量较少，全极化SAR依然能一定程度上保持对其识别最关键的信息，从而取得优于简缩极化SAR的分类精度。

4 总结

本文提出一种基于张量表示的多时相极化SAR农作物分类方法，可同时适用于全极化和简缩极化SAR。该方法采用目标级分类策略，将SLIC方法引入多时相极化SAR实现农作物同质像素的超像素分割；将目标的多时相极化协方差矩阵表示为张量的形式，引入MPCA算法进行多时相极化特征降维；最后用决策树分类器实现农作物分类。对覆盖天津市武清区的4景RADARSAT-2全极化图像进行农作物分类实验，结果表明，相较于直接使用多时相极化协方差矩阵独立元素、使用多时相极化协方差矩阵PCA降维特征、基于其他极化SAR张量表示方法^[41]的MPCA降维特征作为决策树分类器输入特征的情况，本文方法能够取得最优的总体分类结果。此外，将本文方法应用于CTLR模式和π/4模式的简缩极化SAR，分类结果表明简缩极化SAR能取得与全极化SAR相当甚至更优的总体分类精度，但全极化SAR在水稻、荷花等小样本地物中具有更优的表现。

本文方法可以保持多时相极化SAR的时间相干性和极化二阶统计量的数据结构，通过目标分割和张量空间MPCA降维，实现对空间信息、极化散射信息和时序变化信息的综合利用，形成可同时适用于全极化和简缩极化SAR的多时相极化SAR农作物分类统一框架。未来，我们将面向复杂农作物精细分类的需求，构建时间序列数据集，充分考虑时间序列极化SAR多样的图像处理方法和丰富的极化特征，持续开展多模式时间序列极化SAR张量分类方法优化研究。