一种角度优化处理的全极化SAR目标分解方法与流程

本发明属于合成孔径雷达(syntheticapertureradar，sar)图像处理技术领域，更进一步涉及一种角度优化处理的全极化sar目标分解方法。

背景技术：

极化sar通过发射和接收电磁波，能够全天时、全天候进行对地观测，为地物目标的解译提供了丰富的散射信息。通过有效的方法确定和理解散射机制，能够极大地提高极化sar图像在目标检测、地物分类等任务中的能力，使极化sar在应急抢险救灾、城市分析规划等领域发挥更重要的作用。

考虑数据互易性，极化sar数据可表示为后向散射矩阵s2、极化协方差矩阵c3、极化相干矩阵t3，三者具有唯一对应关系，可相互转化。以极化相干矩阵t3矩阵为例，将t3矩阵分解并与标准散射体(如平面、二面角、螺旋体等)的极化相干矩阵进行比较，可以得到目标的散射机理，即利用目标分解方法，可以基于切合实际物理约束地进行目标散射机制的解译，因此在极化sar数据处理与目标解译领域受到广泛关注。

极化sar目标分解方法主要可以分为基于特征值的分解方法和基于模型的分解方法。基于特征值的分解方法具有明确的数学意义，且分解得到的特征值唯一；基于模型的分解方法具有较好的物理可解释性，然而其求解过程基于不同的假设条件和模型，因此对基于模型的分解方法的研究始终是极化sar领域的热点之一。freeman等人于1998年最先提出freeman-durden分解，将目标t3矩阵分解为表面散射、偶次散射和体散射的t3矩阵之和，又称三分量分解。yamaguchi等人于2005年提出yamaguchi分解，在三分量的基础上加入螺旋散射分量，形成四分量分解方法(y4o)。上述两种分解方法是基于模型的分解方法中最为经典的。随后，在yamaguchi四分量的基础上，学者们针对体散射模型的普适性、体散射分量估计过高等问题对模型进行研究和修改，先后提出了去取向的yamaguchi四分量分解(y4r)，扩展体分量的yamaguchi四分量分解(s4r)，两次角度补偿的yamaguchi四分量分解(g4u)，泛化的体分量模型(gsvm)等分解方法。

然而，现有的基于模型的分解方法仍存在三大问题。其一，基于模型的分解方法往往采用求解方程的方式获得各个成分分量的值，由于未知数小于方程总数，往往仅利用了t3矩阵中的9个变量中的5至8个，t3矩阵中的散射信息没有得到充分利用；其二，体散射分量估计过高；其三，分解过程中存在许多假设和条件分支，使处理流程复杂化。因此，如何充分利用极化sar原始数据中的全部信息，有效估计各成分分量，同时简化处理过程中的判断和条件分支，是基于模型的分解方法需要解决的问题。

技术实现要素：

本发明针对现有技术对极化sar数据进行模型分解时信息利用不充分、体散射分量估计过高，以及处理流程条件分支繁多的问题，提出一种角度优化处理的全极化sar目标分解方法。通过构建分解损失，以优化求解的方法代替传统求解方程组的方法获得各成分分量所占比重，充分利用t3矩阵中的全部信息。

本发明的技术方案为一种角度优化处理的全极化sar目标分解方法，包括以下步骤，

步骤1，读入极化sar数据，若源数据为极化相干矩阵t3，进入步骤2；若源数据为后向散射矩阵s2或极化协方差矩阵c3，则将数据转换为极化相干矩阵t3；

步骤2，数据滤波，得到滤波后的极化sar数据；

步骤3，进行去取向的yamaguchi四分量分解，获得初始的表面散射分量fs、偶次散射分量fd、体散射分量fvol、螺旋散射分量fc、偶次散射参数α、表面散射参数β，以及初始的旋转角度θ，将泛化体散射分量fv初始化为当前的体散射分量fvol，并初始化二面角散射贡献的扩展体散射分量fod＝0；

步骤4，利用初始的旋转角度θ，对t3矩阵通过旋转变换获得θ角下的t3矩阵t3(θ)，定义分解模型如下，

将t3(θ)矩阵分解为表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、扩展体散射分量fod、螺旋散射分量fc的加权和，如下式，

t3(θ)＝fstodd+fdtdbl+fvtvol+fodtod+fcthlx+trem

其中，todd为表面散射模型矩阵、tdbl为偶次散射模型矩阵、tvol为泛化体散射模型矩阵、tod为扩展体散射模型、thlx为螺旋散射模型、trem为分解残余量矩阵。

构建并计算分解模型的重构损失如下，

其中，t'ij为残余量矩阵第i行j列的元素：

t'11＝fs+fd|α|²+afv-t11(θ),

t'22＝fs|β|²+fd+bfv+1/2fc+7/15fod-t22(θ),

t'33＝cfv+8/15fod+1/2fc-t33(θ),

re[t'12]＝fsβ^*+fdre[α]+dfv-re[t12(θ)],

im[t'12]＝fdim[α]-im[t12(θ)],

re[t'13]＝re[t13(θ)],im[t'13]＝im[t13(θ)],

re[t'23]＝re[t23(θ)],im[t'23]＝1/2fc-im[t23(θ)]

其中，re[·],im[·]分别表示取实部和取虚部，*表示共轭运算，

重构损失l中，反应了四分量分解后的信息与原始t3矩阵信息损失，反应了四分量分解后各分量功率之和与总功率之差；

角度优化处理的优化问题为表示求解最优的一组参数{fs,fd,fv,fod,fc,α,β,θ}，使得重构损失l最小；

步骤5，利用梯度下降法求解优化问题，得到最优的一组参数{fs,fd,fv,fod,fc,α,β,θ}；

步骤6，利用步骤5获得的表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β、旋转角θ以及二面角散射贡献的扩展体散射分量fod，提取表面散射功率ps、偶次散射功率pd、体散射功率pv、螺旋散射功率pc，获得四分量分解结果，得到sar图像的目标特性。

而且，步骤2中，采用滤波窗口大小为3x3的boxcar滤波器，对原始极化sar数据进行滤波，得到滤波后的极化sar数据。

而且，步骤5的实现方式如下，

计算重构损失l对表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β、旋转角θ以及二面角散射贡献的扩展体散射分量fod的导数；

利用导数值，根据梯度下降算法更新更新表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β、旋转角θ以及二面角散射贡献的扩展体散射分量fod；

判断更新的所有参数是否满足边界条件，若不满足则将其设置为边界值；

利用更新后的参数重新计算分解模型的重构损失l，若达到收敛条件或达到迭代次数，则进入步骤6，否则重复本步骤直至达到收敛条件或达到迭代次数。

而且，步骤6中提取表面散射功率ps、偶次散射功率pd、体散射功率pv、螺旋散射功率pc，实现如下，

ps＝fs(1+|β|²),pd＝fd(1+|α|²),pv＝fv+fod,pc＝fc

本发明的优点是：

1、采用基于优化的方式进行四分量分解，信息利用更加充分。相较于传统基于求解方程组的分解方法，本发明充分利用了原始极化sar数据t3矩阵中的全部9个变量值,即大小为3x3的极化相干矩阵t3的全部9个元素，而传统的yamaguchi分解仅利用了5个变量值，去取向的yamaguchi分解利用了6个变量值，其他分解方法大多也仅利用了t3矩阵中的5至8个变量值，因此本发明的方法信息利用更加充分；

2、采用泛化体散射模型，更加精确地刻画目标体散射情况。传统四分量分解方法中，往往基于某种先验假设采用某一种体散射模型，在实际应用中，有多个体散射模型可以进行选择，难以选择最优模型。针对此问题，本发明采用泛化体散射模型tvol，该模型通引入参数τ，使得体散射模型tvol的各元素根据散射功率的值计算获得，当参数τ为不同值时，对应的体散射模型不同，使整个分解模型更加符合真实的地物散射机理，也能够更好的对复杂场景下的体散射模型进行建模。

3、采用扩展体散射模型，即将体散射分量分为由植被贡献的体散射分量，即泛化体散射分量fv和由二面角散射贡献的体散射分量，即扩展体散射分量fod。传统的四分量分解模型中，由于未对具有较大入射角的二面体散射进行建模，会将具有较大入射角的二面体散射分量计入体散射分量，导致体散射估计过高问题。引入扩展体散射模型，对具有较大入射角的二面体散射分量建模，可有效降低体散射分量的估计量。扩展体散射模型已应用于s4r、g4u等分解方法。然而，在s4r和g4u分解方法中，引入了判决条件判断体散射分量由植被还是二面角散射贡献，在本发明中，利用优化方法进行求解，可自适应的学习到植被与二面角散射贡献的体散射分量的最优比例，不需要引入人工判决条件。

附图说明

图1是本发明实施例的流程图；

图2是美国旧金山地区的光学影像及四种现有的分解方法和本专利方法对airsar数据的分解结果图；

图3是日本新潟地区的光学影像及四种现有的分解方法和本专利方法对pi-sar数据的分解结果图；

图4是美国旧金山地区的y4o伪彩图及四种现有的分解方法和本专利方法对airsar数据三处局部区域的分解结果雷达图；

图5是日本新潟地区的y4o伪彩图及四种现有的分解方法和本专利方法对pi-sar数据两处局部区域的分解结果雷达图；

图6是美国旧金山地区的光学影像及三种现有的分解方法和本专利方法对airsar数据的分解重构损失图；

图7是日本新潟地区的光学影像及三种现有的分解方法和本专利方法对pi-sar数据的分解重构损失图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

本发明实施例提供的一种角度优化处理的全极化sar目标分解方法流程如附图1所示，包括以下步骤：

步骤1，读入极化sar数据，以极化相干矩阵t3＝(tij)3×3为例，tij为t3矩阵的各个元素。

具体实施时，如果原始的极化sar数据是散射矩阵s2、极化协方差矩阵c3，可以先转换成t3矩阵后，按本实施例流程进行处理。

步骤2，数据滤波，得到滤波后的极化sar数据。

为消除相干斑对极化sar图像的影响，实施例优选采用滤波窗口大小为3x3的boxcar滤波器，对原始极化sar数据进行滤波，得到滤波后的极化sar数据。boxcar滤波器为现有技术，其原理为将窗口内变量值的平均值赋给窗口中心点：

其中x(i,j)为原始t3矩阵第i行j列的变量值，x(i+r,j+c)为原始t3矩阵第i+r行j+c列的变量值，ws为空间窗尺寸，此处为3，x'(i,j)为滤波后t3矩阵第i行j列的变量值。

步骤3，利用去取向的yamaguchi四分量分解初始化变量值。

对步骤2中获得的滤波后的极化sar数据t3矩阵进行去取向的yamaguchi四分量分解。矩阵进行去取向的yamaguchi四分量分解为现有技术，可参见yamaguchiy,satoa,boernerwm,etal.four-componentscatteringpowerdecompositionwithrotationofcoherencymatrix[j].ieeetransactionsongeoscienceandremotesensing,2011,49(6):2251-2258.本发明不予赘述。经过去取向的yamaguchi四分量分解后，获得初始的表面散射分量fs、偶次散射分量fd、体散射分量fvol、螺旋散射分量fc、偶次散射参数α、表面散射参数β，以及初始的旋转角度θ，将泛化体散射分量fv初始化为当前的体散射分量fvol，并初始化二面角散射贡献的扩展体散射分量fod＝0；

步骤4，利用初始的旋转角度θ，对t3矩阵通过旋转变换获得θ角下的t3矩阵t3(θ)，定义分解模型如下：将t3(θ)矩阵分解为表面散射fs、偶次散射fd、泛化体散射fv、扩展体散射fod、螺旋散射fc五个分量的加权和形式，即：

t3(θ)＝fstodd+fdtdbl+fvtvol+fodtod+fcthlx+trem

其中todd为表面散射模型矩阵、tdbl为偶次散射模型矩阵、tvol为泛化体散射模型矩阵、tod

为扩展体散射模型、thlx为螺旋散射模型、trem为分解残余量矩阵，其具体形式分别为：

其中，表面散射参数β为一实数，|β|≤1，偶次散射参数α为一复数，|α|≤1,j表示单位虚数，*表示共轭运算，|·|表示取模运算，参数t3矩阵通过旋转变换获得θ角下的t3矩阵t3(θ)的具体方法为：

其中u3(θ)为旋转矩阵，为旋转矩阵的共轭转置矩阵，tij为t3矩阵的各个元素。由上式可得：

t11(θ)＝t11

t22(θ)＝t22cos²2θ+t33sin²2θ+re[t23]sin4θ

t33(θ)＝t33cos²2θ+t22sin²2θ-re[t23]sin4θ

其中re[·],im[·]分别表示取实部和取虚部，*表示共轭运算。

构建并计算分解模型的重构损失：

其中，t'ij为残余量矩阵第i行j列的元素：

t'11＝fs+fd|α|²+afv-t11(θ),

t'22＝fs|β|²+fd+bfv+1/2fc+7/15fod-t22(θ),

t'33＝cfv+8/15fod+1/2fc-t33(θ),

re[t'12]＝fsβ^*+fdre[α]+dfv-re[t12(θ)],

im[t'12]＝fdim[α]-im[t12(θ)],

re[t'13]＝re[t13(θ)],im[t1'3]＝im[t13(θ)],

re[t'23]＝re[t23(θ)],im[t'23]＝1/2fc-im[t23(θ)]

其中，

重构损失l中，反应了四分量分解后的信息与原始t3矩阵信息损失，反应了四分量分解后各分量功率之和与总功率之差。

角度优化处理的优化问题可定义为：

即：求解最优的一组参数{fs,fd,fv,fod,fc,α,β,θ}，使得重构损失l最小。

步骤5，利用梯度下降算法求解优化问题，具体步骤如下：

第一步，计算重构损失l对表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β、旋转角θ以及二面角散射贡献的扩展体散射分量fod的导数。

为了以最小的信息损失获得最优的四分量分解结果，本专利基于步骤4中构建的重构损失l，利用梯度下降算法求解四分量，首先计算重构损失l对各参数的导数：

第二步，利用第一步中计算获得的导数值，根据梯度下降算法更新表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β、旋转角θ以及二面角散射贡献的扩展体散射分量fod：

其中，a←b表示将b的值赋给a，λ为学习率，本实施例中优选取λ＝0.01。

第三步，判断第二步中更新的所有参数是否满足边界条件，若不满足则将其设置为边界值。

设span为目标总功率，c为判决指数，c>0表示体散射功率由泛化体散射分量贡献，c<＝0表示体散射功率由有朝向的二面角散射分量贡献，按照下述流程依次检查参数并修正：

span＝t11+t22+t33，c＝2re[fsβ+fdre[α]]；

1.若fv＜0，则fc＝0，且：

①若fs＜0且fd＜0，则：fs＝fd＝fv＝0,fod＝span,β＝1,α＝-1；

②若fs＜0，fd≥0且c＞0，则fs＝fod＝0,fv＝span-fd(1+|β|²),α＝-1；

③若fs＜0，fd≥0且c≤0，则fs＝fv＝0,fod＝span-fd(1+|β|²),α＝-1；

④若fs≥0，fd＜0且c＞0，则fd＝fod＝0,fv＝span-fs(1+|α|²),β＝1；

⑤若fs≥0，fd＜0且c≤0，则fd＝fv＝0,fod＝span-fs(1+|α|²),β＝1；

⑥若fs≥0，fd≥0且c＞0，则fod＝0，fv＝span-fs(1+|α|²)-fd(1+|β|²)；

⑦若fs≥0，fd≥0且c≤0，则fv＝0，fod＝span-fs(1+|α|²)-fd(1+|β|²)；

2.若fv≥0，则：

①若fc＜0，则fc＝0；

②若fv+fod+fc＞span，则：

(1)若c＞0，则fod＝0,fv＝span-fc,fs＝0,fd＝0,β＝1,α＝-1；

(2)若c≤0，则fv＝0,fod＝span-fc,fs＝0,fd＝0,β＝1,α＝-1；

③若fv+fod+fc≤span，则：

(1)若fs＜0且fd＜0，则fs＝fd＝0,fv＝span-fc-fod,β＝1,α＝-1；

(2)若fs＜0且fd≥0，则fs＝0,fd＝0.5*(span-fc-fv-fod),α＝-1；

(3)若fs≥0且fd＜0，则fd＝0,fs＝0.5*(span-fc-fv-fod),β＝1；

3.若β＞1，则令β＝1；若β＜-1，则令β＝-1；

4.若|α|²＞1且re[α]≥0，则令α＝1；若|α|²＞1且re[α]＜0，则令α＝-1.

第四步，利用第三步中更新后的参数重新计算分解模型的重构损失l，计算方法同步骤4。若达到收敛条件或达到迭代次数，则进入步骤6，否则返回第一步，重复步骤5直至达到收敛条件或达到迭代次数。

本实施例中，优选收敛条件为：当前的重构损失l小于0.001倍初始的重构损失l；迭代次数为：10000次。

步骤6，计算四分量分解结果，获取sar图像的目标特性。

利用步骤5获得的表面散射分量fs、偶次散射分量fd、泛化体散射分量fv、螺旋散射分量fc、偶次散射参数α、表面散射参数β以及二面角散射贡献的扩展体散射分量fod，提取表面散射功率ps、偶次散射功率pd、体散射功率pv、螺旋散射功率pc，获得四分量分解结果：

ps＝fs(1+|β|²),pd＝fd(1+|α|²),pv＝fv+fod,pc＝fc

其中，表面散射功率ps代表了目标布拉格表面散射的强度；偶次散射功率pd代表了目标二面角散射的强度；体散射功率pv代表了体散射的强度；螺旋散射功率pc代表了目标非均匀散射的强度。以角度优化处理获得的四分量分解结果能够减少sar图像信息损失，克服体分量估计过高等技术问题，更加准确地从物理意义出发刻画了sar图像的目标特性，有利于sar图像目标解译技术的发展。

具体实施时，以上流程可采用计算机软件技术实现自动运行，运行流程的系统装置也应当在本发明的保护范围内。

为便于实施参考起见，下面结合airsarl段美国旧金山地区的全极化sar数据和pi-sarl波段日本新潟地区的全极化sar数据对本发明的有效性做进一步说明。

1.实验数据

airsarl段美国旧金山地区的全极化sar数据，图像尺寸为900x1024，地物类型主要为水体、建筑和林地。pi-sarl波段日本新潟地区的全极化sar数据，图像尺寸为1250x2000，地物类型主要为水体、建筑和农田。

2.分解结果

附图2为美国旧金山地区光学图像，y4o、y4r、s4r、g4u、本专利方法对airsarl段美国旧金山地区的全极化sar数据四分量分解中偶次散射分量和体散射分量的图像。可以由附图2中看出，本专利方法中的建筑区域的偶次散射分量有明显提升，特别是在附图2中标识的a区域，由于建筑朝向与雷达照射角度夹角较大，在传统分解方法中体散射分量过高，而在本专利方法的结果中，体散射分量降低。附图3为日本新潟地区的光学图像，y4o、y4r、s4r、g4u、本专利方法对pi-sarl波段日本新潟地区sar数据四分量分解中偶次散射分量和体散射分量的图像。附图4、5分别为y4o、y4r、s4r、g4u、本专利方法对airsarl段美国旧金山地区的全极化sar数据、pi-sarl波段日本新潟地区的全极化sar数据中局部区域四分量分解雷达图。可以看出，本专利方法有效的抑制了体散射分量估计过高的问题。

3.信息损失分析

附图6、7分别为y4o、y4r、s4r、本专利方法对airsarl段美国旧金山地区的全极化sar数据、pi-sarl波段日本新潟地区的全极化sar数据四分量分解后的重构损失l。可以看出，本专利方法的信息损失l更小，原始sar数据的信息利用更加充分。

本发明的创新之处在于，以优化思路求解极化sar四分量分解问题，优化求解去取向角度，充分利用了原始极化sar数据信息，更准确的求解不同散射机制所占比重；采用泛化体散射模型，可以自适应的刻画不同场景目标的体散射情况；采用扩展体分量模型，解构体散射模型，将其区分为由植被贡献的泛化体散射分量和由有朝向的二面角贡献的扩展体散射分量，可以有效地克服体分量估计过高问题；同时，通过优化策略可自适应的学习到植被与二面角散射贡献的体散射分量的最优比例。

以上所述仅为本发明的优选实施例而已，并不用以限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。