流函数算法原理以及实现（在海洋矢量流方面）

流函数算法原理

当海水处于流体静力学平衡，二维的海水流动可以用一系列流线来表征, 流线是具有相同的流函数值的点的轨迹。忽略海水运动在垂直方向上的变化，且认为表层海水密度均匀。在海平面上建立直角坐标系，设x 轴正向朝东， y 轴正向朝北，二维海水运动连续性方程为

\[\nabla\cdot V=0\]

其中V为二维速度矢量，u，v分别为其在x，y方向上的速度分量。根据上面的二维海水运动连续性方程引入流函数（x，y）， 定义为

\[u=-\frac{\partial\psi}{\partial y}\]

\[v=\frac{\partial\psi}{\partial x}\]

并设在边界出\(\psi=0\)

高频地波雷达探测的海流径向分量\(V_{r}\)，用u，v分量表示为：

\[ucos\theta+vsin\theta=V_{r}\]

其中\(V_{r}\)是高频雷达探测到的径向速度，\(\theta\)是极坐标系中的相应方位角。在直角坐标系中可以将上式表示为：

\[-x\frac{\partial\psi}{\partial y}+y\frac{\partial\psi}{\partial x}=\sqrt{x^{2}+y^{2}}V_{r}\]

求解这个偏微分方程就可以得到流函数的表达式。用最小二乘法来求解一个关于流函数泰勒系数的超定方程组。

流函数实现

单站高频雷达探测到的径向流场是等距离等角度离散地分布在雷达覆盖区域中的，因此上述方程的解析解是很困难的，而较易得到数值近似解。利用泰勒展开法来求解线性偏微分方程的近似解。这种方法主要是在一个确定的点周围用系数待定的泰勒序列代替流函数，而这个点附近的径向流场\(V_{r}(x,y)\) 作为已知函数代入方程中通过最小二乘法来求未知系数。这样偏微分方程就转化成关于泰勒系数的线性方程组。在一定的雷达区域中将流函数表示为N阶泰勒多项式:

将其带到原理中所论述的直角坐标系方程得到线性代数方程为

该方程中共有\(N^{2}+2N\)个未知系数，因而至少需要\(N^{2}+2N\)个方程才能确定流函数的表达式。考虑到径向流的探测误差, 要适当增大局部区域的大小使得径向流数目远大于这个值, 从而形成一个超定方程组; 利用最小二乘法来求解该方程组. 设该区域中共有M（\(M>N^{2}+2N\)）个径向流数据，组成一个向量\(V_{r}=(V_{r1},V_{r2}…V_{rm})^{T}\)，相应的直角坐标向量为\(x=(x1,x2,….,xm)^{T}\)和\(y=(y1,y2….ym)^{T}\)。可以得到下面的方程组：

其最小二乘解由下列标准方程得到

\[G^{T}G_{\alpha}=G^{T}V_{r}\]

其中G为\(M\times(N^{2}+2N)\)的系数矩阵，\(\alpha\)是待求系数组成的\(N^{2}+2N\)元列向量。最小二乘法求得流函数的近似表达式后，根据该区域中两速度分量表达式, 二者合成得到矢量流场的大小V和方向
\(\xi\)分别为：

\[V=\sqrt{u^{2}+v^{2}}\]

\[\xi=tan^{-1}(v/u)\]