本发明涉及遥感影像处理领域，更确切地说，它涉及一种面向二维海表流场sar反演的多普勒中心频率估计方法。
背景技术：
1、合成孔径雷达(synthetic aperture radar，sar)散射回波的多普勒频移具有反映海面动态特性的能力，已成为关键的海表动力遥感参数。在中小尺度海流监测和制图中，多普勒频移异常(doppler centroid anomaly)算法是主流方法之一。该算法的本质要求是从sar多普勒中心频率中精确分离出非海流运动导致的多普勒频移贡献。因此，准确估计sar多普勒中心频率对于提高海流sar反演的精度至关重要。
2、目前，仅有sentinel-1卫星和envisat卫星带宽模式的sar数据提供了多普勒中心参数，其他大部分sar数据都未提供这一关键参数，如高分三号、radarsat以及terrasar-x等。对于未提供多普勒中心频率的sar数据，需要先进行多普勒中心估计，才可以用于sar海流反演。多普勒中心频率估计方法包括两大类：一种是根据卫星轨道参数计算，另一种是从sar回波信号中估计。在实际应用中，卫星提供的轨道参数通常是不够精确的，因此，有必要从sar回波数据中估计多普勒中心频率，包括幅度和相位两类方法。
3、幅度方法指利用信号方位功率谱与相关函数之间的傅里叶变换对偶关系，在时域内对信号进行互相关系数计算，从而实现目标多普勒中心频率估计，如相关多普勒法(correlation doppler estimation，cde)和符号多普勒法(sign doppler estimation，sde)。与幅度方法不同的是，相位方法的核心是基于雷达波束的多普勒频谱在方位向上的对称原理，利用不同的加权函数与多普勒功率谱进行卷积运算，从卷积结果中获得多普勒中心频率。如能量均衡法(energy balance，eb)、相关匹配法(match correlation，mc)以及最优估计法(optimal estimation，op)。需要说明的是，采用不同的多普勒中心频率估计方法，其估计结果的平均值和成像质量存在差异。此外，不同多普勒中心频率估计算法对于不同扫描模式和空间分辨率的sar影像成像质量也具有差异。
技术实现思路
1、本发明的目的是针对现有技术的不足，提出了一种面向二维海表流场sar反演的多普勒中心频率估计方法。
2、第一方面，提供了一种面向二维海表流场sar反演的多普勒中心频率估计方法，包括：
3、步骤1、获取sar单视复数影像数据，并对所述sar单视复数影像数据进行分块处理；
4、步骤2、利用幅度方法获取多普勒中心频率估计值；所述幅度方法包括相关多普勒法和符号多普勒法；
5、步骤3、利用相位方法获取多普勒中心频率估计值；所述相位方法包括能量均衡法、匹配相关法和最优估计法；
6、步骤4、计算不同方法所获取多普勒中心频率估计值的信噪比，筛选出成像质量最优的多普勒中心频率估计值作为初估计结果，记为
7、步骤5、基于极大似然原理，对初估计结果进行迭代运算；
8、步骤6、去除卫星姿态引起的预测多普勒频移fgeo，获得多普勒频移异常fdca，再减去风浪偏差贡献的多普勒中心频率，完成径向海流sar反演；
9、步骤7、基于径向海流sar反演结果，利用二维ekman流反演方法，实现二维实时海流矢量反演。
10、作为优选，步骤2包括：
11、步骤2.1、利用相关多普勒法进行多普勒中心频率估计，表示为：
12、
13、
14、其中，x(m)＝x′(mt)定义为两个随机过程，即间隔为η的两个相邻子图像块，其相关系数记为n是每块图像方位向长度，x*(m)是x(m)的共轭复数，η为延迟；arg为相位角函数，t为两个相邻采样信号之间的时间间隔，等于1/prf；prf为脉冲频率；
15、步骤2.2、利用符号多普勒法进行多普勒中心频率估计；
16、步骤2.3、利用滑动窗口的方式，从左往右和从上往下遍历每一个图像块。
17、作为优选，步骤2.2包括：
18、步骤2.2.1、定义符号函数sv，并获取单视复数影像的符号相关性表达式，表示为：
19、
20、
21、其中，ny和nx分别为每块图像的方位向和距离向长度，sx(i+k,j)和sx(i,j)为间隔为k的两个局部窗口/样本；
22、步骤2.2.2、推导出归一化相关系数ρxy(η)和单视复数影像的复相关系数表示为：
23、
24、
25、其中，是实部i和虚部q的相关系数，和rh(k)的幅角相同；
26、步骤2.2.3、将复相关系数带入中，得到多普勒中心频率
27、作为优选，步骤3包括：
28、步骤3.1、对每一块图像进行傅里叶变换，将时域图像转换到频率域图像，求方位向多普勒频谱；
29、步骤3.2、选择不同的加权函数，与方位向多普勒频谱进行卷积运算；
30、步骤3.3、通过积分的方式从一个脉冲发射频率中，搜索卷积运算结果的能量峰值或过零点，其对应的多普勒频率即为多普勒中心频率，每个子图像返回一个多普勒中心频率估计值fdc_block；能量均衡法、匹配相关法和最优估计法的估计结果分别记为和
31、作为优选，步骤3.2包括：
32、步骤3.2.1、选择能量均衡、相关匹配以及最优估计三种加权函数，分别表示为b1(f)，b2(f)和b3(f)，公式为；
33、
34、
35、
36、其中，p(f)为方位向多普勒频谱，e(p(f))则为信号的功率谱密度；
37、步骤3.2.2、利用滑动窗口的方式，从左往右和从上往下遍历每一个图像块，将加权函数与方位向频谱进行卷积运算。
38、作为优选，步骤3.3中，用于搜索卷积运算结果过零点的积分方程为：
39、
40、其中，f(φ)为搜索过零点的积分函数，b(f)对应步骤3.2.1的三个加权函数；由此获得能量均衡法、匹配相关法和最优估计法三种相位方法的估计结果，分别记为和
41、作为优选，步骤4中，信噪比计算公式为：
42、
43、
44、其中，s0和s1分别为sar信号在方位向上的主峰和副瓣；通过对比五种多普勒中心频率估计结果的信噪比，筛选信噪比最优的估计结果作为迭代计算前的初估计值，记为
45、作为优选，步骤5中，迭代计算公式为：
46、
47、其中，fa为真实多普勒中心，f′(φ)为迭代计算中用于搜索方位向频谱过零点积分函数，fdc为迭代计算后的多普勒中心频率，s2/σ2为方差；过零点对应的多普勒中心频率等于
48、作为优选，步骤6包括：
49、步骤6.1、计算卫星姿态贡献的多普勒频移fgeo：
50、fgeo＝d0+d1(ts-t0)+d2(ts-t0)2+d3(ts-t0)3+d4(ts-t0)4
51、fdca＝fdc-fgeo
52、其中，di(i＝0,1,2,3,4)为多普勒系数，ts和t0分别为斜距时间和标准斜距时间；
53、步骤6.2、利用cdop地球物理模式函数估计风浪偏差贡献的多普勒频移fww，并将其从fdca中减去，获得由海流贡献的多普勒频移fosc：
54、
55、其中，u10和分别为风速和相对风向，θ和pol分别为雷达入射角和极化方式；
56、步骤6.3、计算径向流速u：
57、
58、式中，kr为雷达电磁波波数。
59、第二方面，提供了一种面向二维海表流场sar反演的多普勒中心频率估计系统，用于执行第一方面任一所述的方法，包括：
60、第一获取模块，用于获取sar单视复数影像数据，并对所述sar单视复数影像数据进行分块处理；
61、第二获取模块，用于利用幅度方法获取多普勒中心频率估计值；所述幅度方法包括相关多普勒法和符号多普勒法；
62、第三获取模块，用于利用相位方法获取多普勒中心频率估计值；所述相位方法包括能量均衡法、匹配相关法和最优估计法；
63、计算模块，用于计算不同方法所获取多普勒中心频率估计值的信噪比，筛选出成像质量最优的多普勒中心频率估计值作为初估计结果，记为
64、迭代模块，用于基于极大似然原理，对初估计结果进行迭代运算；
65、第一反演模块，用于去除卫星姿态引起的预测多普勒频移fgeo，获得多普勒频移异常fdca，再减去风浪偏差贡献的多普勒中心频率，完成径向海流sar反演；
66、第二反演模块，用于基于径向海流sar反演结果，利用二维ekman流反演方法，实现二维实时海流矢量反演。
67、本发明的有益效果是：本发明提出一种面向二维海表流场sar反演的多普勒中心频率估计方法，能够提升多普勒中心频率估计结果的成像质量和精度；是对现有sar海流反演中多普勒中心频率估计方法的有效补充，能够提升sar海流反演精度，具有重要的实际应用意义。
1.一种面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，包括：
2.根据权利要求1所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤2包括：
3.根据权利要求1所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤2.2包括：
4.根据权利要求3所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤3包括：
5.根据权利要求4所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤3.2包括：
6.根据权利要求5所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤3.3中，用于搜索卷积运算结果过零点的积分方程为：
7.根据权利要求6所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤4中，信噪比计算公式为：
8.根据权利要求7所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤5中，迭代计算公式为：
9.根据权利要求8所述的面向二维海表流场sar反演的多普勒中心频率估计方法，其特征在于，步骤6包括：
10.一种面向二维海表流场sar反演的多普勒中心频率估计系统，其特征在于，用于执行权利要求1至9任一所述的方法，包括：
