本发明涉及高频地波雷达遥感探测技术领域，尤其涉及一种基于三维正压浅海大陆架模型的深层海流反演方法。
背景技术：
为了获得海洋表面以下海流的运动信息，现有技术采用直接接触测量水体的方式，所以获得大面积水域的数据则需要投放相当数量的直接测量仪器。该方法对周围测量环境的要求高，成本高，探测效率低。高频地波雷达遥感技术是近三十年来发展起来的探测海洋表面动力学参数的技术，覆盖范围广，探测精度和效率都较高。但是由于电磁波在水体中的衰减严重，使得高频地波雷达遥感探测技术仅限于探测海洋水体表面，不能直接获得海洋内部水体运动的相关信息是高频地波雷达遥感技术在测量海洋流参数方面的不足。
技术实现要素：
针对现有技术的不足，本发明所解决的问题是利用高频地波雷达探测得到的海洋表面流参数信息即可得到深层流参数信息。
为解决上述技术问题，本发明采用的技术方案是一种基于三维正压浅海大陆架模型的深层海流反演方法，采用归一化的坐标系的三维非线性正压浅海大陆架模式，利用高频地波雷达探测到的海洋表面流数据，考虑复杂地形的边界问题，参数化物理过程，简化海洋控制方程组，再建立三维网格，采用差分法求解离散方程组进而反演得到相应的深层流数据信息，包括如下步骤：
(1)选择σ坐标系；
假设在高频地波雷达所探测的每场数据的累积时间内，海水在垂直方向上处于静力学平衡状态；海水处于正压模式，海水不可压缩，海水的运动速度不随时间变化而变化；为了便于分层，采用归一化的σ坐标系，使水浅处与水深处海水都具有相同的垂向分层数，并且每一层的侧边界不再变化，从海底到海面的高度处处规范化为1；
在σ坐标系下x轴和y轴不变，用σ取代z轴，其中z和σ的转换关系如下：
坐标变换后，设总的海水深度为h＝h+ζ，海洋表面z＝ζ,σ＝0，在海底z＝-h,σ＝-1；
(2)导出归一化坐标系下的控制方程组；
在σ坐标系中，利用dσ和dz之间的微分系数关系
1)连续方程表达式如下：
2)方向动量方程的非线性项采用能量二次守恒的半动量形式，表达式如下：
3)流体静力学方程如下：
表达式中的u、v、w分别代表流速在x、y、σ方向上的分量；h＝h+ζ是总海深；ζ是静止海面状态下的浪高；h是静态水域深度；ρ是海水密度；r是地球半径；
由公式(02)得到垂直流速的表达式如下：
在动量方程中，水平压强梯度力方程如下：
在实际观测资料中发现，海水水平密度变化较小，所以在计算中可以忽略水平压强梯度力第二项，简化的水平压强梯度力方程如下：
(3)参数化物理过程；
海洋表面限制条件由高频地波雷达监测的表面流参数值决定，设定海洋表面垂直速度w＝0，动力学参数限制条件如下：
其中τa表面风应力矢量，τb是海底摩擦力矢量，表达式如下所示：
τa＝ρaca|va|va(13)
τb＝ρbcb|vb|vb(14)
其中ρa是空气密度；va是风速矢量；ca＝(0.73+0.069|va|)×10-3是海水对风的拖曳系数；ρb是海水密度；vb是靠近海底的流速矢量；
其中
(4)建立三维网格，离散化方程组；
采用以下方式分别将u、v、w在x、y、σ方向上进行离散，以u为例：
u的值为同一网格中四个结点的加权平均值，θ是空间加权系数，当θ＜0.5时，上述表达式为一阶，当θ＝0.5时，上述表达式为二阶，采用二阶精度，其格式如下：
(5)迭代法求解方程组得出深层流信息；
在公式(18)中，后一层结点上的值不仅与上一层结点值有关，还和本层已计算出的结点值有关，利用计算机，运用迭代法求解离散化的海洋控制方程组；
(6)实测数据与反演结果对比分析；
分析数值反演角度图发现，随着海流深度的增加，实测海流的流向和反演得出的海流的流向都呈逐渐右偏的趋势；分析数值反演速度图发现，随着海流深度的增加，实测海流的流速和反演得出的海流的流速都呈逐渐减小的趋势，在靠近海底处，由于底部摩擦力的作用，流速都趋近于零。
本发明的有益效果：对比试验结果发现，实测数据曲线图和反演结果曲线图的整体趋势接近，反演曲线可以反映出海洋深层流变化的整体趋势，符合海洋动力学运动规律。
附图说明
图1是本发明流程示意图；
图2是坐标变换前的坐标系示意图；
图3是坐标变换后的坐标系示意图；
图4是探测点在14:10的角度变化图；
图5是探测点在14:30的角度变化图；
图6是探测点在14:10的速度变化图；
图7是探测点在14:30的速度变化图。
具体实施方式
下面结合附图对本发明的具体实施方式作进一步的说明，但不是对本发明的限定。
图1示出了一种基于三维正压浅海大陆架模型的深层海流反演方法，采用归一化的坐标系的三维非线性正压浅海大陆架模式，利用高频地波雷达探测到的海洋表面流数据，考虑复杂地形的边界问题，参数化物理过程，简化海洋控制方程组，再建立三维网格，采用差分法求解离散方程组进而反演得到相应的深层流数据信息，包括如下步骤：
(1)选择σ坐标系；
假设在高频地波雷达所探测的每场数据的累积时间内，海水在垂直方向上处于静力学平衡状态；海水处于正压模式，海水不可压缩，海水的运动速度不随时间变化而变化；为了便于分层，把如图2所示的变换前的坐标系变换为归一化的σ坐标系，如图3所示，使水浅处与水深处海水都具有相同的垂向分层数，并且每一层的侧边界不再变化，从海底到海面的高度处处规范化为1；
在σ坐标系下x轴和y轴不变，用σ取代z轴，其中z和σ的转换关系如下：
坐标变换后，设总的海水深度为h＝h+ζ，海洋表面z＝ζ,σ＝0，在海底z＝-h,σ＝-1；
(2)导出归一化坐标系下的控制方程组；
在σ坐标系中，利用dσ和dz之间的微分系数关系
1)连续方程表达式如下：
2)方向动量方程的非线性项采用能量二次守恒的半动量形式，表达式如下：
3)流体静力学方程如下：
表达式中的u、v、w分别代表流速在x、y、σ方向上的分量；h＝h+ζ是总海深；ζ是静止海面状态下的浪高；h是静态水域深度；ρ是海水密度；r是地球半径；
由公式(02)得到垂直流速的表达式如下：
在动量方程中，水平压强梯度力方程如下：
在实际观测资料中发现，海水水平密度变化较小，所以在计算中可以忽略水平压强梯度力第二项，简化的水平压强梯度力方程如下：
(3)参数化物理过程；
海洋表面限制条件由高频地波雷达监测的表面流参数值决定，设定海洋表面垂直速度w＝0，动力学参数限制条件如下：
其中τa表面风应力矢量，τb是海底摩擦力矢量，表达式如下所示：
τa＝ρaca|va|va(13)
τb＝ρbcb|vb|vb(14)
其中ρa是空气密度；va是风速矢量；ca＝(0.73+0.069|va|)×10-3是海水对风的拖曳系数；ρb是海水密度；vb是靠近海底的流速矢量；
其中
(4)建立三维网格，离散化方程组；
采用以下方式分别将u、v、w在x、y、σ方向上进行离散，以u为例：
u的值为同一网格中四个结点的加权平均值，θ是空间加权系数，当θ＜0.5时，上述表达式为一阶，当θ＝0.5时，上述表达式为二阶，采用二阶精度，其格式如下：
(5)迭代法求解方程组得出深层流信息；
在公式(18)中，后一层结点上的值不仅与上一层结点值有关，还和本层已计算出的结点值有关，利用计算机，运用迭代法求解离散化的海洋控制方程组；
(6)实测数据与反演结果对比分析；
实验结果图中实线代表海流计的实测数据，虚线代表利用本发明方法反演得出的结果；
图4和图5分别为探测点在14:10和14:30两个不同时刻数值反演的角度图。分析数值反演角度图发现，随着海流深度的增加，实测海流的流向和反演得出的海流的流向都呈逐渐右偏的趋势；图6和图7分别为探测点在14:10和14:30两个不同时刻数值反演的速度图。分析数值反演速度图发现，随着海流深度的增加，实测海流的流速和反演得出的海流的流速都呈逐渐减小的趋势，在靠近海底处，由于底部摩擦力的作用，流速都趋近于零；
随着深度的变化，反演海流的整个运动趋势是符合海洋动力学规律的；虽然实测海流数据的流向变化图和流速变化图在小幅度上都出现了反复波动，但并不影响随着深度增加，流向整体右偏和流速逐渐减小的整体趋势；
综上，实测数据曲线图和反演结果曲线图的整体趋势接近，反演曲线可以反映出海洋深层流变化的整体趋势，符合海洋动力学运动规律。
以上结合附图对本发明的实施方式做出了详细说明，但本发明不局限于所描述的实施方式。对于本领域技术人员而言，在不脱离本发明的原理和精神的情况下，对这些实施方式进行各种变化、修改、替换和变型仍落入本发明的保护范围内。