CN106526596A - 合成孔径雷达反演海面风场变分模型的数据处理方法 - Google Patents

Abstract

本发明提供一种合成孔径雷达反演海面风场变分模型的数据处理方法，针对现有技术中存在的EnVar和MCVar求解变分模型过程中运算量大、耗时的缺点问题，通过泰勒展开、线性化代价函数、基于Armijo非精确线搜索技术，实现迭代求取代价函数极小的快速解法，提高运算时效。

Description

合成孔径雷达反演海面风场变分模型的数据处理方法

技术领域

本发明涉及信息技术领域，尤其是指一种合成孔径雷达反演海面风场变分模型的数据处理方法。

背景技术

海面风场资料在气象预报、海浪预报、大气海流系统的长期变化和海洋资源的开发应用等方面起着十分重要的作用；因此对海面风场的探测已成为当前气象学、海洋学研究中的一个热点问题。现有技术中普遍采用合成孔径雷达反演海面风场的方法，其优势主要体现在以下几个方面：

(1)空间分辨率高：现代高科技战争对精细化海面风场提出更高要求，合成孔径雷达的空间分辨率为几十米甚至几米，远高于散射计和辐射计，所以更适合于近岸风速反演；

(2)覆盖范围广：合成孔径雷达可获取全球探测数据，利用反演算法计算得到全球海面风场，所以能够为所有海洋区域渔业提供海面风场数据保障；

(3)数据可靠且多源：随着技术的发展，国内正在开展星载和机载合成孔径雷达的研制工作，即将填补我国高分辨率海洋动力环境卫星的空白，使我国免于受制于国外数据封锁，获取可靠、稳定的数据。欧空局、美国宇航局等组织的合成孔径雷达亦在轨运行，可通过数据共享来进一步提高数据时间、空间分辨率。

合成孔径雷达反演海面风场相关研究始于20世纪70年代末，Weissman等(1979)分析指出合成孔径雷达雷达强度图像与海面风场存在相关性，其中强度图像的风条纹方向与海面风向基本一致，强度大小与风速相关。随后，大量学者开展了卓有成效的研究，主要集中在以下几个研究方向：

1)从合成孔径雷达强度图像的风条纹中提取风向，代表性的方法有傅里叶变换方法(Vachon等，1996)、Sobel算子方法(Koch，2004)、数值微分方 法(姜祝辉，2011)和小波变换方法等(Zecchetto等，2002)；

2)将合成孔径雷达图像中反演得到的海面风向或数值预报、浮标资料中的海面风向当作反演的海面风向，并将其输入到地球物理模型中，直接反演海面风速；

3)借鉴散射计反演海面风场原理，假定相邻风矢量非常接近，利用最大似然法反演海面风场(He等，2005)；

4)基于贝叶斯理论，将后向散射截面与数值预报结果或其他卫星资料反演结果中的海面风场引入到变分方程中，通过求取代价函数极小来确定(Portabella等，2002；Choisnard等，2008)。

Horstmann等(2000)指出，直接用散射计的地球物理模型函数反演合成孔径雷达海面风场时，风向误差将直接导致风速误差。更进一步，由于地球物理模型函数本身的对称性，同样的风向误差导致的风速误差在45°，135°，225°和315°附近最大；低风速情况下，随着入射角增大误差增大；高风速情况下，随着入射角增大误差反而减小；风速精度亦严重依赖于仪器性能。Portabella等(2002)针对上述不确定性，首次提出了利用统计反演方法来反演风矢量的方法，该方法基于贝叶斯理论，将后向散射截面与通过数值预报结果得到的风向先验估计引入到变分方程中，通过求取代价函数极小来确定最优风矢量。Choisnard等(2008)对该方法在不同背景风向情况下的反演结果进行了评估。Dagestad等(2010)将多普勒频移项加入上述代价函数，试验证明在气旋和锋面天气情况下反演结果更符合实际情况。

张毅等在2007年对CMOD4地球物理模型函数进行了仿真，结果表明：后向散射截面在不同程度上依赖于风速、风向、入射角和极化方式。张毅等(2010)、宋贵霆(2007)、杨劲松等(2001)、陈艳玲等(2007)、喻亮等(2005)在合成孔径雷达反演海面风场研究领域开展了创新性研究。在合成孔径雷达反演海面风场变分模型方面，姜祝辉等(2011)给出了考虑正则化参数的变分方法反演海面风场的新思路，姜祝辉等(2014)进一步分析了合成孔径雷达反演海面风场变分模型。

Portabella等(2002)将输入数值预报结果中的海面风场v_B作为背景场(下标B代表背景场)，将合成孔径雷达后向散射截面作为分析场，考虑到 两者均存在误差，构建代价函数，进一步建立了合成孔径雷达反演海面风场变分模型：

其中v为待求海面风场(又称分析场)；σ°为将v和合成孔径雷达其他参数一同输入到CMOD5地球物理模型函数中计算得到的后向散射截面；Δσ°和Δv分别为合成孔径雷达和数值预报结果中的海面风场标准差。

若将v分解为海面风速水平分量u和海面风速垂直分量v，可得代价函数

其中Δσ、Δu和Δv分别为合成孔径雷达后向散射截面标准差、数值预报结果中的海面风场水平分量u和垂直分量v的标准差；其中σ°是(u,v)的函数，即σ°(u,v)，通过求取代价函数J(u,v)的极小值得出的(u,v)就是最终的待求海面风场。

如图1所示的，围绕坐标原点的等值线示意了(2)式右侧第一项，网格中的等值线示意了(2)式右侧第二项和第三项之和，等值线的粗细与最大似然估计的大小成正比。变分模型的最优解是求解J(u,v)的最小值点，图1中间断箭头处就是代价函数极小值点。变分模型的输入输出见图2所示。

其中CMOD5的表达形式为(Hersbach，2007)

σ°＝B₀(1+B₁cosφ+B₂cos2φ)^1.6

其中B₀,B₁,B₂是风速V和入射角θ的函数，φ为风向；

B₀，B₁，B₂定义为

B₂＝(-d₁+d₂v₂)exp(-v₂)

令

其中α＝s₀(1-g(s₀))，

v₂的定义为

以上函数中的a₀,a₁,a₂,γ,s₀,v₀,d₁和d₂为入射角的函数

a₀＝c₁+c₂x+c₃x²+c₄x³

a₁＝c₅+c₆x

a₂＝c₇+c₈x

γ＝c₉+c₁₀x+c₁₁x²

s₀＝c₁₂+c₁₃x

v₀＝c₂₁+c₂₂x+c₂₃x²

d₁＝c₂₄+c₂₅x+c₂₆x²

d₂＝c₂₇+c₂₈x

y₀＝c₁₉

n＝c₂₀。

现有技术中，对合成孔径雷达反演海面风场变分模型进行处理以获得该变分模型的解的方法主要包括枚举法和蒙特卡洛法：

(1)枚举法：Portabella等(2002)给出了枚举法求解变分模型的技术方案(简称EnVar)；其具体步骤包括：

①将待求海面风场在水平方向u上每间隔0.25m/s赋值(假定u的取值范围为-20m/s～20m/s)，即

u＝[-20m/s，-19.75m/s，-19.5m/s，......，-0.25m/s，0m/s，......，19.5m/s，19.75m/s，20m/s]

其维数为1×161；

②将待求海面风场在垂直方向v上每间隔0.25m/s赋值(假定v的取值范围为-20m/s～20m/s)，即

v＝[-20m/s，-19.75m/s，-19.5m/s，......，-0.25m/s，0m/s，......，19.5m/s，19.75m/s，20m/s]

其维数为1×161；

③将以上海面风场水平分量u和垂直分量v组成(u,v)，共形成161×161＝25921组待定海面风场解；图1中的网格就示意了(u,v)的待定海面风场解；

④将25921组待定海面风场解分别代入(2)式中，求得25921个J(u,v)的解；

⑤计算min[J(u,v)]，得出最优(u,v)。

而EnVar的缺点是显而易见的，即需要首先求取多达25921个J(u,v)的解，然后再逐个比较得出最优解，计算速度慢。

(2)蒙特卡洛法

Choisnard等(2008)给出了蒙特卡洛法求解变分模型的技术方案(简称MCVar)；其具体步骤包括：

①假定u、v的取值范围均为-30m/s～30m/s，在背景场上添加高斯误差产生2000个样本，作为待定海面风场；

②将上述2000个待定海面风场解分别代入(2)式中，求得2000个J(u,v)的解；

③计算min[J(u,v)]，得出最优(u,v)。

MCVar的缺点是需要首先求取多达2000个J(u,v)的解，然后再逐个比较得出最优解。其迭代次数虽然较枚举法少，但运算量依旧很大，相当耗时。

由此可见，虽然当前的学者提出了利用变分模型反演合成孔径雷达海面风场的新思路，但变分模型求解速度慢的问题制约了其实用性。

发明内容

针对现有技术中存在的EnVar和MCVar求解变分模型过程中运算量大、耗时的缺点问题，本发明要解决的技术问题是提出一种基于Armijo非精确线搜索的数据处理方法，通过泰勒展开、线性化代价函数、基于Armijo非精确 线搜索技术，实现迭代求取代价函数极小的快速解法，提高运算时效。

为了解决上述问题，本发明实施例提出了一种合成孔径雷达反演海面风场变分模型的数据处理方法，所述合成孔径雷达反演海面风场变分模型为：

其中v为待求海面风场；σ°为将v和合成孔径雷达其他参数一同输入到CMOD5地球物理模型函数中计算得到的后向散射截面；Δσ°和Δv分别为合成孔径雷达和数值预报结果中的海面风场标准差；

将v分解为海面风速水平分量u和海面风速垂直分量v，得到其代价函数

所述处理方法包括：

步骤1、对(1)式在v_B处进行泰勒展开并保留前三项，即

计算(3)式的稳定点

设是非奇异矩阵，得到迭代公式

其中为迭代方向；

步骤2、选取参数β∈(0,1),σ∈(0,0.5)，基于“Armijo condition”，记m_k为满足下列不等式的最小非负整数m：

由此可确定迭代步长在“Armijo condition”下，(5)式可修改为

v＝v_B+α_kd_k (7)；

步骤3、计算代价函数导数，即J(v)的一阶导数和二阶导数包括：

将v展开为(u,v)，将σ°(u,v)在(u_B,v_B)进行泰勒展开，即

将(8)式代入(2)式则

令则

进而

其中(u_G,v_G)可从微分定义出发

这里取δu＝δv＝0.1m/s，则

将(10)式和(11)式代入(9)式，通过迭代即可求出待定海面风场；

步骤4、求解代价函数，包括：

步骤41、选取迭代初始值v_k＝v_B；

步骤42、计算将v_k输入到CMOD5地球物理模型函数中求得进而求得通过(11)式求得v_Gk；计算Δσ、Δu和Δv；将上述量代入(10)式即可求得如果则停算，输出v_a＝v_k；

步骤43、计算将步骤42中求得的v_Gk和中Δσ、Δu和Δv的值代入(11)式即可求得

步骤44、求解线性方程组得出d_k；

步骤45、确定迭代步长α_k；

步骤46、令v_k+1＝v_k+α_kd_k，返回步骤42。

本发明的上述技术方案的有益效果如下：上述技术方案提出了一种合成孔径雷达反演海面风场变分模型的数据处理方法，能够使得合成孔径雷达反演海面风场变分模型求解速度加快，仅通过平均7次左右的迭代即可得到海面风场解，提高了运算效率。

附图说明

图1为现有的变分模型概念图；

图2为图1中的变分模型的输入、输出概念图；

图3为现有的2000次蒙特卡洛方法求解变分模型过程中模拟的u、v分量示意图；

图4为本发明实施例的ArVar算法的流程图；

图5为真实风速为1m/s、背景风速误差为2m/s时，利用ArVar进行求解过程中迭代次数与之间的关系；

图6为真实风速为25m/s、背景风速误差为-1m/s时，利用ArVar进行求解过程中迭代次数与之间的关系；

图7为风速误差随背景场风向误差增大的变化的示意图；

图8为风向反演结果与背景场风向之差随背景场风向误差增大的变化示意图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明实施例提出了一种合成孔径雷达反演海面风场变分模型的数据处理方法，通过泰勒展开、线性化代价函数、基于Armijo非精确线搜索技术，实现迭代求取代价函数极小的快速解法(简称为ArVar)，提高运算时效。

(A)ArVar理论推导

因J(v)的二阶导数(Hesse矩阵)连续，所以前述的(1)式可在v_B处进行泰勒展开并保留前三项，即

求其稳定点，可得

假定是非奇异矩阵，可得迭代公式

其中为迭代方向。

进一步，选取参数β∈(0,1),σ∈(0,0.5)，基于“Armijo condition”，记m_k为满足下列不等式的最小非负整数m：

由此可确定迭代步长在“Armijo condition”下，(5)式可修改为

v＝v_B+α_kd_k (7)。

(B)代价函数导数的计算

现进一步求J(v)的一阶导数和二阶导数

将v展开为(u,v)，将σ°(u,v)在(u_B,v_B)进行泰勒展开，即

将(8)式代入(2)式有

令则

进而

其中(u_G,v_G)可从微分定义出发

这里取δu＝δv＝0.1m/s，则

将(10)式和(11)式代入(9)式，通过迭代即可求出待定海面风场。

(C)代价函数的求解步骤，具体如图4所示的，包括：

①选取迭代初始值v_k＝v_B；

②计算将v_k输入到CMOD5地球物理模型函数中求得进而求得通过(11)式求得v_Gk；其中Δσ、Δu和Δv的值见图2，将上述量代入(10)式即可求得如果则停算，输出v_a＝v_k；

③计算将②中求得的v_Gk和图2中Δσ、Δu和Δv的值代入(11)式即可求得

④求解线性方程组得出d_k；

⑤确定迭代步长α_k；

⑥令v_k+1＝v_k+α_kd_k，返回步骤②。

以下通过一个具体的实例对本发明实施例中提出的合成孔径雷达反演海面风场变分模型的数据处理方法进行模拟实验分析：

令合成孔径雷达后向散射截面无误差，所以观测得到的合成孔径雷达后向散射截面可用真实风速风向代入CMOD5地球物理模型函数中计算得出，合成孔径雷达极化方式为VV极化，入射角θ＝25°，雷达天线方位角为0°，令(6)式中的β＝0.55，σ＝0.4，4.2.3节中ε＝0.01，最大迭代次数为50。令背景场风速＝真实风速+背景场风速误差，背景场风向＝真实风向+背景场风向误差。

迭代次数对比分析

假定真实风向为Fi＝0°，真实风速为V＝[1m/s，3m/s，5m/s，......，23m/s，25m/s]，其维数为1×13，背景场风向误差为eFi＝20°，背景场风速误差为eV＝[-2.5m/s，-2m/s，-1.5m/s，......，2m/s，2.5m/s]T，其维数为11×1。

将11×13组背景场组合转化为(u_B,v_B)，并与一同代入4.2.3节进行计算，得到迭代次数见表1所列(由于本模拟试验定义的误差规则导致前四组数据的背景风速为负值，所以前四组数据不参与迭代，直接将背景场赋值给分析场，表1中用“—”填充)。

表111×13组模拟试验迭代次数

从表1中可见当真实风速为1m/s，背景风速误差分别为0.5m/s、1m/s、1.5m/s、2m/s和2.5m/s时迭代次数为最大迭代次数50次，其余各组试验迭代次数均在13次以下，最小迭代次数仅为1次，并呈现出表1中左下区域迭代次数多，右上区域次之，其余区域迭代次数相对最少的趋势，其原因可能与两个因素相关：①背景风速误差趋近于0m/s时候，背景风速与真实风速接近，迭代次数自然较少；②背景风向为20°不变，当真实风速变大时，背景风速误差的增大可能会促进迭代次数的减少。

下面在真实风速为1m/s，背景风速误差为2m/s和真实风速为25m/s，背景风速误差为-1m/s两种情况下，分析利用ArVar进行求解过程中迭代次数与之间的关系，见图5和图6。

图5代表了迭代次数达最大迭代次数的情况，从图中可见，在第三次迭代之后已经达到了极小值，然而该极小值大于本发明设置的阈值ε＝0.01，所以迭代过程中只能在该极小值附近反复震荡，直至最大迭代次数。

图6代表了一次普通的迭代过程，前三次迭代的值下降速度非常快，这种现象在图5中也有体现，随后缓慢下降，直至达到ε＝0.01后，跳出迭代循环。

由此可见，ArVar提出的技术方案能够实现平均7次迭代即可求出海面风场解。

下面从电脑耗时角度出发，ArVar分别与EnVar和MCVar对比分析。试验选用的电脑配置为：型号：ThinkCenter M8300t，处理器：Intel Core i7-26003.40GHz，内存3GB。

从表2中可见，计算全部上述11×13组试验数据，EnVar和MCVar单独迭代一次耗时0.014s，而两个方法的迭代次数分别为25921次和2000次，不考虑计算的耗时，EnVar耗时362.9s，MCVar耗时28.0s，而本发明平均迭代次数仅为7次，总耗时仅为1.3s。可见本发明大大降低了计算耗时。

表2 ArVar分别与枚举法和蒙特卡洛法对比分析

算法精度对比分析

以下分别利用ArVar和EnVar进行反演计算。

(a)背景场风速无误差，背景场风向有误差情况

假定真实风向为Fi＝45°，真实风速为V＝8m/s，背景场风向误差eFi＝[10°，20°，30°，......，170°，180°]T，背景场风速误差eV＝0m/s。利用ArVar和EnVar进行反演计算得出的结果见图7和图8。

从图7中可见，风速误差随着背景场风向误差增大而呈周期性变化。背景场风向为0°、90°和180°时，风速误差最小，背景场风向为45°和135°时风速误差最大。整体上说，EnVar与ArVar误差相近，EnVar误差大于ArVar误差。

ArVar和EnVar风向反演结果误差与背景场风向之间的关系见图8。ArVar和EnVar风向反演结果误差随背景场风向增大而增大，ArVar风向反演结果误差为55°，EnVar风向反演结果误差为56°，ArVar风向反演结果误差略优于EnVar。

(b)背景场风速有误差，背景场风向无误差情况：

假定真实风向为Fi＝45°，真实风速为V＝8m/s，背景场风向误差eFi＝0°，背景场风速误差eV＝[-5m/s，-4m/s，-3m/s，......，4m/s，5m/s]。利用两种方法进行反演计算得出的结果见图9和图10。

从图9可见，ArVar和EnVar风速反演结果误差随背景场风速的增大而增大，ArVar风速误差略大于EnVar风速误差。

从图10可见，ArVar和EnVar风向反演结果误差随背景场风速的增大而增大，ArVar风向误差略小于EnVar风向误差。

本发明实施例的参考文献如下，本发明实施例中将这些参考文献全文引用于此：

[1]陈艳玲,黄玻,丁晓利等.ERS-2SAR反演海洋风矢量的研究.地球物理学报,2007,50(6):1688-1694.

[2]姜祝辉，黄思训，何然，周晨腾.合成孔径雷达资料反演海面风场的正则化方法研究.物理学报,2011,60(6):068401.

[3]姜祝辉,黄思训,石汉青等.合成孔径雷达图像反演海面风向新方法的研究.物理学报,2011,60(10):108402.

[4]姜祝辉，周晓中，游小宝，易欣，黄为权.合成孔径雷达反演海面风场变分模型分析.物理学报,2014,63(14):148401.

[5]宋贵霆.合成孔径雷达提取海面风、浪信息研究.北京:中国科学院研究生院博士学位论文,2007.

[6]杨劲松,黄韦艮,周长宝.合成孔径雷达图像的近岸海面风场反演.遥感学报,2001,5(1):13-16.

[7]喻亮,丁晓松.利用星载ERS-2SAR进行长江口海面风场反演研究.信息与电子工程,2005,3(3):172-175.

[8]张毅,陈永强,朱敏慧.SAR海面风场反演研究.电子测量技术,2007,30(2):36-38.

[9]张毅,蒋兴伟,林明森等.合成孔径雷达海面风场反演研究进展.遥感技术与应用,2010,25(3):423-429.

[10]Choisnard J,Laroche S.Properties of variational data assimilationfor synthetic aperture radar wind retrieval.J Geophys Res,2008,113,C05006,doi:10.1029/2007JC004534.

[11]Dagestad K F,Mouche A,Collard F,et al.On the use of doppler shiftfor SAR wind retrieval.Proceedings of SeaSAR Workshop,ESRIN Jan,2010.

[12]He Y J,Perrie W,Zou Q P,et al.A new wind vector algorithm for C-band SAR.IEEE Transactions on Geoscience and Remote Sensing,2005,43(7):1453-1458.

[13]Hersbach H,Stoffelen A,de Haan S.An improved C-band scatterometerocean geophysical model function:CMOD5.J Geophys Res,2007,112,C03006,doi:10.1029/2006JC003743.

[14]Horstmann J,Koch W,Lehner S.Wind field retrieval using satellitebased synthetic aperture radars.IEEE International Geoscience and RemoteSensing Symposium,2000,4:1501-1503.

[15]Koch W.Directional analysis of SAR images aiming at winddirection.IEEE Transactions on Geoscience and Remote Sensing,2004,42(4):702-710.

[16]Portabella M,Stoffelen A.Toward an optimal inversion method forsynthetic aperture radar wind retrieval.J Geophys Res,2002,107(C8),3086,doi:10.1029/2001JC000925.

[17]Vachon P W,Dobson F W.Validation of wind vector retrieval fromERS-1SAR images over the ocean.Global Atmos Ocean Syst,1996,5:177-187.

[18]Weissman D E,King D,Thompson T W.Relationship between hurricanessurface winds and L-band radar backscatter from the sea surface.J AppliedMeteorol,1979,18:1023-1034.

[19]Zecchetto S,De Biasio F.Wind field retrieval from SAR imagesusing the continuous wavelet transform.IEEE Transactions on Geoscience andRemote Sensing Symposium,2002,4:1974-1976.

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种合成孔径雷达反演海面风场变分模型的处理方法，其特征在于，所述合成孔径雷达反演海面风场变分模型为：

$J\left(v\right)={\left(\frac{{\mathrm{\σ}}_m^o-{\mathrm{\σ}}^o}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+{\left(\frac{v_B-v}{\mathrm{\Δ}v}\right)}^2---\left(1\right)$

其中v为待求海面风场；σ^o为将v和合成孔径雷达其他参数一同输入到CMOD5地球物理模型函数中计算得到的后向散射截面；Δσ^o和Δv分别为合成孔径雷达和数值预报结果中的海面风场标准差；

将v分解为海面风速水平分量u和海面风速垂直分量v，得到其代价函数

$J(u,v)={\left(\frac{{\mathrm{\σ}}_m^o-{\mathrm{\σ}}^o}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+{\left(\frac{u_B-u}{\mathrm{\Δ}u}\right)}^2+{\left(\frac{v_B-v}{\mathrm{\Δ}v}\right)}^2---\left(2\right)$

所述处理方法包括：

步骤1、对(1)式在v_B处进行泰勒展开并保留前三项，即

$J\left(v\right)=J\left(v_B\right)+{(\▿J(v\left){|}_B\right)}^T(v-v_B)+\frac{1}{2}{(v-v_B)}^T{\▿}^{2}J\left(v\right){|}_B(v-v_B)---\left(3\right)$

计算(3)式的稳定点

$\▿J\left(v\right)=\▿J\left(v\right){|}_B+{\▿}^{2}J\left(v\right){|}_B(v-v_B)=0---\left(4\right)$

设是非奇异矩阵，得到迭代公式

$v=v_B-{\▿}^{2}J\left(v\right){|}_B^{-1}\▿J\left(v\right)---\left(5\right)$

其中为迭代方向；

步骤2、选取参数β∈(0,1),σ∈(0,0.5)，基于Armijo非精确线搜索技术Armijocondition，记m_k为满足下列不等式的最小非负整数m：

$J(v_k+{\mathrm{\β}}^m{d}_k)\≤J\left(v_k\right)+{\mathrm{\σ\β}}^m{(\▿J(v\left){|}_k\right)}^T{d}_k---\left(6\right)$

由此可确定迭代步长在基于Armijo非精确线搜索技术Armijo condition下，(5)式可修改为

v＝v_B+α_kd_k (7)；

步骤3、计算代价函数导数，即J(v)的一阶导数和二阶导数包括：

将v展开为(u,v)，将σ^o(u,v)在(u_B,v_B)进行泰勒展开，即

${\mathrm{\σ}}^o(u,v)\≈{\mathrm{\σ}}_B^o(u_B,v_B)+\frac{\∂{\mathrm{\σ}}^o(u,v)}{\∂u}{|}_B(u-u_B)\frac{\∂{\mathrm{\σ}}^o(u,v)}{\∂v}{|}_B(v-v_B)---\left(8\right)$

将(8)式代入(2)式则

$J(u,v)={\left(\frac{{\mathrm{\σ}}_m^o-\left({\mathrm{\σ}}_B^o\right(u_B,v_B)+\frac{\∂{\mathrm{\σ}}^o(u,v)}{\∂u}{|}_B(u-u_B)+\frac{\∂{\mathrm{\σ}}^o(u,v)}{\∂v}{|}_B(v-v_B))}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+{\left(\frac{u_B-u}{\mathrm{\Δ}u}\right)}^2+{\left(\frac{v_B-v}{\mathrm{\Δ}v}\right)}^2---\left(9\right)$

令则

$J(u,v)={\left(\frac{\mathrm{\σ}*-u_G(u-u_B)-v_G(v-v_B)}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+{\left(\frac{u-u_B}{\mathrm{\Δ}u}\right)}^2+{\left(\frac{v-v_B}{\mathrm{\Δ}v}\right)}^2$

进而

$\▿J(u,v)=\left(\begin{array}{c}-u_G\frac{\mathrm{\σ}*-u_G(u-u_B)-v_G(v-v_B)}{{\left(\mathrm{\Δ}\mathrm{\σ}\right)}^2}+\frac{u-u_B}{{\left(\mathrm{\Δ}u\right)}^2}\\ -v_G\frac{\mathrm{\σ}*-v_G(u-u_B)-v_G(v-v_B)}{{\left(\mathrm{\Δ}\mathrm{\σ}\right)}^2}+\frac{v-v_B}{{\left(\mathrm{\Δ}v\right)}^2}\end{array}\right)---\left(10\right)$

${\▿}^{2}J(u,v)=\left(\begin{array}{cc}{\left(\frac{u_G}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+\frac{1}{{\left(\mathrm{\Δ}u\right)}^2}& \frac{u_G{v}_G}{{\left(\mathrm{\Δ}\mathrm{\σ}\right)}^2}\\ \frac{u_G{v}_G}{{\left(\mathrm{\Δ}\mathrm{\σ}\right)}^2}& {\left(\frac{v_G}{\mathrm{\Δ}\mathrm{\σ}}\right)}^2+\frac{1}{{\left(\mathrm{\Δ}v\right)}^2}\end{array}\right)---\left(11\right)$

其中(u_G,v_G)可从微分定义出发

$\left\{\begin{array}{c}{u}_G=\frac{\∂{\mathrm{\σ}}^o(u,v)}{\mathrm{\δ}u}{|}_B=\underset{\mathrm{\δ}u\→0}{\lim }\frac{{\mathrm{\σ}}^o(u_B+\mathrm{\δ}u,v_B)-{\mathrm{\σ}}^o(u_B,v_B)}{\mathrm{\δ}u}\\ u_G=\frac{\∂{\mathrm{\σ}}^o(u,v)}{\mathrm{\δ}v}{|}_B=\underset{\mathrm{\δ}v\→0}{\lim }\frac{{\mathrm{\σ}}^o(u_B,v_B+\mathrm{\δ}v)-{\mathrm{\σ}}^o(u_B,v_B)}{\mathrm{\δ}V}\end{array}\right.---\left(12\right)$

这里取δu＝δv＝0.1m/s，则

$\left\{\begin{array}{c}{u}_G=\frac{{\mathrm{\σ}}^o(u_B+0.1,v_B)-{\mathrm{\σ}}^o(u_B,v_B)}{0.1}\\ v_G=\frac{{\mathrm{\σ}}^o(u_B,v_B+0.1)-{\mathrm{\σ}}^o(u_B,v_B)}{0.1}\end{array}\right.---\left(13\right)$

将(10)式和(11)式代入(9)式，通过迭代即可求出待定海面风场；

步骤4、求解代价函数。

2.根据权利要求1所述的一种合成孔径雷达反演海面风场变分模型的处理方法，其特征在于，所述步骤4具体包括：

步骤41、选取迭代初始值v_k＝v_B；

步骤42、计算将v_k输入到CMOD5地球物理模型函数中求得进而求得通过(11)式求得v_Gk；计算Δσ、Δu和Δv；将上述量代入(10)式即可求得如果则停算，输出v_a＝v_k；

步骤43、计算将步骤42中求得的v_Gk和中Δσ、Δu和Δv的值代入(11)式即可求得

步骤44、求解线性方程组得出d_k；

步骤45、确定迭代步长α_k；

步骤46、令v_k+1＝v_k+α_kd_k，返回步骤42。