CN105701284A - 时变多尺度电大区域海面电磁散射矢量场的并行计算方法 - Google Patents

Abstract

本发明公开了一种时变多尺度电大区域海面电磁散射矢量场的并行计算方法，包括步骤：建立重力波海面模型；重力波与张力波的截断频率k_c的确定，利用基尔霍夫求解重力波下的相干散射，利用小扰动近似求解张力波的Bragg散射；利用重力波调制张力波谱得到各面元下张力波的矢量场；利用Z-BUFFER方法进行面元遮蔽判断；利用场的矢量叠加原理确定各面元下的矢量场的幅度与相位。本发明基于重力波调制张力波谱以及电大区域粗糙海面的时变多尺度回波特性特点，将张力波利用重力波进行调制，利用基尔霍夫和小扰动近似法结合CUDA并行计算大区域粗糙海面的散射场、和多普勒谱，实现时变多尺度电大区域粗糙海面雷达回波特性分析。

Description

时变多尺度电大区域海面电磁散射矢量场的并行计算方法

技术领域

本发明涉及一种电大区域时变多尺度海面电磁散射矢量场时间序列计算方法，采用基于重力波调制张力波谱的时变多尺度电大区域海面电磁散射矢量场时间序列的并行计算方法，重力波调制张力波谱的高频近似方法可以求解电大区域海面矢量场时间序列，并行计算可以提高计算效率。

背景技术

电磁散射研究在海洋上的遥感、目标探测与追踪以及SAR成像等领域起着极为重要的作用，一直以来传统的高频计算海面电磁散射方法是通过求集平均的方法来统计得到海面总回波散射系数的，在实际的电大尺寸海面电磁散射雷达探测中，通常接受到的信号为一系列时间序列，时间序列中包含了海况等各种海态参数信息，是雷达进行目标探测的主要信息来源，然而通过实验获取这些杂波时间序列成本高昂，所以通过计算机仿真，对动态海面进行电磁散射建模不失为一种有效而快速的研究方法。

针对不同的电磁波参数，电大区域时变多尺度海面对应了不同的散射机制，其中大尺度重力波对应粗糙海面下的相干散射，小尺度张力波对应重力波调制下的布拉格散射，而对于重力波与张力波的截断频率一直是一个难点，早期简单基于电磁波波数，缺乏针对性变化。随着海面区域的增大，海面的剖分面元也随着增大，对时变电大区域海面的散射场矢量时间序列的建模造成制约，根据目前检索国内外资料表明，分析电大区域多尺度粗糙海面散射矢量场的时间序列特性还未见报道。

发明内容

本发明所要解决的技术问题是：提供一种基于重力波调制张力波谱的时变多尺度电大区域海面电磁散射矢量场时间序列的并行计算方法，基于叠加法的大尺度重力波海面并行建模方法，考虑重力波的相干散射，及重力波的倾斜调制和毛细波的相位调制，利用截断频率将粗糙海面分为大尺度重力波与小尺度张力波，，利用Z-Buffer方法对面元进行遮蔽判断，利用矢量叠加原理计算动态海面的杂波时间序列和多普勒谱，并利用CUDA技术进行并行计算。

为解决上述技术问题，本发明的技术方案是：一种基于重力波调制张力波谱的时变多尺度电大区域海面电磁散射矢量场时间序列的并行计算方法，包括以下步骤：

1)基于叠加法的大尺度重力波海面并行建模；

2)基于基尔霍夫计算条件与小扰动计算条件，给出不同频率不同海情下大尺度重力波与小尺度张力波的截断波数k_c

3)基于截断波数k_c，当面元为大尺度重力波时，利用面元基尔霍夫方法求得相干散射场，当面元为小尺度张力波时，利用面元小扰动近似方法得到基于大尺度面元本地坐标系下的张力波Bragg散射场，并利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断。

4)利用重力波调制张力波谱的面元法及场的矢量叠加原理求得电大区域时变多尺度海面的电磁场的时间序列，最终得到场的多普勒谱。。

5)基于GPU利用CUDA对每个海面的面元进行并行，每个kernel计算单个面元的散射场，同时进行计算，并对时间循环，利用异步传输理论，提高并行效率，大大提高计算效率。

所述的步骤1)中基于叠加法的大尺度重力波海面并行建模公式如下：

${\mathrm{\ζ}}_g(x,y,t)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{2S({\mathrm{\ω}}_i,{\mathrm{\θ}}_j-{\mathrm{\θ}}_w)\mathrm{\Δ}k\mathrm{\Δ}\mathrm{\θ}}\cos ({\mathrm{\ω}}_{i}t-k_{i}x{\mathrm{cos\θ}}_j-k_{i}y{\mathrm{sin\θ}}_j+{\mathrm{\ϵ}}_{ij})$

$S(k,\mathrm{\φ})=\frac{1}{k}{S}_{1d}\left(k\right)S_{dir}(k,\mathrm{\φ})$

S_1d(k)＝(B_L+B_H)/k³,S_dir＝[1+Δ(k)cos(2φ)]/2π

其中，S(ω_i,θ_j-θ_w)为Efouhialy海谱。k_i、ω_i、θ_j和ε_ij分别表示组成波的波数、圆频率、方位角和初始相位。初始相位ε_ij在0～2π之间满足均匀分布。M，N分别代表频率和方位角的采样点数。相邻两点间的距离分别是Δx和Δy。充分利用二维海面仿真的内在并行性，采用粗粒度线程的并行策略，即一个GPU线程完成一个采样点的计算，这就意味着我们可以通过更加精细的剖分海面获得更加逼真的海面，但计算量并不会显著增加。

所述的步骤2)基于大尺度重力波分量成立条件与小尺度张力波分量成立条件，给出不同频率不同海情大尺度重力波与小尺度张力波的截断波数k_c。

计算风驱海面的电磁散射中，KA模型适合计算大尺度的重力波，而SPM模型适合计算小尺度的毛细波，这两种模型的计算使用范围和使用条件都不相同，而联系这两种尺度的参量即为截断波数k_cut，即利用k_cut来区分大尺度与小尺度，当k<k_cut时，此时大尺度的重力波为主要影响，应用KA模型计算。当k>k_cut时，这时小尺度的张力波为主要影响，因此用基于重力波调制张力波谱近似方法计算。其中复合表面模型中针对大尺度波分量模型成立条件为

$\frac{2\mathrm{\&pi;}}{k_i}{\mathrm{\&sigma;}}_c<<1,k_i{\mathrm{\&delta;}}_l>\frac{\sqrt{10}}{2\cos \left({\mathrm{\&theta;}}_i\right)}$

式中，k_i为入射电磁波数，θ_i为入射角度，σ_c是平均曲率的标准离差，δ_l为大尺度波面的均方根高度。且有

${\mathrm{\&sigma;}}_c^2={\&Integral;}_0^{k_{cut}}{k}^{4}S\left(k\right)dk,{\mathrm{\&sigma;}}_l^2={\&Integral;}_0^{k_{cut}}S\left(k\right)dk$

复合表面模型中，小尺度波分量成立条件如下：

k_iδ_scos(θ_i)＜＜1

式中，δ_s为小尺度波面的均方根高度

${\mathrm{\&delta;}}_s^2={\&Integral;}_{k_{cut}}^{\mathrm{\&infin;}}S\left(k\right)dk$

优选的，步骤3a)基于截断波数k_c，当k<k_cut时，此时大尺度的重力波为主要影响，应用KA方法进行计算。其公式为：

$\begin{array}{c}{E}_{pq}^{facet}({\hat{k}}_s,{\hat{k}}_i)=ik\frac{e^{ikR}}{4\mathrm{\&pi;}R}{S}_{pq}^{KA}{e}^{-ik(q_x{x}_0+q_y{y}_0+q_z{z}_0)}\\ \&CenterDot;\sin c\&lsqb;k\frac{\mathrm{\&Delta;}x}{2}(q_x+q_z{Z}_x)\&rsqb;\sin c\&lsqb;k\frac{\mathrm{\&Delta;}y}{2}(q_y+q_z{Z}_y)\&rsqb;\mathrm{\&Delta;}x\mathrm{\&Delta;}y\end{array}$

q_x＝k_sx-k_ix,q_y＝k_sy-k_iy,q_z＝k_sz-k_iz

其中，为极化系数，Z_x和Z_y为面元沿x方向与y方向上的斜率。(x₀,y₀,z₀)为每个面元的中心，为入射极化波数矢量，为散射极化波数矢量，p,q＝h,v为入射波与散射波的极化方式。

优选的，步骤3b)中，当k≥k_cut时，这时小尺度的张力波为主要影响，面元为小尺度张力波时，利用面元小扰动近似方法得到基于大尺度面元本地坐标系下的张力波Bragg散射场，其公式为：

$\begin{array}{c}{S}_{pq}^{SPM}({\hat{k}}_s,{\hat{k}}_i)=\frac{k^2(1-\mathrm{\&epsiv;})\mathrm{\&Delta;}S}{8{\mathrm{\&pi;n}}_z}{e}^{-i\hat{q}\&CenterDot;r}{\tilde{F}}_{pq}^{SPM}\&lsqb;B\left(k_c^{+}\right)\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-i)}^n{J}_n\&lsqb;q_{z}B\left(k_c^{+}\right)\&rsqb;I_0\left(k_c^{+}\right)\\ +B\left(k_c^{-}\right)\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-i)}^n{J}_n\&lsqb;q_{z}B\left(k_c^{-}\right)\&rsqb;I_0\left(k_c^{-}\right)\&rsqb;\end{array}$

$\begin{array}{c}{I}_0\left(k_c\right)=e^{-i(1+n){\mathrm{\&omega;}}_{c}t}\sin c\{\frac{{\mathrm{\&Delta;x}}_g}{2}\&lsqb;(1+n)k_{cx}-q_x-q_z{Z}_x\&rsqb;\}\&CenterDot;\sin c\{\frac{{\mathrm{\&Delta;y}}_g}{2}\&lsqb;(1+n)k_{cy}-q_y-q_z{Z}_y\&rsqb;\}\\ +e^{-i(1-n){\mathrm{\&omega;}}_{c}t}\sin c\{\frac{{\mathrm{\&Delta;x}}_g}{2}\&lsqb;(1-n)k_{cx}+q_x+q_z{Z}_x\&rsqb;\}\&CenterDot;\sin c\{\frac{{\mathrm{\&Delta;y}}_g}{2}\&lsqb;(1-n)k_{cy}+q_y+q_z{Z}_y\&rsqb;\}\end{array}$

其中，ε为粗糙海面相对介电常数，为海面单个面元的法线矢量，为本地坐标系下的极化系数，为毛细波幅度，k_c为Bragg波长对应的波数，ΔS为面元大小，J_n(x)为n阶贝塞尔函数。

优选的，步骤3c)利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断。

3c1)、将海面三角面元离散，并将不同的面元渲染不同的颜色，因为每个颜色分量取值为0-255，因此最多可以渲染的面元数为255*255*255≈1600W面元。

3c2)、将海面投影在于视线方向垂直的平面上，因为在OpenGL的绘图过程中，即使是后绘制的物体，但如果后绘制的物体在先绘制物体的后面，那么它将不会被绘制。在实际操作中，调用OpenGL函数，改变视线方向。

3c3)、读取屏幕中的颜色值，因为每一个面元都顺序分配了唯一的RGB颜色值，如果在读取的颜色值中有该面元的颜色值，那么该面元不被遮挡，反之则被其他面元遮挡。

优选的，步骤4)利用重力波调制张力波谱的面元法及场的矢量叠加原理求得电大区域时变多尺度海面的电磁场的时间序列，最终得到场的多普勒谱。

假设二维海面模拟样本在x和y方向上的长度分别是L_x和L_y，面积为A＝L_xL_y，等间隔离散点数为M和N，相邻两点间的距离分别为Δx和Δy，忽略各小面元之间的多次散射作用，可以写出海面总散射矢量场为

$E_{pq}^{total}({\hat{k}}_s,{\hat{k}}_i,t)=\underset{M}{\&Sigma;}\underset{N}{\&Sigma;}\&lsqb;E_{pq}^{SPM}({\hat{k}}_s,{\hat{k}}_i,t)+E_{pq}^{KA}({\hat{k}}_s,{\hat{k}}_i,t)\&rsqb;$

根据上式可以得到一组对应不同时刻海面样本的散射矢量场时间序列，这样时变海面的多普勒谱即可通过标准谱估计方法得到，公式如下

$S_0\left(f\right)=\frac{1}{T}{\left|{\&Integral;}_0^T{E}_{PQ}^{scatt}(\widehat{k_i},\widehat{k_s},t)\exp \left(i2\mathrm{\&pi;}ft\right)dt\right|}^2$

上式可以借助快速傅里叶变换实现。这样，我们容易得到该多普勒的中心频率f_c和展宽γ，公式如下

$f_c=\frac{\&Integral;{\mathrm{fS}}_0\left(f\right)df}{\&Integral;S_0\left(f\right)df}$

${\mathrm{\&gamma;}}^2=\frac{\&Integral;{(f-f_c)}^2{S}_0\left(f\right)df}{\&Integral;S_0\left(f\right)df}$

优选的，步骤5)利用CUDA并行计算将粗糙海面不同面元进行并行，提高计算效率。

为了提高时变多尺度海面后向电磁散射时间序列建模的运算速度，我们通过在CPU-GPU异构平台上设计并实现了一个并行方案。在并行方案中，每一个海面中划分的海面面元为一个thread，最终对时间进行循环，对每个面元进行同时并行，得到单个面元对应的散射场。最后返还到主机端，利用场的矢量叠加原理，得到每个时间点的单个角度的散射场。

采用了上述技术方案，本发明的有益效果为：本发明采用叠加法对时变多尺度海面进行并行建模，利用大尺度重力波分量成立条件与小尺度张力波分量成立条件得到多尺度海面的截断频率，，利用重力波调制张力波谱的高频近似方法得到每个面元的散射场，结合散射场的矢量叠加原理求得总的海面的散射场，最终得到多普勒谱。利用GPU并行原理对时变多尺度海面进行并行计算。本发明所述的方法既有清晰的物理意义，并利用高频方法与并行思想解决了大区域海面电磁散射计算重大技术难题，对分析大区域海面雷达回波特性具有重要的意义。

附图说明

图1是本发明电大区域多尺度时变粗糙海面回波并行计算流程图；

图2是多尺度海面建模

图3通过Z-Buffer判断遮蔽结果

图4是大尺度重力波电磁散射示意图；

图5是微粗糙面电磁散射示意图；

图6是全局直角坐标系；

图7是本地坐标系。

具体实施方式

下面结合实施例对本发明进一步说明。

参照图1，本发明的具体实现步骤如下：

步骤1，电大区域时变多尺度海面并行建模。

利用叠加法基于海谱建立时变多尺度电大区域海面模型，根据不同的雷达参数，根据入射波频率，剖分为不同的面元，结合CUDA高性能并行计算方法，对海面每个点的高度进行并行建模，每个kernel对应每个点，提高建模效率。

叠加法模拟粗糙海面公式如下：

${\mathrm{\&zeta;}}_g(x,y,t)=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\sqrt{2S({\mathrm{\&omega;}}_i,{\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_w)\mathrm{\&Delta;}k\mathrm{\&Delta;}\mathrm{\&theta;}}\cos ({\mathrm{\&omega;}}_{i}t-k_{i}x{\mathrm{cos\&theta;}}_j-k_{i}y{\mathrm{sin\&theta;}}_j+{\mathrm{\&epsiv;}}_{ij})$

$S(k,\mathrm{\&phi;})=\frac{1}{k}{S}_{1d}\left(k\right)S_{dir}(k,\mathrm{\&phi;})$

S_1d(k)＝(B_L+B_H)/k³,S_dir＝[1+Δ(k)cos(2φ)]/2π

其中，S(ω_i,θ_j-θ_w)为Efouhialy海谱。k_i、ω_i、θ_j和ε_ij分别表示组成波的波数、圆频率、方位角和初始相位。初始相位ε_ij在0～2π之间满足均匀分布。M，N分别代表频率和方位角的采样点数。相邻两点间的距离分别是Δx和Δy。根据线性叠加法，假设某时刻t，海上一个固定点的水面波动可以用多个随机余弦波叠加来描述，并假定只在平面内产生波浪，且波浪沿固定方向传播，我们模拟出不同风速与风向条件下的海面模型，但是模拟离散点为512×512的海面模型仍需要耗时532.0s。当涉及到更大范围的海面模拟或剖分网格要求较细密的情况时，离散点数会更高，使得这种简单的嵌套叠加方法生成海面效率更低，将不能满足实时模拟的需求，这种情况下就需要GPU高性能并行计算来提高模拟效率。在我们的并行方案中，充分利用二维海面仿真的内在并行性，采用粗粒度线程的并行策略，即一个GPU线程完成一个采样点的计算，这就意味着我们可以通过更加精细的剖分海面获得更加逼真的海面，但计算量并不会显著增加。

步骤2，利用大尺度重力波分量与小尺度张力波分量成立条件，求得不同风速不同雷达波频率下对应的截断频率k_c。

计算风驱海面的电磁散射中，KA模型适合计算大尺度的重力波，而SPM模型适合计算小尺度的毛细波，这两种模型的计算使用范围和使用条件都不相同，而联系这两种尺度的参量即为截断波数k_cut，即利用k_cut来区分大尺度与小尺度，当k<k_cut时，此时大尺度的重力波为主要影响，应用KA模型计算。当k>k_cut时，这时小尺度的张力波为主要影响，因此用基于重力波调制张力波谱近似方法计算。其中复合表面模型中针对大尺度波分量模型成立条件为

$\frac{2\mathrm{\&pi;}}{k_i}{\mathrm{\&sigma;}}_c<<1,k_i{\mathrm{\&delta;}}_l>\frac{\sqrt{10}}{2\cos \left({\mathrm{\&theta;}}_i\right)}$

式中，k_i为入射电磁波数，θ_i为入射角度，σ_c是平均曲率的标准离差，δ_l为大尺度波面的均方根高度。且有

${\mathrm{\&sigma;}}_c^2={\&Integral;}_0^{k_{cut}}{k}^{4}S\left(k\right)dk,{\mathrm{\&sigma;}}_l^2={\&Integral;}_0^{k_{cut}}S\left(k\right)dk$

复合表面模型中，小尺度波分量成立条件如下：

k_iδ_scos(θ_i)＜＜1

式中，δ_s为小尺度波面的均方根高度

${\mathrm{\&delta;}}_s^2={\&Integral;}_{k_{cut}}^{\mathrm{\&infin;}}S\left(k\right)dk$

步骤3，利用重力波调制张力波谱的高频近似方法求解电大区域时变多尺度海面矢量场时间序列，并利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断。

当k<k_cut时，此时大尺度的重力波为主要影响，应用KA方法进行计算，如图2所示，其公式如下：

$\begin{array}{c}{E}_{pq}^{facet}({\hat{k}}_s,{\hat{k}}_i)=ik\frac{e^{ikR}}{4\mathrm{\&pi;}R}{S}_{pq}^{KA}{e}^{-ik(q_x{x}_0+q_y{y}_0+q_z{z}_0)}\\ \&CenterDot;\sin c\&lsqb;k\frac{\mathrm{\&Delta;}x}{2}(q_x+q_z{Z}_x)\&rsqb;\sin c\&lsqb;k\frac{\mathrm{\&Delta;}y}{2}(q_y+q_z{Z}_y)\&rsqb;\mathrm{\&Delta;}x\mathrm{\&Delta;}y\end{array}$

其中，为极化系数，Z_x和Z_y为面元沿x方向与y方向上的斜率。(x₀,y₀,z₀)为每个面元的中心，为入射极化波数矢量，为散射极化波数矢量，p,q＝h,v为入射波与散射波的极化方式。

当k≥k_cut时，这时小尺度的张力波为主要影响，面元为小尺度张力波时，利用面元小扰动近似方法得到基于大尺度面元本地坐标系下的张力波Bragg散射场，如图3与图4所示，其公式如下：

$\begin{array}{c}{S}_{pq}^{SPM}({\hat{k}}_s,{\hat{k}}_i)=\frac{k^2(1-\mathrm{\&epsiv;})\mathrm{\&Delta;}S}{8{\mathrm{\&pi;n}}_z}{e}^{-i\hat{q}\&CenterDot;r}{\tilde{F}}_{pq}^{SPM}\&lsqb;B\left(k_c^{+}\right)\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-i)}^n{J}_n\&lsqb;q_{z}B\left(k_c^{+}\right)\&rsqb;I_0\left(k_c^{+}\right)\\ +B\left(k_c^{-}\right)\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-i)}^n{J}_n\&lsqb;q_{z}B\left(k_c^{-}\right)\&rsqb;I_0\left(k_c^{-}\right)\&rsqb;\end{array}$

$\begin{array}{c}{I}_0\left(k_c\right)=e^{-i(1+n){\mathrm{\&omega;}}_{c}t}\sin c\{\frac{{\mathrm{\&Delta;x}}_g}{2}\&lsqb;(1+n)k_{cx}-q_x-q_z{Z}_x\&rsqb;\}\&CenterDot;\sin c\{\frac{{\mathrm{\&Delta;y}}_g}{2}\&lsqb;(1+n)k_{cy}-q_y-q_z{Z}_y)\}\\ +e^{-i(1-n){\mathrm{\&omega;}}_{c}t}\sin c\{\frac{{\mathrm{\&Delta;x}}_g}{2}\&lsqb;(1-n)k_{cx}+q_x+q_z{Z}_x\&rsqb;\}\&CenterDot;\sin c\{\frac{{\mathrm{\&Delta;y}}_g}{2}\&lsqb;(1-n)k_{cy}+q_y+q_z{Z}_y)\}\end{array}$

其中，ε为粗糙海面相对介电常数，为海面单个面元的法线矢量，为本地坐标系下的极化系数，为毛细波幅度，k_c为Bragg波长对应的波数，ΔS为面元大小，J_n(x)为n阶贝塞尔函数。

利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断，得到电大区域时变海面电磁散射矢量场。

1)、将海面三角面元离散，并将不同的面元渲染不同的颜色，因为每个颜色分量取值为0-255，因此最多可以渲染的面元数为255*255*255≈1600W面元。

2)、将海面投影在于视线方向垂直的平面上，因为在OpenGL的绘图过程中，即使是后绘制的物体，但如果后绘制的物体在先绘制物体的后面，那么它将不会被绘制。在实际操作中，调用OpenGL函数，改变视线方向。

3)、读取屏幕中的颜色值，因为每一个面元都顺序分配了唯一的RGB颜色值，如果在读取的颜色值中有该面元的颜色值，那么该面元不被遮挡，反之则被其他面元遮挡。

步骤4，利用重力波调制张力波谱的面元法及场的矢量叠加原理求得电大区域时变多尺度海面的电磁场的时间序列，最终得到场的多普勒谱。

根据散射场的矢量叠加原理得到整个海面散射矢量场的时间序列，其公式为：

$E_{pq}^{total}({\hat{k}}_s,{\hat{k}}_i,t)=\underset{M}{\&Sigma;}\underset{N}{\&Sigma;}\&lsqb;E_{pq}^{SPM}({\hat{k}}_s,{\hat{k}}_i,t)+E_{pq}^{KA}({\hat{k}}_s,{\hat{k}}_i,t)\&rsqb;$

时变多尺度海面电磁散射场多普勒谱公式为：

$S_0\left(f\right)=\frac{1}{T}{\left|{\&Integral;}_0^T{E}_{PQ}^{scatt}(\widehat{k_i},\widehat{k_s},t)\exp \left(i2\mathrm{\&pi;}ft\right)dt\right|}^2$

步骤5，利用CUDA并行计算将粗糙海面不同面元进行并行，提高计算效率。

为了提高时变多尺度海面后向电磁散射时间序列建模的运算速度，我们通过在CPU-GPU异构平台上设计并实现了一个并行方案。在并行方案中，每一个海面中划分的海面面元为一个thread，最终对时间进行循环，对每个面元进行同时并行，得到单个面元对应的散射场。最后返还到主机端，利用场的矢量叠加原理，得到每个时间点的单个角度的散射场。

本发明通过重力波调制张力波谱的高频近似方法，可以很好的解决时变多尺度海面的电磁散射时间序列及多普勒谱特性问题，能够正确的实现不同海况条件下的电大区域海面电磁散射时间序列及多普勒谱的建模。但该算法非常耗时，海面上任意一个面元每一个时间点都需要进行大量的电磁散射运算。为了提高时变多尺度海面后向电磁散射时间序列建模的运算速度，我们通过在CPU-GPU异构平台上设计并实现了一个并行方案。在并行方案中，每一个海面中划分的海面面元为一个thread，最终对时间进行循环，对每个面元进行同时并行，得到单个面元对应的散射场。最后返还到主机端，利用场的矢量叠加原理，得到每个时间点的单个角度的散射场。

本发明不局限于上述具体的实施方式，本领域的普通技术人员从上述构思出发，不经过创造性的劳动，所作出的种种变换，均落在本发明的保护范围之内。

Claims (6)

1.一种时变多尺度电大区域海面电磁散射矢量场的并行计算方法，其特征在于，包括以下步骤：

1)利用叠加法基于海谱建立时变多尺度电大区域海面模型，根据不同的雷达参数，剖分为不同的面元，结合CUDA高性能并行计算方法，对海面每个点的高度进行并行建模，每个kernel对应每个点，提高建模效率；

2)基于大尺度重力波分量成立条件与小尺度张力波分量成立条件，给出不同频率不同海情大尺度重力波与小尺度张力波的截断波数k_c；

3)基于截断波数k_c，当面元为大尺度重力波时，利用面元基尔霍夫方法求得相干散射场，当面元为小尺度张力波时，利用面元小扰动近似方法得到基于大尺度面元本地坐标系下的张力波Bragg散射场，并利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断；

4)利用重力波调制张力波谱的面元法及场的矢量叠加原理求得电大区域时变多尺度海面的电磁场的时间序列，最终得到场的多普勒谱；

5)基于GPU利用CUDA对每个海面的面元进行并行，每个kernel计算单个面元的散射场，同时进行计算，并对时间循环，利用异步传输理论，提高并行效率，大大提高计算效率。

2.根据权利要求1所述的时变多尺度电大区域海面电磁散射矢量场的并行计算方法，其特征在于，所述的步骤1)中给出基于叠加法的时变多尺度海面并行建模，时变多尺度海面大尺度重力波的建模，并对其进行并行计算提高效率：

S_1d(k)＝(B_L+B_H)/k³,S_dir＝[1+Δ(k)cos(2φ)]/2π\*MERGEFORMAT(3)

其中，S_1d(k)为一维Efouhaily谱，S_dir为方向谱，B_L为长波曲率谱，B_H为短波曲率谱，k对应海面的波数，ω为角频率，ω²＝gk，ε_ij为海面的初始相位，在0～2π内满足均匀分布，φ为雷达波与风向的夹角；M，N代表频率和方向角的采样点数；基于线性叠加法的二维海面的模拟，本质就是对海面剖分网格计算，海面剖分网格越密，海面的采样点数就越多，那么得到的海面就更接近真实海面，但是这意味着更大的计算量，这对于采用串行方案进行二维海面模拟是非常耗时的。在我们的并行方案中，充分利用二维海面仿真的内在并行性，采用粗粒度线程的并行策略，即一个GPU线程完成一个采样点的计算，这就意味着我们可以通过更加精细的剖分海面获得更加逼真的海面，但计算量并不会显著增加。

3.根据权利要求1所述的时变多尺度电大区域海面电磁散射矢量场的并行计算方法，其特征在于，步骤2)利用基尔霍夫方法与小扰动近似方法，求得不同风速，不同雷达波频率下对应的截断频率k_c。

电大区域海面不同区域对应不同散射机理，大尺度重力波对应小角度镜面反射，小尺度张力波对应张力波的布拉格散射。随着海情增大，出现泡沫覆盖区域，出现粒子与粗糙面耦合散射机理。组成海面的面元在不同的风速与不同的雷达频率照射下对应不同的雷达频率，因此需将得到大尺度重力波与小尺度张力波的截断频率k_c，粗糙海面模型中针对大尺度杂波分量模型的成立的条件为：

其中，k_i为入射电磁波数，θ_i为入射角度，σ_c是海面平均曲率的标准离差，δ_l为大尺度重力波的均方根高度，S(k)为海谱函数。

复合表面模型中，小尺度波分量成立条件为：

k_iδ_scos(θi)＜＜1\*MERGEFORMAT(6)

δ_s为小尺度波面的均方根高度。

4.根据权利要求1所述的时变多尺度电大区域海面电磁散射矢量场的并行计算方法，其特征在于，将步骤3)基于入射频率与风速求得截断频率k_c，面元为大尺度重力波时，利用面元基尔霍夫方法求得相干散射场：

其中，Z_x和Z_y为面元沿x方向与y方向上的斜率。(x₀,y₀,z₀)为每个面元的中心，为极化幅度系数，为入射极化波数矢量，为散射极化波数矢量，p,q＝h,v为入射波与散射波的极化方式。当面元为小尺度张力波时，利用面元小扰动近似方法得到基于大尺度面元本地坐标系下的张力波Bragg散射场：

其中，ε为粗糙海面相对介电常数，为海面单个面元的法线矢量，为本地坐标系下的极化系数，为毛细波幅度，k_c为Bragg波长对应的波数，ΔS为面元大小，J_n(x)为n阶贝塞尔函数；

利用Z-Buffer方法对不同雷达照射下海面进行面元遮挡判断。

5.根据权利要求1所述的时变多尺度电大区域海面电磁散射矢量场的并行计算方法，其特征在于，步骤4)利用场的矢量叠加原理得到散射场的时间序列，并得到多普勒谱：

6.根据权利要求1所述的基于重力波调制张力波谱的时变多尺度电大区域海面电磁散射矢量场时间序列的并行计算方法，其特征在于，将步骤5)利用CUDA并行计算将粗糙海面不同面元进行并行，提高计算效率。