基于多重分形谱的SAR图像二次成像方法与流程

本发明涉及sar图像处理技术领域，具体涉及一种基于多重分形谱的sar图像二次成像方法。

背景技术：

合成孔径雷达(syntheticapertureradar，sar)具有全天候、全天时的成像特点、潜在的高分辨率和巨大的应用前景。从理论上来讲，sar成像可以理解为：载有天线的系统平台在空间中沿着某个运动轨迹运动，在运动过程中，发射某个载频的微波信号，然后地面反射回来的回波由系统平台的接收机接收，当平台走过一个合成孔径的距离后，将积累的回波信号在处理器中进行成像处理形成雷达图像。总体来说sar成像就是天线和观测区域之间的相对运动，利用在运动方向上合成的虚拟孔径产生合成孔径图像，从而获得雷达观察区域的图像。sar成像基本原理一定程度限制了对雷达电磁回波所蕴含信息的提取能力。sar图像中待检测目标的种类繁多，且目标所处的环境较为复杂，噪声干扰强。复杂场景、极低信噪比条件下sar图像解译是雷达领域面临的世界性难题，受到学术界的广泛关注。

在过去的几十年里，复杂的自然系统中出现的尺度不变性和分形行为已经被实验所证实。分形维数(fractaldimension，fd)和多重分形谱(multifractalspectrum，mfs)是分形几何分支中的基础性概念，也是随机分形信号分析和处理的核心问题和概念之一。多重分形是定义在分形结构上的有多个标度指数所组成的集合，对信号进行多重分形特性分析是为了从系统的局部出发来研究其整体的特征，并借助统计物理学的方法来讨论特征参量的概率测度的分布规律。本申请提出基于多重分形谱的sar图像二次成像方法，有望为复杂场景、极低信噪比条件下sar图像解译提供新的途径和方法。

技术实现要素：

针对复杂场景、极低信噪比条件下，sar图像解译存在的问题，本申请提供一种基于多重分形谱的sar图像二次成像方法，包括步骤：

对sar图像像元进行二维局部延迟自相关分析，获得四维数组，所述四维数组包括原始自变量和二维延迟分量或伪魏格纳时频展开；

保留所述原始自变量，对所述二维延迟分量或二维频率分量进行二维多重分形谱分析，获得多重分形谱；

遍历原始sar图像各空间点，对所述多重分形谱进行处理，得到所述sar图像的多重分形谱sar二次成像结果。

进一步优选的，所述对sar图像像元进行二维局部延迟自相关分析之前包括：建立一个m×n×n_1×n_2的四维数组，其中，m×n为sar图像的大小，n_1、n_2分别为方位向和距离向的延迟尺度，每一个像元一一对应一幅n_1×n_2大小的子图像。

进一步优选的，采用瞬时自相关算法对所述sar图像像元进行二维局部延迟自相关分析，具体包括：

r(x，r，x′，r′)＝i(x+x′，r+r′)×i^*(x-x′，r-r′)

其中，i(x，r)为sar图像的某一像元，x＝1，2，...，m，r＝1，2，...n，x′和r′分别代表方位向和距离向的延迟，x′＝1，2，...，n1，r′＝1，2，...n2，*代表共轭，r(x，r，x′，r′)为经瞬时自相关算法分析获取的四维数组r(x，r，x′，r′)的大小为m×n×n1×n2，n1、n2分别为方位向和距离向的延迟尺度。

进一步优选的，采用2d-pwvd算法对所述sar图像像元进行二维局部延迟自相关分析，具体包括：

其中，i(x，r)为sar图像的某一像元，x＝1，2，...，m，r＝1，2，...n，*代表共轭，r(x，r，fx′，fr′)为i(x，r)的二维维格纳分布，k和l分别为方位向和距离向的延迟分量，(fx′，fr′)为频率分量。

进一步优选的，所述方位向和距离向的延迟尺度n1、n2均为奇数，且n1＝n2。

进一步优选的，采用基于配分函数法的二维多重分形分析法、二维多重分形降趋波动分析法或二维多重分形降趋移动平均法对所述二维延迟分量或二维频率分量进行二维多重分形谱分析。

进一步优选的，采用二维多重分形降趋波动分析法对所述二维延迟分量或二维频率分量进行二维多重分形谱分析，具体包括步骤：

采用二维阵列x(i，j)表示待处理的二维分量所构成的矩阵(对应于原始sar图像中某个局部点)，其中，i＝1，2，...，n1′，j＝1，2，...，n2′，x(i，j)被划分成n1s×n2s个互不相交的有着相同大小的s×s的方块区间，n1s＝[n′1/s]，n2s＝[n′2/s]，每个方块用xv，w表示，xv，w(i，j)＝x(l1+i，l2+j)，1≤i，j≤s，其中，l1＝(v-1)s1，v＝1，2，...，n1s，l2＝(w-1)s2，w＝1，2，...，n2s；

对每一个方块xv，w进行累加计算，获取uv，w(i，j)：

其中，i，j＝1，2，...，s，uv，w是一个曲面；

采用二变量多项式降趋拟合函数(i，j)对曲面uv，w的趋势进行拟合；

计算残差矩阵：

根据残差矩阵εv，w(i，j)的样本方差定义方块xv，w的降趋函数fv，w(s)：

对所有的方块计算全局的降趋波动函数：

改变尺度s的数值，范围从smin≈6到之间变化，确定尺度s和波动函数fq(s)之间的尺度关系：

fq(s)～s^h(q)

对每一个q值，获得相应的传统的质量指数函数τ(q)表达式：

τ(q)＝qh(q)-df

其中df为拓扑维数，对于二维信号序列，df＝2；

通过勒让德变换计算出奇异性指数函数α(q)和多重分形谱f(α)：

α(q)＝dτ(q)/dq；

f(α)＝infq[qα(q)-τ(q)]。

进一步优选的，所述二变量多项式降趋拟合函数的表达式为：

或者，

或者，

或者，

或者，

其中1≤i，j≤s，a，b，c，d，e，f是待定的自由参数。

本发明提供的基于多重分形谱的sar图像二次成像方法，具有的有益效果是：使得sar图像可以从两个不同角度进行描述：图像角度fα(x，r)与多重分形谱角度fx，r(α)，从而，为sar图像分析和特征提取提供更丰富的细节，为sar图像解译提供一种新思路。

附图说明

图1为基于多重分形谱的sar图像二次成像方法流程图；

图2为本实施例中使用的sar图像，(a)纯海面图像，(b)舰船目标图像；

图3为纯海面的多重分形谱sar二次成像结果，(a)(b)(c)(d)分别对应不同奇异性指数强度；

图4为舰船目标图像的多重分形谱sar二次成像结果，(a)(b)(c)(d)分别对应不同奇异性指数强度；

图5为sar图像从多重分形谱角度得到的处理结果，(a)纯海面图像多重分形谱分布情况，(b)舰船目标图像多重分形谱分布情况。

具体实施方式

下面通过具体实施方式，结合附图对本发明作进一步详细说明。

本发明提供一种基于多重分形谱的sar图像二次成像方法，其流程图如图1所示，具体包括以下步骤。

s100：对sar图像像元进行二维局部延迟自相关分析，获得四维数组，四维数组包括原始自变量和二维延迟分量或伪魏格纳时频展开。

s200：保留所述原始自变量，对二维延迟分量或二维频率分量进行二维多重分形谱分析，获得多重分形谱。

s300：遍历原始sar图像各空间点，对多重分形谱进行处理，得到sar图像的多重分形谱sar二次成像结果。

上述方法针对复杂场景、极低信噪比条件下，为sar图像解译提供一种新的思路，使得sar图像可以从两个不同角度进行描述：图像角度fα(x，r)与多重分形谱角度fx，r(α)，从而，为sar图像分析和特征提取提供更丰富的细节。

下面对上述各步骤进行详细说明。

在步骤s100中，对sar图像像元进行二维局部延迟相关分析之前，需要先建立一个m×n×n1×n2的四维数组，其中，m×n为sar图像的大小，n1、n2分别为方位向和距离向的延迟尺度，每一个像元一一对应一幅n1×n2大小的子图像。

本实施例提供两种方式实现二维局部延迟自相关分析，该两种实现方式具体如下：

一种实现方式是，通过时域实现二维局部延迟自相关分析，具体的，采用瞬时自相关算法对sar图像像元进行二维局部延迟自相关分析，该种实现方式的具体过程如下：

具体的，瞬时自相关算法可以定义为：

r(x，r，x′，r′)＝i(x+x′，r+r′)·i^*(x-x′，r-r′)；

其中，i(x，r)为sar图像的某一像元，x＝1，2，...，m，r＝1，2，...n，x′和r′分别代表方位向和距离向的延迟，x′＝1，2，...，n1，r′＝1，2，...n2，*代表共轭，r(x，r，x′，r′)为瞬时自相关算法的结果，也即是，r(x，r，x′，r′)为经瞬时自相关算法分析获取的四维数组。

优选的，方位向和距离向的延迟尺度n1、n2均为奇数，且n1＝n2，同时，n1、n2越大，处理计算量越大，且多元成分的干扰越严重；但n1、n2过小，多重分形谱估计受图像信噪比的影响越明显，更易出现误差。本实施例中，经过实验综合论证，n1＝n2＝17。

另一种实现方式是，通过频域实现二维局部延迟自相关分析，具体的，采用2d-pwvd算法(即图像瞬时自相关的傅立叶变换)对所述sar图像像元进行二维局部延迟自相关分析，具体包括：

其中，i(x，r)为sar图像的某一像元，x＝1，2，...，m，r＝1，2，...n，*代表共轭，r(x，r，fx′，fr′)为i(x，r)的二维维格纳分布，k和l分别为方位向和距离向的延迟分量，(fx′，fr′)为频率分量。

在步骤s200中，可以采用基于配分函数法的二维多重分形分析法、二维多重分形降趋波动分析法或二维多重分形降趋移动平均法对所述二维延迟分量或二维频率分量进行二维多重分形谱分析。

优选的，本实施例采用二维多重分形降趋波动分析法对所述二维延迟分量或二维频率分量进行二维多重分形谱分析，具体包括步骤：

步骤1：采用二维阵列x(i，j)表示待处理的二维分量所构成的矩阵(对应于原始sar图像中某个局部点)。

其中，i＝1，2，...，n1′，j＝1，2，...，n2′，x(i，j)被划分成n1s×n2s个互不相交的有着相同大小的s×s的方块区间，n1s＝[n′1/s]，n2s＝[n′2/s]，每个方块用xv，w表示，xv，w(i，j)＝x(l1+i，l2+j)，1≤i，j≤s，其中，l1＝(v-1)s1，v＝1，2，...，n1s，l2＝(w-1)s2，w＝1，2，...，n2s。

步骤2：对每一个方块xv，w进行累加计算，获取uv，w(i，j)。

该步骤中的具体计算公式如下：

其中，i，j＝1，2，...，s，uv，w是一个曲面；

步骤3：采用二变量多项式降趋拟合函数对曲面uv，w的趋势进行拟合。

该步骤中，二变量多项式降趋拟合函数的表达式可以采用如下任一种表达式：

或者，

或者，

或者，

或者，

其中1≤i，j≤s，a，b，c，d，e，f是待定的自由参数，这些自由参数可以利用最小二乘法得到具体数值。

步骤4：计算残差矩阵。

该步骤中的残差矩阵可以采用如下公式计算：

步骤5：根据残差矩阵εv，w(i，j)的样本方差定义方块xv，w的降趋函数fv，w(s)。

步骤6：对所有的方块计算全局的降趋波动函数。

步骤6：改变尺度s的数值，范围从smin≈6到之间变化，确定尺度s和波动函数fq(s)之间的尺度关系。

fq(s)～s^h(q)

对每一个q值，获得相应的传统的质量指数函数τ(q)，其表达式如下：

τ(q)＝qh(q)-df

其中df为拓扑维数，对于二维信号序列，df＝2。

步骤7：通过勒让德变换计算出奇异性指数函数α(q)和多重分形谱f(α)：

α(q)＝dτ(q)/dq；

f(α)＝infq[qα(q)-τ(q)]。

在步骤s300中，对步骤s200中获得的奇异性指数函数α(q)和多重分形谱f(α)进行处理，得到多重分形谱sar二次成像f(x，r，α)。

通过对f(x，r，α)进行处理，获取奇异性强度函数的最大值αmax与最小值αmin，将其构成的区间[αmin，αmax]均分成9个子区间，按照子区间将奇异性强度函数划分，并计算对应多重分形谱的平均值，可得到不同奇异性指数对应sar图像fα(x，r)。

例如，纯海面图像和舰船目标图像的sar图如图2所示，纯海面sar图像部分子区间的成像效果如图3所示，其对应的奇异性指数范围分别为[0.2559，0.8826]、[0.8826，1.5092]、[1.5092，2.1359]、[2.7625，3.3891]。舰船目标sar图像部分子区间的成像效果如图4所示，其对应的奇异性指数范围分别为[0.5671，1.3501]、[1.3501，2.1330]、[2.1330，2.9159]、[2.9159，3.6989]。奇异性指数与分形维数有非常紧密的联系。而对于自然表面，其分形维数与表面粗糙度有直观的直接关系：一个近乎光滑的分形表面其分形维数略大于2，因为它在经典模型表面上趋于平坦；相反地，一个非常粗糙的表面有一个接近3的分形维数。从图3、图4可以分析得到，海面区域的奇异性指数跨度较小，而舰船区域的奇异性指数跨度较大，具有较明显的区分度。

通过对f(x，r，α)进行处理，保留x，r变量，获取δα(δα＝max(fx，r(α))-min(fx，r(α)))作为z轴变量，x轴与y轴分别对应x，r变量，得到多重分形谱的分布情况，可得到不同散射点的多重分形特征谱fx，r(α)，如图5所示。

因此，sar图像处理结果f(x，r，α)可以从两个不同角度进行描述：不同奇异性指数对应sar图像fα(x，r)和不同散射点的多重分形特征谱fx，r(α)。

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。