在进行极化SAR图像的解译和研究时，变化检测是非常重要的一个环节。经典的变化检测流程通常分为三个步骤：图像预处理、差异图提取和差异图分析。经过三个步骤，输入的两个时相多极化影像之间的差异部分将会被标记出来，以便进一步研究。

对于最近正在做的洪涝灾害评估，变化检测是标记洪涝前后淹没范围的最主要方法。由于网络上目前的资料鲜有代码，不方便理解，故把自己写的Matlab代码放上来。其中有理解错误的地方还请谅解。

与普通的SAR影像不同，极化SAR要考虑到多个极化波段的信息，故不能简单的对灰度图进行差值或比值运算——因此，在提取差异图的时候，需要对传统方法进行改进。目前对于极化SAR影像差异图提取的主流方法有以下几种：一是基于主成份分析的原理，二是基于极化距离的，三是基于极化特征的。本文采用的是基于Wishart距离的提取方法。

对于哨兵-1号的双极化数据，其提取流程较全极化简单，这也是本文使用的数据。下面将分步介绍提取流程：

预处理

首先要对两个时相的图像VV和VH波段分别进行配准，目的是使两个图像的每一个像元（像素）都彼此对齐。该部分参考较多，不再赘述 : )

然后是影像的读取、多视和滤波——SAR影像的距离向和方位向分辨率并不相同，因此需要进行多视；而滤波可以降低图像的噪声——本数据已经经过了滤波，故不再演示。下面是Matlab代码：

clc;clear;
%读取第一个时相SAR数据文件
file12_1 = "/20210428/s12.bin";
file22_1 = "/20210428/s22.bin";
s12_1 = multibandread(file12_1,[3822,15021*2,1],'float32',0,'bsq','ieee-le');
s22_1 = multibandread(file22_1,[3822,15021*2,1],'float32',0,'bsq','ieee-le');
[a1,b1] = ComplexSeparate(s12_1);
[a2,b2] = ComplexSeparate(s22_1);
S12_1 = complex(a1,b1);
S22_1 = complex(a2,b2);

T1 = S22_1.*conj(S22_1);
T2 = S12_1.*conj(S12_1);
T3 = S22_1.*conj(S12_1);

Tm1_1 = multilookProcessing(T1,1,4);
Tm2_1 = multilookProcessing(T2,1,4);
Tm3_1 = multilookProcessing(T3,1,4);

%读取第二个时相SAR数据文件
file12 = "/20210522/s12.bin";
file22 = "/20210522/s22.bin";
s12 = multibandread(file12,[3822,15021*2,1],'float32',0,'bsq','ieee-le');
s22 = multibandread(file22,[3822,15021*2,1],'float32',0,'bsq','ieee-le');
[a1,b1] = ComplexSeparate(s12);
[a2,b2] = ComplexSeparate(s22);
S12 = complex(a1,b1);
S22 = complex(a2,b2);

T1 = S22.*conj(S22);
T2 = S12.*conj(S12);
T3 = S22.*conj(S12);

Tm1_2 = multilookProcessing(T1,1,4);
Tm2_2 = multilookProcessing(T2,1,4);
Tm3_2 = multilookProcessing(T3,1,4);

function [a,b]=ComplexSeparate(m)
A=size(m);
a=m(:,1:2:A(1,2));
b=m(:,2:2:A(1,2));
end

经过上面的操作，可以得到两个时相多视后的图像矩阵Tm1、Tm2、Tm3，分别代表每一个像元的极化相干矩阵的一个元素，下面将会进行组合。

差异图提取

原理

根据论文（张永红，等:一种基于极化距离测度的 SAR 变化检测方法），变化检测与图像分类其实是相似的。极化分类是将极化特性相似的像元归为同一类别，而变化检测则是要从前后时相影像上找出极化特性相差较大的像元。因此，极化分类中使用的最大似然分类原理也可以被用在变化检测上。

又因为多视处理的极化 SAR影像的相干矩阵T服从复Wishart分布，此时，根据最大似然分类原理，可以推导出如下的样本相干矩阵与类中心相干矩阵之间的 Wishart距离\(d (T，V_m )=ln｜Vm｜ +tr (V_m^{-1}T)\)，任一像元可根据其与所有类别中心相干矩阵的 Wishart距离最小值判断其归 属 于哪 一类。

而在变化检测中，该公式还存在一定的不足，其中存在与变化无关的信息。因此，需要提出一种经过变形的公式来表征适用于变化检测的极化距离测度 ，即：\(D(i,j)={1\over 2}tr(T_1^{-1} T_2 + T_2^{-1} T_1)-3\)。式中，T1、T2分别表示前后时相像元的极化相干矩阵；减3是为了保证同质像元之间的距离为零，以满足距离测度的非负性。

实现

按照多视后图像的分辨率设置矩阵D用于存放差异图的像素

通过循环遍历每一个像元，组成像元的极化相干矩阵，带入公式\(D(i,j)={1\over 2}tr(T_1^{-1} T_2 + T_2^{-1} T_1)-3\)计算D的每一个元素。

%生成极化距离图像D
D = zeros(3822,3755);
for i = 1:3822
    for j = 1:3755
        T_1 = [Tm1_1(i,j),Tm3_1(i,j);Tm3_1(i,j)',Tm2_1(i,j)]; %第一个时相ij像元的极化相干矩阵
        T_2 = [Tm1_2(i,j),Tm3_2(i,j);Tm3_2(i,j)',Tm2_2(i,j)]; %第二个时相ij像元的极化相干矩阵
        D(i,j) = 1/2 * trace(T_1^-1 * T_2 + T_2^-1 * T_1)-3;
    end
end

imtool(D)

D就是所求差异图

效果展示