Wishart距离简介

深度学习自诞生之后，就被各个领域视为全能方法。而传统机器学习因为要手动提取特征，因此对数据本身的研究要远远强于深度学习方法，因此若可以将特定数据集的特征应用到深度学习中，想必是一个很好的方向。

对于极化SAR数据而言，Wishart距离是一个不可避免的概念，在论文[2]中详细介绍了Wishart距离的由来，下面结合这个论文对Wishart距离进行总结。

Single-look极化SAR图像分类器

在给定入射角时，一个介质的完全散射复矩阵$S$可以由下面式子给出：

$$\begin{eqnarray} {\rm S} = \left[ \begin{array}{c} S_{hh} & S_{hv} \\ S_{vh} & S_{vv} \end{array} \right] \end{eqnarray}$$

对于一个互易介质，满足$S_{hv}=S_{vh}$，复数向量可以使用三个不同的元素定义：

$$\begin{eqnarray} {\rm u} = \left[ \begin{array}{c} S_{hh} \\ S_{hv} \\ S_{vv} \end{array} \right] \end{eqnarray}$$

如果要求$\rm u$满足$|\rm u|=span(总功率)$，则$S_{hv}$的系数应该为$\sqrt 2$。但是该系数对于下面的分类结果没有任何影响。因此可以省略。

当雷达照射到一个由许多基本散射体组成的随机表面的区域时，$\rm u$满足多元复高斯分布：

$$\begin{eqnarray} p({\rm u})= \frac{1}{\pi ^3 |{\rm C}|}e^{-{\rm u}^T {\rm C}^{-1} {\rm u}} \end{eqnarray}$$

其中，${\rm C}=E\left[{\rm u} {\rm u}^T \right]$，$T$为共轭转置，$|\rm {C}|$为$\rm C$的行列式。复协方差矩阵$\rm C$为厄米矩阵（Hermitian，对角线为实数，其余元素共轭对称），即满足${\rm C}={\rm C}^T$。我们定义$\rm CN(0,C)$为均值为0，协方差为$\rm C$的复高斯分布。令${\rm u}_j=x_j+iy_j$为$\rm u$的一个元素，则$\rm u$满足复高斯分布的条件为$E[x_j]=E[y_j]=0$，$E[x_j^2]=E[y_j^2]$，$E[x_jy_j]=0$，$E[x_jx_k]=E[y_jy_k]$，$E[y_jx_k]=-E[x_jy_k]$。

定义$w_m$为第$m$类。每一个类均可以由自己的协方差矩阵表征，我们称之为特征协方差矩阵。第$m$类（$w_m$）的特征协方差矩阵${\rm C}_m$可由第$w_m$类的训练样本评估得到。也就是说在上面公式$({3})$中的$\rm C$理论上是要求${\rm u} {\rm u}^T $的期望，但是期望通常是求不出来的。所以通常使用第$w_m$类的训练样本评估得到第$w_m$类的${\rm C}_m$。在已知第$w_m$类的${\rm C}_m$后就可知道第$w_m$类$\rm u$的概率密度函数。

根据贝叶斯最大似然分类器，如果满足下式，则$\rm u$属于第$w_m$类：

$$\begin{eqnarray} P(w_m|{\rm u}) \geq P(w_j|{\rm u}), for \ all \ j \neq m \end{eqnarray}$$

运用贝叶斯公式，可以得到

$$\begin{eqnarray} P(w_m|{\rm u}) = \frac{p({\rm u}|w_m) P(w_m)}{p(\rm u)} \end{eqnarray}$$

进一步可以得到

$${\rm u} \in w_m, \ if \ p({\rm u}|w_m) P(w_m) > p({\rm u}|w_j) P(w_j), for \ all \ j \neq m$$

其中，$p({\rm u}|w_m)$满足${\rm CN(0,C}_m)$，$P(w_m)$为第$w_m$类的先验概率。

现在是一个最大化问题，而对于距离是一个最小化问题。对公式$({5})$等号右边的分子部分取负对数，可以得到：

$$\begin{eqnarray} \begin{split} -{\rm ln}p({\rm u}|w_m) P(w_m)&=&-{\ln}p({\rm u}|w_m) -{\ln}P(w_m) \\ &=&-{\rm ln} \frac{1}{\pi ^3 |{\rm C}_m|}e^{-{\rm u}^T {\rm C}_m^{-1} {\rm u}} -{\ln}P(w_m) \\ &=& {\rm u}^T{\rm C}_m^{-1}{\rm u} + {\rm ln}|{\rm C}_m| + {\rm ln \pi^3} -{\ln}P(w_m) \end{split} \end{eqnarray}$$

其中，${\rm ln \pi^3}$是常数，对于距离的衡量没有作用，故舍去，得到最终的距离度量公式：

$$\begin{eqnarray} d_1({\rm u}, w_m)= {\rm u}^T{\rm C}_m^{-1}{\rm u} + {\rm ln}|{\rm C}_m|-{\ln}P(w_m) \end{eqnarray}$$

若$ d_1({\rm u}, w_m) \leq d_1({\rm u}, w_j), \ for \ all \ j \neq m$，则${\rm u} \in w_m$。

Multi-look极化SAR图像分类器

为了斑点抑制和数据压缩，极化SAR数据通常采用multi-look处理。这种处理要求对几个独立的single-look的协方差矩阵进行平均，也就是：

$$\begin{eqnarray} Z=\frac{1}{n}\sum_{k=1}^n {\rm u}(k){\rm u}(k)^T \end{eqnarray}$$

其中，$n$为look的数量，向量${\rm u}(k)$是第$k$个single-look样本。

我们让：

$$\begin{eqnarray} {\rm A}=nZ=\sum_{k=1}^n {\rm u}(k){\rm u}(k)^T \end{eqnarray}$$

则矩阵$ {\rm A}$满足复Wishart分布。Goodman（1963）根据特征函数和傅里叶变换导出了复Wishart分布。Srivastava（1965）给出了它的简化版本，概率密度函数为：

$$\begin{eqnarray} p({\rm A})=\frac{|{\rm A}|^{n-q}e^{\left[-Tr(C^{-1}{\rm A})\right]}}{K(n,q)|C|^n}\end{eqnarray}$$

其中，$Tr(C^{-1}{\rm A})$为$C^{-1}{\rm A}$的迹，$q$为向量$\rm u$的维度，对于在互易介质的单极化SAR数据，$q=3$，并且

$$\begin{eqnarray} K(n,q)=\pi ^{\frac{1}{2}q(q-1)}\Gamma(n)...\Gamma(n-q+1) \end{eqnarray}$$

同上面Single-look极化SAR图像分类器一样，$C$理论上是要求${\rm u} {\rm u}^T $的期望，但是期望通常是求不出来的。所以通常使用第$w_m$类的训练样本评估得到第$w_m$类的$C_m$。在已知第$w_m$类的$C_m$后就可知道第$w_m$类$A$的概率密度函数。所以使用${\rm C}_m$代替$C$，我们可以将$({11})$式重写为$p({\rm A}|w_m)$。类似于公式$({6})$，我们需要最大化$p({\rm A}|w_m)p(w_m)$，对其取负自然对数可以得到：

$$\begin{eqnarray} \begin{split} d({\rm A}, w_m)&=&-{\rm ln}p({\rm A}|w_m) P(w_m) \\ &=&-{\ln}p({\rm A}|w_m) -{\ln}P(w_m) \\ &=&-{\rm ln} \frac{|{\rm A}|^{n-q}e^{\left[-Tr({\rm C}_m^{-1}{\rm A})\right]}}{K(n,q)|{\rm C}_m|^n} -{\ln}P(w_m) \\ &=& n {\rm ln}|{\rm C}_m| + Tr({\rm C}_m^{-1}{\rm A}) -{\ln}P(w_m) - (n-q) {\rm ln |A|}+{\rm ln}[K(n,q)] \end{split} \end{eqnarray}$$

因为后面两项不是关于$w_m$的函数，对于分类没有贡献，所以可以省略。带入${\rm A}=nZ$，用于$n$-look极化SAR图像分类的距离度量为：

$$\begin{eqnarray} d(Z, w_m)= n \left[ {\rm ln}|{\rm C}_m| + Tr({\rm C}_m^{-1}{\rm Z}) \right]-{\ln}P(w_m) \end{eqnarray}$$

从该式子可以看出，随着looks数量$n$的增加，先验概率$P(w_m)$重要性越低。对于不知道每个类别所占比例的极化SAR数据，假设$P(w_m)$相等，此时距离度量不是关于$n$的函数。为了实现这个度量，${\rm C}_m$从第$m$类的训练区域中评估得到，然后对其余像素点进行分类。

Wishart距离公式分析

乍一看，上面的距离度量公式$({14})$还是比较的复杂的。为了更好理解，下面对该公式的运算进行详细的分析。

因为是来自不同论文的分析，所以数学符号可能不太相同。具体而言，下面的$\rm \langle T \rangle$对应上面的$Z$，下面的$\rm \langle T_m \rangle$对应上面的${\rm C}_m^{-1}$，下面的$d\left(\langle \rm{T} \rangle | \langle \rm{T}_m \rangle \right) $对应上面的$d(Z, w_m)$。

从上面的介绍可知，极化SAR数据常用的相干矩阵以及协方差矩阵是具有复Wishart分布的（现在所用的数据大多是multi-look的？？）。而相干矩阵和协方差矩阵之间可以使用线性变换互相转换，所以这里只以复相干矩阵（Complex coherency matric $\rm \langle T \rangle$）为例，而它是共轭对称（conjugate symmetric）的，即：

$$\begin{eqnarray} \rm \langle T \rangle = \left[ \begin{array}{c} T_{11} & T_{12} & T_{13} \\ \overline{T_{12}} & T_{22} & T_{23} \\ \overline{T_{31}} & \overline{T_{32}} & T_{33} \end{array} \right] \end{eqnarray}$$

其中，$T_{11}, T_{22}, T_{33}$为实值的，剩下的元素是复值的，$\overline{\circ}$表示共轭。

根据上面Multi-look极化SAR图像分类器小节的分析，multi-look POLSAR像素$\rm \langle T \rangle$能够根据Wishart距离$d\left(\langle \rm{T} \rangle | \langle \rm{T}_m \rangle \right)$分类，正如下面公式所示：

$$\begin{eqnarray} d\left(\langle \rm{T} \rangle | \langle \rm{T}_m \rangle \right) = Trace \left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right) + ln| \langle \rm{T}_m \rangle| \end{eqnarray}$$

其中，$Trace(\cdot)$是矩阵的迹，$\cdot ^{-1}$是矩阵的逆，$|\cdot|$是矩阵的行列式，$ \langle \rm{T}_m \rangle$是由训练集的第$m$类评估得到，可以视为第$m$类的聚类中心。通过计算一个像素点与所有类聚类中心的距离，距离最近的类即为该像素点所属的类别。

对于无监督训练，可以先使用聚类的方法得到不同类别的数据；而对于有监督训练，训练集中已经包含了不同类别的训练数据。因此$\langle \rm{T}_m \rangle$可以使用如下公式得到：

$$\begin{eqnarray} \langle \rm{T}_m \rangle = E \left[ \langle \rm{T} \rangle | \langle \rm{T} \rangle \in \Omega_m \right ] = \frac{1}{|\Omega_m|} \sum_{ \langle \rm{T} \rangle \in \Omega_m} \langle \rm{T} \rangle \end{eqnarray}$$

其中，$\Omega_m$是第$m$类的像素集合，而$|\Omega_m|$是$\Omega_m$的像素总数。

为了下面表述方便，我们记$Trace \left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right) $为第一项，记做$ ln| \langle \rm{T}_m \rangle| $为第二项。

第一项分析

根据论文[1]中的说明，由于$\langle \rm{T}_m \rangle^{-1}$为共轭矩阵，我们定义其为如下形式：

$$\begin{eqnarray} \rm \langle T_m \rangle^{-1} = \left[ \begin{array}{c} a_{11} & a_{12}+jb_{12} & a_{13}+jb_{13} \\ a_{12}-jb_{12} & a_{22} & a_{23}+jb_{23} \\ a_{31}-jb_{13} & a_{32}-jb_{23} & a_{33} \end{array} \right] \end{eqnarray}$$

而$ \langle \rm{T} \rangle$也为共轭矩阵，我们定义其为如下形式：

$$\begin{eqnarray} \rm \langle T \rangle = \left[ \begin{array}{c} c_{11} & c_{12}+jd_{12} & c_{13}+jd_{13} \\ c_{12}-jd_{12} & c_{22} & c_{23}+jd_{23} \\ c_{31}-jd_{13} & c_{32}-jd_{23} & c_{33} \end{array} \right] \end{eqnarray}$$

则$\langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle $可以表示成下式，注意这里的$T$为斜体，即$\langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle=\it T$。

$$\begin{eqnarray} \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle = \left[ \begin{array}{c} T_{11} & & \\ & T_{22} & \\ & & T_{33} \end{array} \right] \end{eqnarray}$$

其中，

$$\begin{eqnarray} \begin{split} T_{11}&=&a_{11}c_{11}+(a_{12}+jb_{12})(c_{12}-jd_{12})+(a_{13}+jb_{13})(c_{13}-jd_{13})\\ T_{22}&=&a_{22}c_{22}+(a_{12}-jb_{12})(c_{12}+jd_{12})+(a_{23}+jb_{23})(c_{23}-jd_{23})\\ T_{33}&=&a_{33}c_{33}+(a_{13}-jb_{13})(c_{13}+jd_{13})+(a_{23}-jb_{23})(c_{23}+jd_{23}) \end{split} \end{eqnarray}$$

由该公式可得，$(a_{12}+jb_{12})(c_{12}-jd_{12})$与$(a_{12}-jb_{12})(c_{12}+jd_{12})$共轭，$(a_{13}+jb_{13})(c_{13}-jd_{13})$与$(a_{13}-jb_{13})(c_{13}+jd_{13})$共轭，$(a_{23}+jb_{23})(c_{23}-jd_{23})$与$(a_{23}-jb_{23})(c_{23}+jd_{23})$共轭。所以$\left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right)$的计算公式为：

$$\begin{eqnarray} \begin{split} Trace \left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right)&=&T_{11}+T_{22}+T_{33}\\ &=&a_{11}c_{11}+ a_{22}c_{22}+ a_{33}c_{33} \\ &+&(a_{12}+jb_{12})(c_{12}-jd_{12}) + (a_{12}-jb_{12})(c_{12}+jd_{12}) \\ &+&(a_{23}-jb_{23})(c_{23}+jd_{23}) + (a_{23}+jb_{23})(c_{23}-jd_{23}) \\ &+&(a_{13}-jb_{13})(c_{13}+jd_{13}) + (a_{13}+jb_{13})(c_{13}-jd_{13}) \end{split} \end{eqnarray}$$

可知两个共轭复数相加，则消去虚部，实部相加，可得最终的最简结果：

$$\begin{eqnarray} \begin{split} Trace \left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right) &=&a_{11}c_{11}+ a_{22}c_{22}+ a_{33}c_{33} \\ &+&2(a_{12}c_{12}+b_{12}d_{12}) \\ &+&2(a_{13}c_{13}+b_{13}d_{13}) \\ &+&2(a_{23}c_{23}+b_{23}d_{23}) \\ &=&a_{11}c_{11}+ a_{22}c_{22}+ a_{33}c_{33} \\ &+&2(a_{12}c_{12}+b_{12}d_{12} + a_{13}c_{13}+b_{13}d_{13} + a_{23}c_{23}+b_{23}d_{23}) \end{split} \end{eqnarray}$$

由上面公式$({23})$可得，Wishart距离中的$Trace \left( \langle \rm{T}_m \rangle^{-1} \langle \rm{T} \rangle \right)$为实数值。

第二项分析

因为$ \langle \rm{T}_m \rangle$为复共轭对称矩阵，我们记为：

$$\begin{eqnarray} \rm \langle T_m \rangle = \left[ \begin{array}{c} a_{11} & a_{12}+jb_{12} & a_{13}+jb_{13} \\ a_{12}-jb_{12} & a_{22} & a_{23}+jb_{23} \\ a_{31}-jb_{13} & a_{32}-jb_{23} & a_{33} \end{array} \right] \end{eqnarray}$$

计算$ \langle \rm{T}_m \rangle$的行列式，如下：

$$\begin{eqnarray} \begin{split} | \langle \rm{T}_m \rangle| &=& a_{11} \left[a_{22}a_{23}-(a_{23}+jb_{23})(a_{23}-jb_{23})\right] \\ &\ &- (a_{12}+jb_{12})\left[(a_{12}-jb_{12})a_{33} -(a_{23}+jb_{23})(a_{13}-jb_{13}) \right]\\ &\ & + (a_{13}+jb_{13})\left[(a_{12}-jb_{12})(a_{23}-jb_{23}) -a_{22}(a_{13}-jb_{13}) \right] \\ &=& a_{11}a_{22}a_{33}-a_{11}(a_{23}^2+b_{23}^2) \\ &\ & -(a_{12}^2+b_{12}^2)a_{33}+(a_{12}+jb_{12})(a_{23}+jb_{23})(a_{13}-jb_{13}) \\ &\ & +(a_{13}+jb_{13})(a_{12}-jb_{12})(a_{23}-jb_{23})-a_{22}(a_{13}^2+b_{13}^2) \end{split} \end{eqnarray}$$

上面式子中只有第四项和第五项有虚部，将这两项单独拉出来分析如下：

$$\begin{eqnarray} \begin{split} &=&(a_{12}+jb_{12})(a_{23}+jb_{23})(a_{13}-jb_{13})+(a_{13}+jb_{13})(a_{12}-jb_{12})(a_{23}-jb_{23}) \\&=&(a_{12}a_{13}-ja_{12}b_{13}+ja_{13}b_{12}+b_{12}b_{13})(a_{23}+jb_{23}) \\&\ &+ (a_{12}a_{13}+ja_{12}b_{13}-ja_{13}b_{12}+b_{12}b_{13})(a_{23}-jb_{23}) \\ &=&a_{12}a_{13}a_{23}+ja_{12}a_{13}b_{23}-ja_{12}b_{13}a_{23}+a_{12}b_{13}b_{23}+ja_{13}b_{12}a_{23}-a_{13}b_{12}b_{23}+b_{12}b_{13}a_{23}+jb_{12}b_{13}b_{23} \\ &\ &+a_{12}a_{13}a_{23}-ja_{12}a_{13}b_{23}+ja_{12}b_{13}a_{23}+a_{12}b_{13}b_{23}-ja_{13}b_{12}a_{23}-a_{13}b_{12}b_{23}+b_{12}b_{13}a_{23}-jb_{12}b_{13}b_{23} \\ &=&2a_{12}a_{13}a_{23}+2a_{12}b_{13}b_{23}-2a_{13}b_{12}b_{23}+2b_{12}b_{13}a_{23} \end{split} \end{eqnarray}$$

可以看到，最终两项的虚部相互抵消，只剩下了实部。综合上面两个式子，可以得到：

$$\begin{eqnarray} \begin{split} | \langle \rm{T}_m \rangle| &=& a_{11}a_{22}a_{33}-a_{11}(a_{23}^2+b_{23}^2) -a_{33}(a_{12}^2+b_{12}^2)-a_{22}(a_{13}^2+b_{13}^2) \\ &\ &+ 2a_{12}a_{13}a_{23}+2a_{12}b_{13}b_{23}-2a_{13}b_{12}b_{23}+2b_{12}b_{13}a_{23} \end{split} \end{eqnarray}$$

$ln| \langle \rm{T}_m \rangle|$是一个实数值。

结论

Wishart距离本身为实数值。

我们使用MATLAB对其验证这个结论。已知矩阵matrix是一个厄米矩阵，矩阵的形式如下：

矩阵a也是一个厄米矩阵，矩阵形式如下：

代码如下：

K>> format longg;
K>> det(matrix)

ans =

       0.000242440443638616 -  2.30476891505805e-23i

K>> log(det(matrix))

ans =

          -8.32475447099251 -  9.50653645269501e-20i
          
K>> trace(matrix\a)

ans =

          0.413019547470957 -  1.30104260698261e-18i
 
K>> trace(matrix\a) + log(det(matrix))

ans =

          -7.91173492352156 -  1.39610797150956e-18i

可以看到，矩阵a与矩阵matrix的Wishart矩阵基本上是只有实部有值，虚部的值可以忽略不记，这里之所以虚部有很小很小的值，是因为存在计算误差的原因。

有的时候，因为矩阵的行列式出现了负值，因此在对其取对数时，会出现虚部，甚至出现虚部前面的系数比实部更大的情况。

K>> value = -0.11248527642128878-1.474514954580286e-17j

value =

         -0.112485276421289 -  1.47451495458029e-17i

K>> log(value)

ans =

          -2.18493294215792 -      3.14159265358979i

K>> log(-0.11248527642128878)

ans =

          -2.18493294215792 +      3.14159265358979i

K>> exp(-2.18493294215792+3.14159265358979i)

ans =

         -0.112485276421288 +  3.63449956051925e-16i

K>> exp(-2.18493294215792+0i)

ans =

         0.112485276421288

出现了这种情况，一开始我一直在纠结Wishart距离为复数值怎么比较大小呢？为什么所有的代码都是只取其实部呢？为什么所有的论文以及参考书了里面都不提怎么处理复数值呢？自己手动推到了一遍之后才发现原来Wishart矩阵理论上就是实数，而之所以出现虚部，是因为矩阵行列式出现了负值的原因。

总结

Wishart距离有如下几个性质：

• 本身为实数
• 不满足齐次性、对称性与三角不等式，因此不算是严格意义上的距离。关于这个的讨论可以在[3]中找到
• Wishart距离有可能是负数，关于这个的讨论可以在Is Wishart distance always positive?找到

关于Wishart的更多信息可以参考下面的书籍[4]，该书籍的下载链接在这里

参考

[1] 基于散射能量和Wishart的深度学习极化SAR分类_刘永坤
[2] J. S. Lee, M. R. Grunes, and R. Kwok, “Classification of multi-look polarimetric SAR imagery based on complex Wishart distribution,” Int. J. Remote Sens., vol. 15, no. 11, pp. 2299–2311, 1994.
[3] Hänsch R, Jäger M, Hellwich O. Clustering by deterministic annealing andWishart based distance measures for fully-polarimetric SAR-data[C]//7th European Conference on Synthetic Aperture Radar. VDE, 2008: 1-4.
[4] Lee J S, Pottier E. Polarimetric radar imaging: from basics to applications[M]. CRC press, 2009.
Is Wishart distance always positive?
Wishart分布简介
维希特分布 (Wishart) 的分布密度