如何基于单极化SAR强度图构建距离测度进行变化检测

Astrophel2024-10-172024-10-25

概述

由于其不易受天气影响的特性，极化SAR在越来越多的对地观测场景下发挥着重大的作用，尤其是洪涝灾害的变化检测。洪涝灾害期间，常伴随着大雨或阴天，传统的光学影像难以发挥作用，此时仅能依靠SAR影像实现变化检测。
与服从简单正态分布的光学影像不同，SAR影像的概率分布比较复杂，原因是其涉及到各类平方及乘除法运算。已有的SAR变化检测方法，无论是多极化还是单极化，除先分类后检测的分类方法外，基本都依赖极化SAR影像的概率分布构建距离测度。距离测度最早应用于SAR影像的分类中，后许多学者对其加以改良，将此类方法应用于变化检测。这是很自然的，因为变化检测实质上就是一种最简单的二分类问题。通过理解变化检测中距离测度的构建，可以更容易理解分类场景下距离测度的构建。
本着有难到易，“先学走，后学跑”的原则，本文将介绍单极化SAR影像距离测度的构建，选取了同一区域不同时相的两景Sentinel-1影像作为例子进行演示。下文的第一部分，将介绍一些需要了解的预备知识；第二部分将介绍如何构建一种适用于单极化强度影像的极化距离测度；第三部分将使用构建的测度进行一次变化检测的尝试，定性地验证其效果。

预备知识

下面介绍阅读本文所需的一些预备知识，包括几种常见的概率分布，以及假设检验的一些知识。这里只进行基本的介绍，具体内容需要自行参阅其他材料。

正态分布

正态分布（Normal distribution），又称高斯分布（Gaussian distribution），是最重要的分布之一，其对整个数学和科学领域都非常重要。
如果随机变量$X$的概率密度函数是：

$f_X(x)={1\over \sqrt{2\pi \sigma ^2}}e^{-(x-\mu)^2/2\sigma ^2}$

那么$X$就服从均值为$\mu$，方差为$\sigma ^2$的正态分布，并记做$X\sim N(\mu,\sigma ^2)$，其形状是一个钟形曲线。

正态分布的特殊性在于它是一个“稳定分布”，即服从正态分布的相互独立的两个随机变量之和依然服从正态分布。

卡方分布

既然服从正态分布的相互独立的两个随机变量之和依然服从正态分布，那么服从正态分布的随机变量的平方服从何种分布呢？
答案是：标准正态分布的平方服从自由度为1的卡方分布，即：
如果$X\sim N(0,1)$，那么$X^2\sim \chi ^2(1)$.
由于任意正态分布都可以从原始数据中减去均值并除以标准差得到标准正态分布，因此该定理的应用范围非常广泛。
卡方分布定义为：
如果随机变量$X$服从自由度为$\nu > 0$的卡方分布，那么$X$的概率密度函数为：

$f(x)= \begin{cases} {1\over {2^{\nu/2}\Gamma(\nu/2)}}x^{(\nu/2-1)e^{-x/2}},x\geq 0 \\ 0,Others \end{cases}$

记做$X\sim \chi^2(\nu)$，$\nu$为自由度.
与正态分布相同，卡方分布也是一种“稳定的”分布，意味着一组独立随机卡方分布的和依然服从卡方分布，且新卡方分布的自由度等于一组卡方分布自由度的和。

指数分布

如果随机变量$X$满足:

$f(x)= \begin{cases} {1\over \lambda}e^{-x/\lambda},x\geq 0 \\ 0,Others \end{cases}$

那么说$X$服从参数为$\lambda > 0$的指数分布，并记做$X\sim Exp(\lambda)$.

Gamma分布

如果随机变量$X$满足:

$f(x)= \begin{cases} {1\over \Gamma(k)\sigma^k}x^{k-1}e^{-x/\sigma},x\geq 0 \\ 0,Others \end{cases}$

那么说$X$服从形状参数为$k$，尺度参数为$\sigma$的Gamma分布，并记做$X\sim Gamma(k,\sigma)$.
Gamma分布与指数分布间有一个有趣的定理：如果若干个随机变量是独立的，并且各自服从指数分布，那么其和服从伽马分布。

F-分布

在遥感影像处理中，常用到两个波段的比。那么，两个卡方分布的和依然为卡方分布，两个卡方分布的比是什么分布呢？答案是F-分布。
F分布是1924年英国统计学家Ronald.A.Fisher爵士提出，并以其姓氏的第一个字母命名的。它是两个服从卡方分布的独立随机变量各除以其自由度后的比值的抽样分布。如果两个卡方分布的自由度相等且为m，那么其比值就遵循自由度为m,m的F-分布。
更一般地，两个独立Gamma分布的比也服从F-分布。

中心极限定理

通俗的说，中心极限定理(CLT)是指，对于一些相互独立的“好”的随机变量，伴随着变量数量的不断增加，它们的和将收敛于高斯分布，而高斯分布的均值和方差由这些独立变量的均值和方差决定。
准确定义如下：设$X_1,…,X_N$是独立同分布的随机变量，它们的前3阶矩都是有限的，而且概率密度函数的表达式衰减的足够快。此时，令：

$\bar {X_N}={X_1 + ... + X_N\over N}$

并设：

$Z_N={\bar{X_N} - \mu \over \sigma / \sqrt N}$

那么，当N趋近于无穷大时，$Z_N$的分布将收敛于标准正态分布。

似然比及似然比检验

似然性是与概率相对应的一个概念。概率用于在已知一些参数的情况下，预测接下来的观测所得到的结果，而似然性则是用于在已知某些观测所得到的结果时，对有关事物的性质的参数进行估计。
似然函数是观测值$x$与其概率分布参数θ的函数，其在数值上等于等于每一个$x$取值的概率乘积。最大似然就是其实就是求似然函数最大时对应的θ值，方法就是对似然函数求导，找到其导数为0时的点。
与其他假设检验方法类似，似然比检验也要先构建两个假设，称为原假设和备择假设，也可以称为0假设和1假设。用于决定是否接受原假设的统计量称为检验统计量，通常用字母Q表示。顾名思义，在似然比检验中，似然比就是Q，或者说是原假设和备择假设最大似然值之比。
得到似然比后，我们可以通过似然比所服从的的概率密度函数，在某一个显著性水平$\alpha$下找到一个阈值，决定是否拒绝原假设。
需要注意的是，似然比检验并不是适用于本文场景的唯一一种检验方式，其他方法，例如KL散度等，也可以用于两种概率分布情况的相似度检验。

原理推导

在推导距离测度前，让我们回顾一下极化SAR影像的成像机理：

SAR影像的成像机理

在时空中传播的电磁波能够在到达一个特定目标后与之相互作用。在该过程中，一部分入射波能量被目标本身吸收，而其余部分则被再辐射形成电磁波。由于波与目标相互作用，可能会使再辐射波的性质与入射波不同。为了描述这种效应，引入散射矩阵$S$的概念。
将入射电磁波和散射电磁波分别表示为$E_1$和$E_s$（事实上是其琼斯矢量），包含观测对象影响的散射过程可由下式表达：

$E_S = {e^{jkr}\over r}SE_1 = {e^{jkr}\over r}\begin{pmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \\ \end{pmatrix}E_1$

式中，$S_{ij}$表示复散射系数或复散射幅度，对角线元素表示相同极化方式的入射场和散射场关系，即同极化项；非对角线元素代表正交极化方式的入射场和散射场关系，称为交叉极化项，${e^{jkr}\over r}$表示电磁波传播本身所引起的幅度和相位变化。
对于全极化SAR，一个天线交替发射水平极化和垂直极化的电磁波，两个天线同时接收水平极化和垂直极化回波。对于双极化SAR，仅有一个发射天线。因此，对于Sentinel-1 VV+VH双极化数据，以上公式可以被描述为：

$\begin{pmatrix} E_v^b \\ E_h^b \\ \end{pmatrix}= {e^{jkr}\over r} S \begin{pmatrix} E_v^i \\ 0 \\ \end{pmatrix} ={e^{jkr}\over r} \begin{pmatrix} S_{vv} & S_{vh} \\ S_{hv} & S_{hh} \\ \end{pmatrix} \begin{pmatrix} E_v^i \\ 0 \\ \end{pmatrix}$

其中，S矩阵可以通过传感器测量获得其中的两项，即：

$S=\begin{pmatrix} S_{vv} \\ S_{vh} \\ \end{pmatrix}$

这也是我们下载的SLC影像经辐射定标等预处理过程所获得的值。
对于中低分辨率SAR下森林、农田、裸地等“分布式目标”，引入极化相干矩阵$T$与极化协方差矩阵$C$来进行描述。事实上，在双极化情况下，二者是等价的。C被定义为：

$C_2=SS^H=\begin{pmatrix} S_{vv} \\ S_{vh} \\ \end{pmatrix}\begin{pmatrix} S_{vv}^*~ S_{vh}^* \\ \end{pmatrix}= \begin{pmatrix} |S_{vv}|^2 & S_{vv}^*S_{vh} \\ S_{hv}^*S_{vv} & |S_{hh}|^2 \\ \end{pmatrix}$

式中，H表示共轭转置，$*$表示共轭。
对于一个极化通道的发射和接收，（如vv），接收（测量）到的 SLC 信号可模拟为以下形式：

$S_{vv}={|S_{vv}^a|\over \sqrt n}\Sigma_{k=1}^ne^{i\phi_k}$

式中，$|S_{vv}^a|$为单个像素覆盖区域散射信号的总振幅（下文也可能使用“基本信号振幅”表示），n为随机散射体数量，求和项说明了随机散射体相干效应引起的接收信号相位变化。这个式子说明了SAR的散射系数（测量到的信号）并不等于单个像素覆盖区域散射信号的总振幅（基本信号振幅），而是还与地面散射体响应有关。这也是SAR影像斑点的来源。

概率分布推导

在上面的工作中，我们回顾了SAR的成像机理和对其的常见数学描述。下面，我们将推导SAR传感器接收到的信号服从何种概率分布，进而对其实施似然比检验，度量两景影像的相似性。
将测量到的SLC信号模拟式指数项展开为三角函数形式：

$S_{vv}={|S_{vv}^a|\over \sqrt n}\Sigma_{k=1}^ne^{i\phi_k}={|S_{vv}^a|\over \sqrt n}(x+iy)$

其中：

$x=\Sigma_k cos(\phi_k),y=\Sigma_k sin(\phi_k)$

由于相移$\phi k$是随机均匀分布的，因此变量x和y分别是同分布余弦项和正弦项的和。根据统计学的中心极限定理，x和y在大量散射体n的极限情况下，将服从均值为0且方差为$\sigma ^2 = n/2$的正态分布。
我们日常看到的SAR影像通常是强度图，它有很多种计算方法且都是等价的，在这里，我们使用$|S_{vv}|^2$表示其强度，这也是极化协方差矩阵C的其中一个主对角线元素。此时有：

$|S_{vv}|^2=S_{vv}S_{vv}^*={|S_{vv}|^2 \over n}(x^2+y^2)$

将$x,y$化为标准正态分布以便后续推导，即除以其方差$n/2$：

$|S_{vv}^2|={|S_{vv}|^2 \over n}{n\over 2}({x^2\over n/2}+{y^2\over n/2})={|S_{vv}|^2 }{u\over 2}$

式中，

$u={x^2\over n/2}+{y^2\over n/2}$

u是服从标准正态分布的变量的平方和，根据正态分布与卡方分布的关系，u是自由度m=2的卡方分布。
为简化符号，令:

$s = |S_{vv}|^2,a = {\lvert S_{vv}^a\rvert}^2$

那么：

$s=a{u\over 2}$

为了从u服从的卡方分布推导观测信号强度s的分布，应用标准变化公式（standard transformation formula）：

$p_s(s)=p_u(u)|{du \over ds}|={1\over 2}e^{-u/2}{2\over a}={1\over a}e^{-s/a}$

该式与指数分布的定义式相符。因此，最终得到s服从指数分布，其均值和方差等于基本信号强度a。
考虑到通常情况下，极化SAR影像需经过多视平均处理，根据指数分布与Gamma分布之间的关系，多个指数分布的和服从Gamma分布，因此$m-look$多视后的SAR影像服从Gamma分布，其参数为参数为$\alpha = m,\beta = a/m$。多视影像的均值依然等于a.

似然比检验

现在可以设置原假设和备择假设：我们的目的是度量两景影像是否相似，进而找出哪里发生了变化——原假设通常与我们的目的相反，因此原假设$H_0$为假设在两景影像中，信号的功率没有发生变化，那么备择假设$H_1$就是两景影像中，信号强度发生了变化，即:

$\begin{equation} \begin{cases} H_0: a_1 = a_2 = a \\ H_1: a_1 \neq a_2 \end{cases} \end{equation}$

式中，用$a$指代单个像素覆盖区域散射信号强度（基本信号强度）。
由于上文已经推导出，多视平均后的强度影像服从Gamma分布，参数为$\alpha = m,\beta = a/m$，其概率密度为：

$p(s_i,a_i)={1\over {(a_i/m)^m}\Gamma(m)}s_i^{m-1}e^{-s_im/a_i},i=1,2$

下面计算似然函数。
当原假设为真时：由于两时相单极化通道可以被视为独立的，所以其似然函数就是其概率密度的乘积：

$L_0(a)=p(s_1,a_1)p(s_2,a_2)=p(s_1,a)p(s_2,a)$

设其最大似然为：

$L_0(\hat a)=p(s_1,\hat a)p(s_2,\hat a),\hat a=arg max L_0(a)$

令$L_0$导数为0，求得其极大似然估计值：

$\hat a = {s_1+s_2\over 2}$ $L_0(\hat a)={1\over ((s_1+s_2)/2m)^{2m}\Gamma(m)^2}(s_1s_2)^{m-1}e^{-2m}$

同理可得，在备择假设成立时：

$\hat a_1 = s_1,\hat a_2 = s_2$ $L_1(\hat a_1,\hat a_2)={m^{2m}\over\Gamma(m)^2}s_1s_1e^{-2m}$

计算似然比Q，若：

$Q={L_0(\hat a)\over L_1(\hat a_1,\hat a_2)}=[2^2({s_1s_2\over (s_1+s_2)^2})]^m\leq k$

成立，则拒绝原假设，转而接受备择假设——即认为图像发生了变化。k为某个适当的阈值。
经过一系列代数运算，可以得到，当：

$Q_1={s_1\over s_2}<c_1~or~Q_2={s_2\over s_1}<c_2$

上式成立。$c_1,c_2$均是与k相关的一个常数。现在只要确定常数的取值，就能找出变化区域。
其实，Q可以被姑且认为是一种距离测度。但是距离测度还要满足非负性、对称性和最小恒定性，Q是否是一种距离测度还需要进一步验证。
那么如何确定常数的取值呢？
$s_1,s_2$服从Gamma分布，令$z=2sm/a$，容易证明z服从自由度为2m的卡方分布。因此变量$2s_1m/a$和$2s_2m/a$服从卡方分布。由于卡方分布之比服从F分布，对其求比值，得到$s1/s2$和$s2/s1$服从自由度为2m，2m的F分布。
现在就可以对Q（或者说分别对$Q_1,Q_2$）进行检验了。设显著性系数$\alpha=0.05$，由于这是一个双侧检验，我们利用F分布表来找到检验统计量远离0的概率，如果概率小于$\alpha /2$，则拒绝原假设转而支持备择假设。（事实上，临界值就是概率值等于$\alpha /2$时自由度为（2m，2m）的F分布的ppf值）

实验验证

选取Sentinel-1鄱阳湖流域两景影像实施基于似然比检验的单极化通道强度变化检测：

使用F-分布对比其单极化通道时相比直方图：

下面，我们以VH极化通道的变化检测为例（VH结果相对较好，适合举例）：使用Matlab分别读取两张影像的VH极化通道，然后计算其似然比$Q_1,Q_2$，最后使用F分布概率密度函数，以0.001为显著性系数进行分割，得到变化区域。具体代码如下：

% 数据读取，其中read_ENVI函数为自建函数，实质是读取二进制数据，可以使用Matlab自带函数替代
file1 = '/Volumes/Astrophel/DATA/PolSAR/Sentinel-1/poyangBIG/20210428_LEE/C2/C11.bin';
file2 = '/Volumes/Astrophel/DATA/PolSAR/Sentinel-1/poyangBIG/20210428_LEE/C2/C22.bin';
file3 = '/Volumes/Astrophel/DATA/PolSAR/Sentinel-1/poyangBIG/20210522_LEE/C2/C11.bin';
file4 = '/Volumes/Astrophel/DATA/PolSAR/Sentinel-1/poyangBIG/20210522_LEE/C2/C22.bin';
vv1 = NormalizeCommon(OutlierElimination(read_ENVI(file2)','maxmin'),0,1);
vh1 = NormalizeCommon(OutlierElimination(read_ENVI(file1)','maxmin'),0,1);
vv2 = NormalizeCommon(OutlierElimination(read_ENVI(file4)','maxmin'),0,1);
vh2 = NormalizeCommon(OutlierElimination(read_ENVI(file3)','maxmin'),0,1);
% 似然比Q的计算
Q1vh = vh1./vh2;
Q2vh = vh2./vh1;
% F分布PPF计算
x = finv(0.0005,8,8);
% 两个变化方向的提取
Q1vh1 = Q1vh;
Q1vh1(Q1vh < 10*x) = 1;
Q1vh1(Q1vh >= 10*x) = 0;

Q2vh1 = Q2vh;
Q2vh1(Q2vh < 10*x) = 1;
Q2vh1(Q2vh >= 10*x) = 0;