CN104297751A - 一种分布式卫星sar系统的主辅sar图像融合方法 - Google Patents

Abstract

本发明属于SAR图像方位向分辨率提升技术领域，特别涉及一种分布式卫星SAR系统的主辅SAR图像融合方法。其具体步骤为：多普勒谱包络校正；雷达几何法图像配准；地形相位反演；(SAR图像最大对比度法估计配准误差及残留相位偏差；主辅SAR图像精配准及相位偏差校正；多普勒谱拼接。本发明具有在沿航向长基线分布式卫星SAR系统的多普勒谱重叠度近似为零情况下，能够对主辅SAR图像对进行有效的图像融合以提高方位向分辨率，满足高质量InSAR处理的性能要求。

Description

一种分布式卫星SAR系统的主辅SAR图像融合方法

技术领域

本发明属于SAR图像方位向分辨率提升技术领域，特别涉及一种分布式卫星SAR系统的主辅SAR图像融合方法，当沿航向长基线分布式卫星SAR系统的多普勒谱重叠度近似为零，方位预滤波后SAR图像的分辨率非常低，传统处理估计得到的图像配准偏移量及补偿相位估计精度难以满足后续处理需求时，本发明可以通过对主辅SAR图像对进行有效的图像融合以提高方位向分辨率。

背景技术

SAR图像方位分辨率与雷达观测目标的合成孔径长度(对应方位角带宽)密切相关，它随合成孔径长度的增大而提高。SAR图像融合技术利用多幅具有不同视角范围的SAR图像融合处理得到一幅具有更高分辨率的SAR图像，该技术的处理性能与各幅SAR图像获取时间间隔内的地物变化(即时间去相干)及方位角范围密切相关。分布式卫星SAR系统单次航过即可获取同一地区的多幅SAR图像，这些SAR图形间的时间去相干可忽略不计。

此外，通过有效调节卫星编队构型并利用波束同步技术，分布式卫星SAR系统在沿航向长基线构型下可获取方位视角重叠度近似为零的多幅SAR图像，此时利用图像融合技术科有效提升SAR图像方位分辨率。因此该技术可视为高分辨率宽测绘带(High-Resolution and Wide-Swath,HRWS)SAR系统对地观测的一种有效实现。Prats P.,Scheiber R.and Mittermayer J.在文献“TanDEM-X experiments in pursuit monostatic configuration”(9th EuropeanConference on Synthetic Aperture Radar(EUSAR),2012:159-162.)中提出了一种分布式卫星SAR图像融合提升方位分辨率的方法。该方法的具体步骤为：1)对主辅SAR图像多普勒包络进行校正；2)通过方位向预滤波截取具有公共多普勒的主辅SAR图像；3)利用雷达几何和图像相关法对方位预滤波后的主辅SAR图像进行配准；4)对配准后的主辅SAR图像进行干涉相位提取及相位偏差校正；5)进行多普勒谱拼接。该方法的性能依赖于公共多普勒谱分量的大小。当沿航向长基线分布式卫星SAR系统的多普勒谱重叠度近似为零时，方位预滤波后SAR图像的分辨率非常低时，此时传统处理估计得到的图像配准偏移量及补偿相位估计精度难以满足后续融合处理要求，从而造成融合后SAR图像散焦(分辨率损失、旁瓣抬高等问题)。

发明内容

本发明的目的在于提出一种分布式卫星SAR系统的主辅SAR图像融合方法。本发明首先利用外部先验DEM信息及主辅雷达成像几何反演观测场景的地形干涉相位，然后利用最大对比度准则估计主辅SAR图像配准误差及残留干涉相位偏差，最后利用上述估计值对主辅SAR图像进行相干相加得到一高分辨率的SAR图像。

为实现上述技术目的，本发明采用如下技术方案予以实现。

一种分布式卫星SAR系统的主辅SAR图像融合方法包括以下步骤：

步骤1，利用分布式卫星SAR系统获取主SAR图像数据和辅SAR图像数据；

步骤2，分别对主SAR图像数据和辅SAR图像数据进行多普勒谱包络校正，得出多普勒谱包络校正后的主SAR图像数据和多普勒谱包络校正后的辅SAR图像数据；

步骤3，对多普勒谱包络校正后的辅SAR图像数据进行雷达几何配准，得出几何配准后的辅SAR图像数据；

步骤4，得出多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的地形相位估计值；

步骤5，在时域将几何配准后的辅SAR图像数据和多普勒谱包络校正后的主SAR图像数据进行相干叠加，获得初步融合后的SAR图像；以初步融合后的SAR图像的对比度最大为目标，建立优化模型，通过求解所述优化模型得出方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res；

步骤6，针对多普勒谱包络校正后的辅SAR图像数据，按照步骤5得出的方位向配准误差和距离向配准误差进行精配准，得出精配准后的辅SAR图像数据；根据步骤5得出的残留相位偏差对精配准后的辅SAR图像数据进行相位偏差校正，得到相位偏差校正后的辅SAR图像数据；

步骤7，分别对多普勒谱包络校正后的主SAR图像数据和相位偏差校正后的辅SAR图像数据进行过采样，得出过采样后的主SAR图像数据和过采样后的辅SAR图像数据；分别对过采样后的主SAR图像数据和过采样后的辅SAR图像数据进行傅里叶变换，得到主SAR图像多普勒域数据和辅SAR图像多普勒域数据；对主SAR图像多普勒域数据和辅SAR图像多普勒域数据在多普勒域进行多普勒谱拼接，得到多普勒谱拼接后的SAR图像多普勒域数据；将多普勒谱拼接后的SAR图像多普勒域数据进行傅里叶逆变换，得到融合后的SAR图像数据。

本发明的有益效果为：本发明首先利用先验DEM及主辅SAR图像雷达成像几何模拟计算地形干涉相位，然后利用融合后SAR图像最大对比度估计残留相位偏差及配准偏移量，克服了现有技术分布式卫星SAR数据处理方法中的提升方位分辨率方法在主辅SAR图像多普勒谱重叠近零时，处理性能不能满足工程应用要求的缺陷，使得分布式SAR产品可以达到高精度测绘产品的要求。

附图说明

图1为本发明的一种分布式卫星SAR系统的主辅SAR图像融合方法的流程图；

图2a为实测数据实验中观测区域的主SAR图像；

图2b为实测数据实验中观测区域的辅SAR图像；

图2c为实测数据实验中采用传统谱拼接提升方位分辨率算法(参见“LiuY.,Li Z.and Suo Z.et.al..Azimuth resolution improvement of spaceborne SARimages with nearly non-overlapped Doppler bandwidth.IET International RadarConference,2013:1-4.”)融合SAR图像；

图2d为实测数据实验中采用本发明得出的融合SAR图像；

图3为实测数据实验得出的角反射器的方位向包络示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的一种分布式卫星SAR系统的主辅SAR图像融合方法的流程图，其具体步骤为：

步骤1，利用分布式卫星SAR系统获取主SAR图像数据(主SAR图像数据通过分布式卫星合成孔径雷达的主合成孔径雷达获取)和辅SAR图像数据(辅SAR图像数据通过分布式卫星合成孔径雷达的辅合成孔径雷达获取)。主SAR图像数据和辅SAR图像数据要满足干涉处理的在相干性和成像质量等方面的质量要求。

设置SAR成像处理参数，本发明实施例中，SAR成像处理参数包括：原始回波信号的多普勒中心频率、最近斜距、斜距上的采样间隔、方位向采样间隔等。

步骤2，分别对主SAR图像数据和辅SAR图像数据进行多普勒谱包络校正，得出多普勒谱包络校正后的主SAR图像数据和多普勒谱包络校正后的辅SAR图像数据。

其具体子步骤为：

(2.1)根据分布式卫星SAR系统的主合成孔径雷达的天线方向图、分布式卫星SAR系统的主合成孔径雷达的波束指向(发射波束指向)、以及设定的成像处理窗函数，得出主SAR图像数据的多普勒中心频率f_dc1。根据分布式卫星SAR系统的辅合成孔径雷达的天线方向图、分布式卫星SAR系统的辅合成孔径雷达的波束指向(发射波束指向)、以及设定的成像处理窗函数，得出辅SAR图像数据的多普勒中心频率f_dc2。本发明中，得出主SAR图像数据的多普勒中心频率和辅SAR图像数据的多普勒中心频率的方法为现有方法。

然后，得出主SAR图像数据的多普勒包络校正函数H₁和辅SAR图像数据的多普勒包络校正函数H₂，H₁＝exp(j2π*f_a*f_dc1)，H₂＝exp(j2π*f_a*f_dc2)，其中，f_a表示方位频率。

(2.2)将主SAR图像数据的多普勒包络校正函数在多普勒域与主SAR图像数据相乘，进行主SAR图像多普勒谱包络校正，得出多普勒谱包络校正后的主SAR图像数据。对辅SAR图像数据的多普勒包络校正函数在多普勒域与辅SAR图像数据相乘，进行辅SAR图像多普勒谱包络校正，得出多普勒谱包络校正后的辅SAR图像数据。在进行主辅SAR图像多普勒谱包络校正之后，获得了多普勒带宽内包络平坦的主辅SAR图像数据。

步骤3，对多普勒谱包络校正后的辅SAR图像数据进行雷达几何配准，得出几何配准后的辅SAR图像数据。

其具体子步骤为：

(3.1)得出多普勒谱包络校正后的主SAR图像数据中每个数据点的雷达斜距，多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点的雷达斜距表示为r_M(m,n)。得出多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的主天线相位中心位置、以及多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的速度矢量；多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点的方位时刻表示为t_M(m)，多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_M(m)对应的主天线相位中心位置表示为p_M(m,n；t_M)，多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_M(m)对应的速度矢量表示为v_M(m,n；t_M)。

(3.2)根据多普勒谱包络校正后的主SAR图像数据中每个数据点的雷达斜距、多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的主天线相位中心位置、以及多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的速度矢量，并基于先验数字高程模型(DigitalElevation Model,简称DEM)，进行目标点定位，得出多普勒谱包络校正后的主SAR图像数据中每个数据点对应的目标高程。

在子步骤(3.2)中，得出多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标高程的过程包括如下子步骤：

(3.2.1)预设目标高程为h₀；令h＝h₀

(3.2.2)求解SAR几何定位方程，得到多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标点空间位置p_t：

p_t＝(p_t,x,p_t,y,p_t,z)^T

其中，p_t,x表示对应数据点在WGS-84坐标系下(World Geadetic system-1989Coordinate System)的x轴坐标，p_t,y表示对应数据点在WGS-84坐标系下的y轴坐标，p_t,z表示对应数据点在WGS-84坐标系下的z轴坐标，上标T表示矩阵或向量的转置。

在子步骤(3.2.2)中，SAR几何定位方程为：

$\left\{\begin{array}{c}\frac{2v_M(m,n;t_M)\·[p_t-p_M(m,n;t_M)]}{\mathrm{\λ}{r}_M(m,n)}=0\\ |p_t-p_M(m,n;t_M)|=r_M(m,n)\\ \frac{p_{t,x}^2+p_{t,y}^2}{{(R_e+h)}^2}+\frac{p_{t,z}^2}{{(1-f)}^2{(R_e+h)}^2}=1\end{array}\right.$

在上述SAR几何定位方程中，λ为分布式卫星SAR系统的主合成孔径雷达的发射信号波长，|·|表示取模值，R_e为赤道半径，f为地球扁率因子，h表示目标高程。

(3.2.3)将目标点空间位置p_t转换至先验数字高程模型所在坐标系，得出坐标转换后的目标点空间位置。

(3.2.4)在先验数字高程模型内插值出该目标点水平位置对应的高程h₁；

(3.2.5)如果|h₁-h₀|<δ，则将h₁作为多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标高程，然后转至子步骤(3.3)，|h₁-h₀|表示h₁和h₀的差的绝对值，δ为预设门限；如果|h₁-h₀|≥δ，则返回至子步骤(3.2.1)。

(3.3)根据多普勒谱包络校正后的主SAR图像数据中每个数据点对应的目标高程，将对应的数据点反定位到多普勒谱包络校正后的辅SAR图像数据上，得出多普勒谱包络校正后的辅SAR图像数据中数据点对应的方位时刻、以及多普勒谱包络校正后的辅SAR图像数据中数据点对应的雷达斜距。本发明实施例中，多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的方位时刻表示为t_S(m)，多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的雷达斜距表示为r_S(m,n)。

具体地说，根据多普勒谱包络校正后的主SAR图像数据中每个数据点在地球椭球模型下的坐标、以及分布式卫星SAR系统的辅合成孔径雷达记录的辅合成孔径雷达的天线相位中心坐标参数，并采用多普勒中心频率比较法得到多普勒谱包络校正后的辅SAR图像数据中对应数据点对应的方位时刻；再根据SAR辅图像成像几何关系，得出多普勒谱包络校正后的辅SAR图像数据中对应数据点对应的雷达斜距。

(3.4)计算多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的斜距点、以及多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的方位点。

多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的斜距点i'和多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的方位点j的计算公式为：

t_S(m)＝t_{S_start}+(j-1)/PRF

r_S(m,n)＝r_{S_near}+(i'-1)*r_{S_sample}

其中，t_{S_start}为多普勒谱包络校正后的辅SAR图像数据的方位向起始时刻，PRF为多普勒谱包络校正后的辅SAR图像数据的方位向采样率，r_{S_near}为多普勒谱包络校正后的辅SAR图像数据中所有数据点的斜距的最近斜距，r_{S_sample}为多普勒谱包络校正后的辅SAR图像数据的斜距采样间隔，t_{S_start}、PRF、r_{S_near}、r_{S_sample}为幅图像数据已知参数。

对多普勒谱包络校正后的辅SAR图像数据按照斜距点和方位点进行插值，将多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据表示为S'_s(m,n)，将多普勒谱包络校正后的辅SAR图像数据在斜距点i'和方位点j处的插值数据表示为S_s(j,i')。在多普勒谱包络校正后的辅SAR图像数据中，将每个数据点的数据值替换为对应斜距点和对应方位点的插值数据，从而得出几何配准后的辅SAR图像数据。其中，在多普勒谱包络校正后的辅SAR图像数据中，将方位向坐标为m距离向坐标为n的数据点的数据S'_s(m,n)替换为S_s(j,i')。

步骤4，得出多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的地形相位估计值。

其具体子步骤为：

(4.1)求取多普勒谱包络校正后的主SAR图像数据中每个数据点在主SAR图像成像几何下的回波延时、以及多普勒谱包络校正后的辅SAR图像数据中每个数据点在辅SAR图像成像几何下的回波延时。

多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在主SAR图像成像几何下的回波延时τ₁、以及多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点在辅SAR图像成像几何下的回波延时τ₂分别为：

τ₁＝2|p_M(m,n；t_M)-p_t|/c

τ₂＝(|p_M(m,n；t_M)-p_t|+|p_S(m,n；t_S)-p_t|)/c

其中，|·|表示取模值，p_S(m,n；t_S)表示多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_S(m)对应的辅天线相位中心位置，p_M(m,n；t_M)表示多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_M(m)对应的主天线相位中心位置，c为光速。

(4.2)计算多普勒谱包络校正后的主/辅SAR图像数据中每个数据点对应的地形相位估计值，多普勒谱包络校正后的主/辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的地形相位估计值为：

${\mathrm{\&phi;}}_{\mathrm{topo}}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}c({\mathrm{\&tau;}}_1-{\mathrm{\&tau;}}_2)$

其中，λ表示分布式卫星SAR系统的主合成孔径雷达的发射信号波长，c表示光速。

步骤5，在时域将步骤3几何配准后的辅SAR图像数据和步骤2得出的多普勒谱包络校正后的主SAR图像数据进行相干叠加，获得初步融合后的SAR图像；)以初步融合后的SAR图像的对比度最大为目标，建立优化模型，通过求解所述优化模型得出方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res。

其具体子步骤为：

(5.1)初始阶段，将(Δt_a,Δτ,φ_res)的初始值设定为(0,0,)，其中，Δt_a为方位向配准误差、Δτ为距离向配准误差，φ_res为残留相位偏差。

(5.2)在时域将步骤3几何配准后的辅SAR图像数据和步骤2得出的多普勒谱包络校正后的主SAR图像数据进行相干叠加，获得初步融合后的SAR图像。

(5.3)以初步融合后的SAR图像的对比度最大为目标，建立优化模型，通过求解所述优化模型得出方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res。

子步骤(5.3)中，优化模型为：

$[{\mathrm{\&Delta;}\hat{t}}_a,\mathrm{\&Delta;}\hat{t},{\widehat{\mathrm{\&phi;}}}_{\mathrm{res}}]=\arg \underset{\mathrm{\&Delta;t},\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}}}{\max }F(\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})$

$F({\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})=\underset{(m,n)\&Element;A}{\mathrm{\&Sigma;}}{\left|S_{\mathrm{syn}}(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})\right|}^4$

S_syn(m,n；Δt_a,Δτ,φ_res)＝S'_m(m,n)+S'_s(m,n；Δt_a,Δτ,φ_res)

其中，F(Δt_a,Δτ,φ_res)为目标函数，F(Δt_a,Δτ,φ_res)表示初步融合后的SAR图像的对比度；|S_syn(m,n；Δt_a,Δτ,φ_res)|为S_syn(m,n；Δt_a,Δτ,φ_res)的模值，A表示主/辅SAR图像数据中所有数据点的坐标位置(包括方位向坐标和距离向坐标)的集合；S'_m(m,n)表示多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据，S'_s(m,n；Δt_a,Δτ,φ_res)表示多普勒谱包络校正后方位向配准误差为Δt_a距离向配准误差为Δτ且残留相位偏差为φ_res的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据。这里需要指出的是，当雷达斜视角较小时，距离向配准误差Δτ引起的相位误差可忽略不计，但距离向配准误差将使主辅SAR图像包络失配，降低合成后SAR图像的对比度。因此，通过上述目标函数可较好地估计距离向配准误差。

在求解方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res时，选择方位向配准误差Δt_a的多种不同样本值、距离向配准误差Δτ的多种不同样本值以及残留相位偏差φ_res的多种不同样本值，来计算初步融合后的SAR图像的对比度F(Δt_a,Δτ,φ_res)，最终确定使F(Δt_a,Δτ,φ_res)取最大值时的方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res。本发明实施例中，方位向配准误差Δt_a的样本选取步长为0.1个SAR图像采样间隔，距离向配准误差Δτ的样本选取步长为0.1个SAR图像采样间隔，残留相位偏差φ_res的样本选取步长可利用牛顿迭代方法确定：

${\mathrm{\&Delta;\&phi;}}_{\mathrm{res},i}=\frac{{\frac{\&PartialD;F({\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})}{{\&PartialD;\mathrm{\&phi;}}_{\mathrm{res}}}|}_{{\mathrm{\&phi;}}_{\mathrm{res}}={\mathrm{\&phi;}}_{\mathrm{res},i-1}}}{{\frac{{\&PartialD;}^{2}F({\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})}{{\&PartialD;\mathrm{\&phi;}}_{\mathrm{res}}^2}|}_{{\mathrm{\&phi;}}_{\mathrm{res}}={\mathrm{\&phi;}}_{\mathrm{res},i-1}}}$

其中，

$\begin{array}{c}{\frac{\&PartialD;F(\mathrm{\&Delta;t},\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})}{{\&PartialD;\mathrm{\&phi;}}_{\mathrm{res}}}|}_{{\mathrm{\&phi;}}_{\mathrm{res}}={\mathrm{\&phi;}}_{\mathrm{res},i-1}}\\ =4\underset{(m,n)\&Element;A}{\mathrm{\&Sigma;}}{\left|S_{\mathrm{syn}}(m,n;{\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res},i-1})\right|}^2\&times;\\ \left\{\mathrm{Re}\right[S_1(m,n)\&CenterDot;S_2(m,n;{\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res},i-1})]\sin (-{\mathrm{\&phi;}}_{\mathrm{res},i-1})+\\ \mathrm{Im}[S_1(m,n)\&CenterDot;S_2(m,n;{\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res},i-1})]\cos {\mathrm{\&phi;}}_{\mathrm{res},i-1}\}\end{array}$

$\begin{array}{c}\frac{{\&PartialD;}^{2}F(\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})}{\&PartialD;{\mathrm{\&phi;}}_{\mathrm{res}}^2}{|}_{{\mathrm{\&phi;}}_{\mathrm{res}}={\mathrm{\&phi;}}_{\mathrm{res},i-1}}\\ =8\underset{(m,n)\&Element;A}{\mathrm{\&Sigma;}}\left\{\mathrm{Re}\right[S_1(m,n)\&CenterDot;S_2(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{tes},i-1})]\sin (-{\mathrm{\&phi;}}_{\mathrm{res},i-1})\\ +\mathrm{Im}[S_1(m,n)\&CenterDot;S_2(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res},i-1})]\cos {\mathrm{\&phi;}}_{\mathrm{res},i-1}{\}}^2\\ -4\underset{(m,n)\&Element;A}{\mathrm{\&Sigma;}}{\left|S_{\mathrm{syn}}(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{tes},i-1})\right|}^2\&CenterDot;\left\{\mathrm{Re}\right[S_1(m,n)\&CenterDot;S_2(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res},i-1})]\\ \&times;\cos {\mathrm{\&phi;}}_{\mathrm{res},i-1}+\mathrm{Im}[S_1(m,n)\&CenterDot;S_2(m,n;\mathrm{\&Delta;}{t}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{tes},i-1})]\sin {\mathrm{\&phi;}}_{\mathrm{res},i-1}\}\end{array}$

其中，Δφ_res,i为残留相位偏差φ_res的第i次步长，φ_res,i为残留相位偏差φ_res的第i次取值，Re[·]表示取实部，Im[·]表示取虚部；S₁(m,n)为多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据，S₂(m,n；Δt_a,Δτ,φ_res,i-1)表示多普勒谱包络校正后方位向配准误差为Δt_a距离向配准误差为Δτ且残留相位偏差为φ_res,i-1的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据。

步骤6，针对多普勒谱包络校正后的辅SAR图像数据，按照步骤5得出的方位向配准误差和距离向配准误差进行精配准，得出精配准后的辅SAR图像数据；根据步骤5得出的残留相位偏差对精配准后的辅SAR图像数据进行相位偏差校正，得到相位偏差校正后的辅SAR图像数据。

其具体子步骤为：

(6.1)针对多普勒谱包络校正后的辅SAR图像数据中每个数据点的数据，根据步骤5得出的对应的方位向配准误差和距离向配准误差进行插值，得出对应的插值数据；然后将多普勒谱包络校正后的辅SAR图像数据中每个数据点的数据替换为对应的插值数据，得出精配准后的辅SAR图像数据。其中，对于多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据，根据方位向配准误差Δt_a和距离向配准误差Δτ进行插值，得出对应的插值数据S'_s(m+Δt_a,n+Δτ)；然后，将多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据S'_s(m,n)替换为用S'_s(m+Δt_a,n+Δτ)。

(6.2)将精配准后的辅SAR图像数据中每个数据点的数据与步骤5得出的对应的残留相位偏差相乘，得到相位偏差校正后的辅SAR图像数据。其中，对于精配准后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据S'_s(m+Δt_a,n+Δτ)，将其与步骤5得出的对应的残留相位偏差φ_res进行相乘，将相乘后的结果替换S'_s(m+Δt_a,n+Δτ)。

步骤7，分别对多普勒谱包络校正后的主SAR图像数据和相位偏差校正后的辅SAR图像数据进行过采样，得出过采样后的主SAR图像数据和过采样后的辅SAR图像数据；分别对过采样后的主SAR图像数据和过采样后的辅SAR图像数据进行傅里叶变换，得到主SAR图像多普勒域数据和辅SAR图像多普勒域数据；对主SAR图像多普勒域数据和辅SAR图像多普勒域数据在多普勒域进行多普勒谱拼接，得到多普勒谱拼接后的SAR图像多普勒域数据；将多普勒谱拼接后的SAR图像多普勒域数据进行傅里叶逆变换，得到融合后的SAR图像数据。

其具体子步骤为：

(7.1)对步骤2得出的多普勒谱包络校正后的主SAR图像数据进行两倍过采样，得出过采样后的主SAR图像数据，对过采样后的主SAR图像数据进行傅里叶变换，将其变换至多普勒域，得到主SAR图像多普勒域数据。

对相位偏差校正后的辅SAR图像数据进行两倍过采样，得出过采样后的辅SAR图像数据，对过采样后的辅SAR图像数据进行傅里叶变换，将其变换至多普勒域，得到辅SAR图像多普勒域数据。

(7.2)对主SAR图像多普勒域数据和辅SAR图像多普勒域数据在多普勒域进行多普勒谱拼接，得到多普勒谱拼接后的SAR图像多普勒域数据。

多普勒谱拼接的过程可以用以下公式进行表示：

S_c(f_d)＝H₁(f_d)S₁(f_d)+H₂(f_d)S'₂(f_d)

式中，f_d为多普勒频率，S_c(f_d)表示多普勒谱拼接后的SAR图像多普勒域数据；S₁(f_d)表示主SAR图像多普勒域数据，S'₂(f_d)表示辅SAR图像多普勒域数据，H₁(f_d)为主SAR图像多普勒域数据的多普勒谱拼接函数，H₂(f_d)表示辅SAR图像多普勒域数据的多普勒谱拼接函数。

本发明实施例中，f_dc,1表示主SAR图像数据的多普勒中心频率，f_dc,2表示主SAR图像数据的多普勒中心频率；当f_dc,1<f_dc,2时，

当f_dc,1≥f_dc,2时，

其中，B_a,1表示主SAR图像多普勒域数据的多普勒带宽，B_a,2表示辅SAR图像多普勒域数据的多普勒带宽，B_a,c为公共多普勒带宽，B_a,c＝(f_dc,1+B_a,1/2)-(f_dc,2-B_a,2/2)。

(7.3)将多普勒谱拼接后的SAR图像多普勒域数据进行傅里叶逆变换，将其转换至时域，得到融合后的SAR图像数据。

此外，作为本发明实施例的一种变形，在子步骤(7.2)中，在得出多普勒谱拼接后的SAR图像多普勒域数据S_c(f_d)之后，还可以对多普勒谱拼接后的SAR图像多普勒域数据S_c(f_d)进行加窗处理，得到加窗处理后的SAR图像多普勒域数据，加窗处理可以抑制SAR图像旁瓣。在子步骤(7.3)中，对加窗处理后的SAR图像多普勒域数据进行傅里叶逆变换，将其转换至时域，得到融合后的SAR图像数据。

下面结合实测数据实验对本发明的效果做进一步的说明。

这里利用德国TerraSAR-X系统重复航过获取的两幅SAR图像模拟仿真多普勒带宽近零重叠的两幅沿航向长基线分布式卫星SAR图像。在实测数据实验中，采用两幅TerraSAR-X实测SAR图像(分别于2008年2月25日和2008年3月8日获取)，地面观测区域为我国兰州地区，其原始回波多普勒中心分别约为-54.4Hz和4.1Hz，成像处理采用的多普勒处理带宽均为2765Hz，垂直有效基线为6.851m，对应的系统模糊高度为602.8m。为模拟仿真多普勒带宽近零重叠的两幅沿航向长基线分布式卫星SAR图像，这里将两幅SAR图像进行多普勒谱截取，带宽为1331Hz，重叠带宽为14.6Hz。这对应于沿航向基线长度约为1533.88m。参见表1，为实测数据实验使用的TerraSAR-X系统参数。

TerraSAR-X 系统参数

参照图2a，为实测数据实验中观测区域的主SAR图像，参照图2b，为实测数据实验中观测区域的辅SAR图像。图2a和图2b中，圆圈标注处为角反射器(对应主/辅SAR图像中的数据点)。在实测数据实验中，分别采用传统谱拼接提升方位分辨率算法和本发明进行主辅SAR图像融合，参照图2c，为实测数据实验中采用传统谱拼接提升方位分辨率算法得出的融合SAR图像；参照图2d，为实测数据实验中采用本发明得出的融合SAR图像。在图2a至图2d中，横轴表示距离向坐标，单位为点数，纵轴表示方位向坐标，单位为点数。由图2c和图2d的对比可见，与采用传统谱拼接提升方位分辨率算法得出的融合SAR图像相比，采用本发明得出的融合SAR图像在方位向(纵向)更清晰，对比度更高，说明与采用传统谱拼接提升方位分辨率算法得出的融合SAR图像相比，采用本发明得出的融合SAR图像提升了SAR图像的方位向分辨率。

参照图3，为实测数据实验得出的角反射器的方位向包络示意图。图中三条线分别表示主SAR图像的角反射器的方位向包络示意图(对应图3中的SAR主图像)、采用传统谱拼接提升方位分辨率算法得出的融合SAR图像的角反射器的方位向包络示意图(对应图3中的传统方法处理结果)、以及采用本发明得出的融合SAR图像的角反射器的方位向包络示意图(对应图3中的新算法处理结果)。图3中，横轴表示距离向坐标，单位为m，纵轴表示强度，单位为dB。在实测数据实验中，主SAR图像的方位向分辨率为4.74m，峰值旁瓣比分别为-12.58dB；采用传统谱拼接提升方位分辨率算法得出的融合SAR图像的方位向分辨率为2.56m，峰值旁瓣比分别为-4.15dB；采用本发明得出的融合SAR图像的方位向分辨率为2.40m，峰值旁瓣比分别为-11.47dB。可以看出，在实测数据实验中(主辅SAR图像多普勒带宽近零重叠)，本发明可以将主SAR图像的方位向分辨率提升接近一倍。本发明的SAR图像融合方法的方位向分辨率和峰值旁瓣比的性能均高于传统谱拼接提升方位分辨率算法。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (7)

1.一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，包括以下步骤：

步骤1，利用分布式卫星SAR系统获取主SAR图像数据和辅SAR图像数据；

步骤2，分别对主SAR图像数据和辅SAR图像数据进行多普勒谱包络校正，得出多普勒谱包络校正后的主SAR图像数据和多普勒谱包络校正后的辅SAR图像数据；

步骤3，对多普勒谱包络校正后的辅SAR图像数据进行雷达几何配准，得出几何配准后的辅SAR图像数据；

步骤4，得出多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的地形相位估计值；

步骤5，在时域将几何配准后的辅SAR图像数据和多普勒谱包络校正后的主SAR图像数据进行相干叠加，获得初步融合后的SAR图像；以初步融合后的SAR图像的对比度最大为目标，建立优化模型，通过求解所述优化模型得出方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res；

步骤6，针对多普勒谱包络校正后的辅SAR图像数据，按照步骤5得出的方位向配准误差和距离向配准误差进行精配准，得出精配准后的辅SAR图像数据；根据步骤5得出的残留相位偏差对精配准后的辅SAR图像数据进行相位偏差校正，得到相位偏差校正后的辅SAR图像数据；

步骤7，分别对多普勒谱包络校正后的主SAR图像数据和相位偏差校正后的辅SAR图像数据进行过采样，得出过采样后的主SAR图像数据和过采样后的辅SAR图像数据；分别对过采样后的主SAR图像数据和过采样后的辅SAR图像数据进行傅里叶变换，得到主SAR图像多普勒域数据和辅SAR图像多普勒域数据；对主SAR图像多普勒域数据和辅SAR图像多普勒域数据在多普勒域进行多普勒谱拼接，得到多普勒谱拼接后的SAR图像多普勒域数据；将多普勒谱拼接后的SAR图像多普勒域数据进行傅里叶逆变换，得到融合后的SAR图像数据。

2.如权利要求1所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，所述步骤2的具体子步骤为：

(2.1)根据分布式卫星SAR系统的主合成孔径雷达的天线方向图、分布式卫星SAR系统的主合成孔径雷达的波束指向、以及设定的成像处理窗函数，得出主SAR图像数据的多普勒中心频率f_dc1，根据分布式卫星SAR系统的辅合成孔径雷达的天线方向图、分布式卫星SAR系统的辅合成孔径雷达的波束指向、以及设定的成像处理窗函数，得出辅SAR图像数据的多普勒中心频率f_dc2；得出主SAR图像数据的多普勒包络校正函数H₁和辅SAR图像数据的多普勒包络校正函数H₂，H₁＝exp(j2π*f_a*f_dc1)，H₂＝exp(j2π*f_a*f_dc2)，其中，f_a表示方位频率；

(2.2)将主SAR图像数据的多普勒包络校正函数在多普勒域与主SAR图像数据相乘，得出多普勒谱包络校正后的主SAR图像数据；对辅SAR图像数据的多普勒包络校正函数在多普勒域与辅SAR图像数据相乘，得出多普勒谱包络校正后的辅SAR图像数据。

3.如权利要求1所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，所述步骤3的具体子步骤为：

(3.1)得出多普勒谱包络校正后的主SAR图像数据中每个数据点的雷达斜距，得出多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的主天线相位中心位置、以及多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的速度矢量；

(3.2)根据多普勒谱包络校正后的主SAR图像数据中每个数据点的雷达斜距、多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的主天线相位中心位置、以及多普勒谱包络校正后的主SAR图像数据中每个数据点在其方位时刻对应的速度矢量，基于先验数字高程模型进行目标点定位，得出多普勒谱包络校正后的主SAR图像数据中每个数据点对应的目标高程；

(3.3)根据多普勒谱包络校正后的主SAR图像数据中每个数据点对应的目标高程，将对应的数据点反定位到多普勒谱包络校正后的辅SAR图像数据上，得出多普勒谱包络校正后的辅SAR图像数据中数据点对应的方位时刻、以及多普勒谱包络校正后的辅SAR图像数据中数据点对应的雷达斜距；

(3.4)计算多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的斜距点、以及多普勒谱包络校正后的辅SAR图像数据中每个数据点对应的方位点；对多普勒谱包络校正后的辅SAR图像数据按照斜距点和方位点进行插值，在多普勒谱包络校正后的辅SAR图像数据中，将每个数据点的数据值替换为对应斜距点和对应方位点的插值数据，得出几何配准后的辅SAR图像数据。

4.如权利要求3所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，在子步骤(3.2)中，得出多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标高程的过程包括如下子步骤：

(3.2.1)预设目标高程为h₀；令h＝h₀

(3.2.2)求解SAR几何定位方程，得到多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标点空间位置p_t：

p_t＝(p_t,x,p_t,y,p_t,z)^T

其中，p_t,x表示对应数据点在WGS-84坐标系下的x轴坐标，p_t,y表示对应数据点在WGS-84坐标系下的y轴坐标，p_t,z表示对应数据点在WGS-84坐标系下的z轴坐标，上标T表示矩阵或向量的转置；

(3.2.3)将目标点空间位置p_t转换至先验数字高程模型所在坐标系，得出坐标转换后的目标点空间位置；

(3.2.4)在先验数字高程模型内插值出目标点空间位置对应的高程h₁；

(3.2.5)如果|h₁-h₀|<δ，则将h₁作为多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的目标高程，然后转至子步骤(3.3)，|h₁-h₀|表示h₁和h₀的差的绝对值，δ为预设门限；如果|h₁-h₀|≥δ，则返回至子步骤(3.2.1)。

5.如权利要求1所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)求取多普勒谱包络校正后的主SAR图像数据中每个数据点在主SAR图像成像几何下的回波延时、以及多普勒谱包络校正后的辅SAR图像数据中每个数据点在辅SAR图像成像几何下的回波延时；

多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在主SAR图像成像几何下的回波延时τ₁、以及多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点在辅SAR图像成像几何下的回波延时τ₂分别为：

τ₁＝2|p_M(m,n；t_M)-p_t|/c

τ₂＝(|p_M(m,n；t_M)-p_t|+|p_S(m,n；t_S)-p_t)|/c

其中，|·|表示取模值，p_S(m,n；t_S)表示多普勒谱包络校正后的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_S(m)对应的辅天线相位中心位置，p_M(m,n；t_M)表示多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点在方位时刻t_M(m)对应的主天线相位中心位置，c为光速；

(4.2)计算多普勒谱包络校正后的主/辅SAR图像数据中每个数据点对应的地形相位估计值，多普勒谱包络校正后的主/辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点对应的地形相位估计值为：

${\mathrm{\&phi;}}_{\mathrm{topo}}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}c({\mathrm{\&tau;}}_1-{\mathrm{\&tau;}}_2)$

其中，λ表示分布式卫星SAR系统的主合成孔径雷达的发射信号波长。

6.如权利要求1所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，所述步骤5的具体子步骤为：

(5.1)设置方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res的初始值；

(5.2)在时域将步骤3几何配准后的辅SAR图像数据和步骤2得出的多普勒谱包络校正后的主SAR图像数据进行相干叠加，获得初步融合后的SAR图像；

(5.3)以初步融合后的SAR图像的对比度最大为目标，建立优化模型，通过求解所述优化模型得出方位向配准误差Δt_a、距离向配准误差Δτ以及残留相位偏差φ_res；所述优化模型为：

$[\mathrm{\&Delta;}{\hat{t}}_a,\mathrm{\&Delta;}\widehat{\mathrm{\&tau;}},{\widehat{\mathrm{\&phi;}}}_{\mathrm{res}}]=\arg \underset{\mathrm{\&Delta;t},\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}}}{\max }F({\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})$

$F({\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})=\underset{(m,n)\&Element;A}{\mathrm{\&Sigma;}}{\left|S_{\mathrm{syn}}(m,n;{\mathrm{\&Delta;t}}_a,\mathrm{\&Delta;\&tau;},{\mathrm{\&phi;}}_{\mathrm{res}})\right|}^4$

S_syn(m,n；Δt_a,Δτ,φ_res)＝S′_m(m,n)+S′_s(m,n；Δt_a,Δτ,φ_res)

其中，F(Δt_a,Δτ,φ_res)表示初步融合后的SAR图像的对比度；|S_syn(m,n；Δt_a,Δτ,φ_res)|为S_syn(m,n；Δt_a,Δτ,φ_res)的模值，A表示主/辅SAR图像数据中所有数据点的坐标位置的集合；S′_m(m,n)表示多普勒谱包络校正后的主SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据，S′_s(m,n；Δt_a,Δτ,φ_res)表示多普勒谱包络校正后方位向配准误差为Δt_a距离向配准误差为Δτ且残留相位偏差为φ_res的辅SAR图像数据中方位向坐标为m距离向坐标为n的数据点的数据。

7.如权利要求1所述的一种分布式卫星SAR系统的主辅SAR图像融合方法，其特征在于，所述步骤7的具体子步骤为：

(7.1)对步骤2得出的多普勒谱包络校正后的主SAR图像数据进行两倍过采样，得出过采样后的主SAR图像数据，对过采样后的主SAR图像数据进行傅里叶变换，得到主SAR图像多普勒域数据；对相位偏差校正后的辅SAR图像数据进行两倍过采样，得出过采样后的辅SAR图像数据，对过采样后的辅SAR图像数据进行傅里叶变换，得到辅SAR图像多普勒域数据；

(7.2)对主SAR图像多普勒域数据和辅SAR图像多普勒域数据在多普勒域进行多普勒谱拼接，得到多普勒谱拼接后的SAR图像多普勒域数据；多普勒谱拼接后的SAR图像多普勒域数据S_c(f_d)为：

S_c(f_d)＝H₁(f_d)S₁(f_d)+H₂(f_d)S′₂(f_d)

式中，f_d为多普勒频率，S₁(f_d)表示主SAR图像多普勒域数据，S′₂(f_d)表示辅SAR图像多普勒域数据，H₁(f_d)为主SAR图像多普勒域数据的多普勒谱拼接函数，H₂(f_d)表示辅SAR图像多普勒域数据的多普勒谱拼接函数；

当主SAR图像数据的多普勒中心频率f_dc,1小于主SAR图像数据的多普勒中心频率f_dc,2时，

当主SAR图像数据的多普勒中心频率f_dc,1大于或等于主SAR图像数据的多普勒中心频率f_dc,2时，

其中，B_a,1表示主SAR图像多普勒域数据的多普勒带宽，B_a,2表示辅SAR图像多普勒域数据的多普勒带宽，B_a,c＝(f_dc,1+B_a,1/2)-(f_dc,2-B_a,2/2)；

(7.3)将多普勒谱拼接后的SAR图像多普勒域数据进行傅里叶逆变换，将其转换至时域，得到融合后的SAR图像数据。