CN102693355A - 稀少验潮站数据控制的高精度水位推算技术 - Google Patents

Abstract

本发明涉及一种稀少验潮站数据控制的高精度水位推算技术，首先利用资料同化技术同化验潮站数据，获取高精度潮汐调和常数，通过建立潮汐预报模型，计算瞬时潮高，然后由验潮站观测数据检测非潮汐水位，并根据余水位良好的空间一致性特点，将非潮汐水位推估至所需水位计算点，实现实际水位的高精度恢复与再现，达到为海岛礁测量提供精确瞬时水位信息的目的。

Description

稀少验潮站数据控制的高精度水位推算技术

技术领域

本发明属于海洋监测领域，具体属于一种稀少验潮站潮汐预报及水位推算方法。 

背景技术

潮汐预报通过潮汐参数仅能计算潮汐作用驱动的海面高度，事实上，海面的实际升降虽主体上变现为潮汐变化，但气象等因素引起的短期信号性质的变化也具有不可忽略的量级。本方法首先利用资料同化技术同化验潮站数据，获取高精度潮汐调和常数，通过建立潮汐预报模型，计算瞬时潮高，然后由验潮站观测数据检测非潮汐水位，并根据余水位良好的空间一致性特点，将非潮汐水位推估至所需水位计算点，实现实际水位的高精度恢复与再现，达到为海岛礁测量提供精确瞬时水位信息的目的。 

发明内容

为了实现稀少验潮站控制下的水位推算，包括数值模拟、资料同化和潮汐快速预报步骤。本方法主要借助数值模拟方法构建目标海域潮波数值模式，利用资料同化技术提高模式模拟精度，模拟得到各主要分潮调和常数，然后利用潮汐快速预报方法，实现潮汐预报及水位方法推算。 

附图说明：

图1潮站数据控制的高精度水位推算流程图。 

实施例：

(一)数值模拟方法 

本方法利用σ坐标下的水动力数值模式构建目标海域潮波数值模式。模式的控制方程可表达为： 

$\frac{\∂\mathrm{DU}}{\∂x}+\frac{\∂\mathrm{DV}}{\∂y}+\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}+\frac{\∂\mathrm{\ζ}}{\∂x}=0,---\left(1\right)$

$\frac{\∂\mathrm{UD}}{\∂t}+\frac{\∂U^{2}D}{\∂x}+\frac{\∂\mathrm{UVD}}{\∂y}+\frac{\∂\mathrm{U\ω}}{\∂\mathrm{\σ}}-\mathrm{fVD}+\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂x}$

$+\frac{{\mathrm{gD}}^2}{{\mathrm{\ρ}}_0}{\∫}_{\mathrm{\σ}}^0[\frac{{\∂\mathrm{\ρ}}^{\′}}{\∂x}-\frac{{\mathrm{\σ}}^{\′}}{D}\frac{\∂D}{\∂x}\frac{\∂{\mathrm{\ρ}}^{\′}}{\∂{\mathrm{\σ}}^{\′}}]d{\mathrm{\σ}}^{\′}=\frac{\∂}{\∂\mathrm{\σ}}\left[\frac{K_M}{D}\frac{\∂U}{\∂\mathrm{\σ}}\right]+F_x,---\left(2\right)$

$\frac{\∂\mathrm{VD}}{\∂t}+\frac{\∂\mathrm{UVD}}{\∂x}+\frac{\∂V^{2}D}{\∂y}+\frac{\∂\mathrm{V\ω}}{\∂\mathrm{\σ}}+\mathrm{fUD}+\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂y}$

$+\frac{{\mathrm{gD}}^2}{{\mathrm{\ρ}}_0}{\∫}_{\mathrm{\σ}}^0[\frac{{\∂\mathrm{\ρ}}^{\′}}{\∂y}-\frac{{\mathrm{\σ}}^{\′}}{D}\frac{\∂D}{\∂y}\frac{\∂{\mathrm{\ρ}}^{\′}}{\∂{\mathrm{\σ}}^{\′}}]d{\mathrm{\σ}}^{\′}=\frac{\∂}{\∂\mathrm{\σ}}\left[\frac{K_M}{D}\frac{\∂U}{\∂\mathrm{\σ}}\right]+F_y,---\left(3\right)$

为了使本方法适用于深水区域，在潮波模式中考虑了引潮力。为考虑引潮力的作用，将模式控制方程(2)和(3)中的ζ应改写为ζ′，后者考虑引潮势、体潮、负载潮及自吸效应。其原理阐述如下。 

卫星高度计观测得到的海面高度是海面与地球参考椭球面的距离，称作地心潮(geocentric tide)，记为ζ_g，它等于通常的海潮ζ和固体地球潮之和；而后者又由体潮(body tide)ζ_b和海洋负载潮(ocean-load tide)ζ_l之和所组成。因而有(Ray，1998；黄祖珂、黄磊，2005) 

ζ_g＝ζ+ζ_b+ζ_l.      (4) 

这些潮汐中的任何一项均可展开为球调和级数，例如， 

$\mathrm{\ζ}=\underset{n}{\mathrm{\Σ}}{\mathrm{\ζ}}_n.---\left(5\right)$

体潮对引潮力的响应近似为平衡态，故有 

${\mathrm{\ζ}}_b=\underset{n}{\mathrm{\Σ}}{h}_n{\mathrm{\Φ}}_n/g,---\left(6\right)$

其中Φ_n为引潮势，h_n为Love数。由于主要的半日和全日分潮均属于n＝2的二阶球调和项，故可取 

ζ_b＝h₂Φ₂/g.     (7) 

类似地，负载潮ζ_l亦可写作 

${\mathrm{\ζ}}_l=\underset{n}{\mathrm{\Σ}}{h_n}^{\′}{U}_n/g,---\left(8\right)$

其中，h_n′为负载Love数，U_n为ζ_n引起的重力势： 

U_n＝3g(ρ_w/ρ_e)ζ_n/(n+1)，    (9) 

式中ρ_w和ρ_e分别为海水和地球平均密度。式(8)也写作 

${\mathrm{\ζ}}_l=\underset{n}{\mathrm{\Σ}}{h^{\′}}_n{\mathrm{\α}}_n{\mathrm{\ζ}}_n,---\left(10\right)$

其中 

α_n＝3(ρ_w/ρ_e)/(n+1).    (11) 

体潮、海潮和负载潮均引起重力势的改变，故总的作用在海水上的重力势Γ为 

$\mathrm{\Γ}=(1+k_2){\mathrm{\Φ}}_2+\underset{n}{\mathrm{\Σ}}(1+{k^{\′}}_n)U_n$

$=(1+k_2){\mathrm{\Φ}}_2+\underset{n}{\mathrm{\Σ}}(1+{k^{\′}}_n)g{\mathrm{\α}}_n{\mathrm{\ζ}}_n.---\left(12\right)$

因而可得 

${\mathrm{\ζ}}^{\′}={\mathrm{\ζ}}_g-\mathrm{\Γ}/g$

$=\mathrm{\ζ}-(1+k_2-h_2){\mathrm{\Φ}}_2/g-\underset{n}{\mathrm{\Σ}}(1+{k^{\′}}_n-{h^{\′}}_n){\mathrm{\α}}_n{\mathrm{\ζ}}_n.---\left(13\right)$

在海洋潮汐研究中常令 

$\stackrel{\‾}{\mathrm{\ζ}}=(1+k_2-h_2){\mathrm{\Φ}}_2/g,---\left(14\right)$

及 

${\mathrm{\ζ}}_{\mathrm{SAL}}=\underset{n}{\mathrm{\Σ}}(1+{k^{\′}}_n-{h^{\′}}_n){\mathrm{\α}}_n{\mathrm{\ζ}}_n.---\left(15\right)$

称为修正(即考虑体潮效应)后的平衡潮高，ζ_SAL称为自吸负载(self-attraction/load)项。于是， 

${\mathrm{\ζ}}^{\′}=\mathrm{\ζ}-\stackrel{\‾}{\mathrm{\ζ}}-{\mathrm{\ζ}}_{\mathrm{SAL}}.---\left(16\right)$

我国近海早期的潮波数值模式中大部分都已考虑了 

项。方国洪等(1994)对南海潮波的模拟表明， 项对南海M₂的影响可达20％左右。但自吸负载项的作用在以往的模式中均未考虑。为了提高本方法的适用范围，需进一步考虑自吸负载项。 

自吸负载潮在大洋潮波数值模拟中常作以下近似： 

ζ_SAL≈bζ，  (17) 

其中b取作常数，因而ζ-ζ_SAL写作 

ζ-ζ_SAL≈βζ，     (18) 

其中β＝1-b。实际上，取b(或β)为常数的作法在大洋效果较好，近海则不应取作常数。本方法给出b为水深h的函数，表示为： 

b＝a₀+a₁lnh+a₂(lnh)²，  (19) 

其中a₀，a₁，a₂为待定常数。 

根据Ray(1999)，近似有 

ζ_SAL≈-1.8ζ_l.     (20) 

从而可以先计算ζ_l和ζ之间的关系： 

ζ_l≈-b′ζ，     (21) 

b′＝a′₀+a′₁lnh+a′₂(lnh)²，  (22) 

其中 

(b′，a′₀，a′₁，a′₂)＝(b，a₀，a₁，a₂)/1.8     (23) 

按式(21)， 

b′＝-(H_l/H)cos(g_l-g)，   (24) 

式中(H_l，g_l)和(H，g)分别为负载潮ζ_l和海潮ζ的振幅和迟角。 

(二)资料同化方法 

在同化方法方面，海底地形的底摩擦效应对潮波运动有着至关重要的作用，在潮汐潮流的数值计算中，底摩擦效应主要由底摩擦系数来刻画，例如潮能底边界耗散的计算公式为：d＝C_dρ<u³>(W/m²)，其中C_d为底摩擦系数。本技术 方案将伴随方法应用于示范区潮波数值模拟，利用验潮站资料，反演得到底摩擦系数。 

同化模式采用的潮汐模型为： 

$\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;t}+\frac{\&PartialD;\left[(h+\mathrm{\&zeta;})u\right]}{\&PartialD;x}+\frac{\&PartialD;\left[(h+\mathrm{\&zeta;})v\right]}{\&PartialD;y}=0---\left(25\right)$

$\frac{\&PartialD;u}{\&PartialD;t}+u\frac{\&PartialD;u}{\&PartialD;x}+v\frac{\&PartialD;u}{\&PartialD;y}-\mathrm{fv}+g\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;x}+\frac{c_d\sqrt{u^2+v^2}u}{h+\mathrm{\&zeta;}}-A(\frac{{\&PartialD;}^{2}u}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}u}{\&PartialD;y^2})=0---\left(26\right)$

$\frac{\&PartialD;v}{\&PartialD;t}+u\frac{\&PartialD;v}{\&PartialD;x}+v\frac{\&PartialD;v}{\&PartialD;y}+\mathrm{fu}+g\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;y}+\frac{c_d\sqrt{u^2+v^2}v}{h+\mathrm{\&zeta;}}-A(\frac{{\&PartialD;}^{2}v}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}v}{\&PartialD;y^2})=0---\left(27\right)$

其中，t代表时间；x，y是Cartesian坐标，分别取向东和向北为正；u(x，y，t)和v(x，y，t)分别代表x，y方向的流速分量；h(x，y)表示(x，y)处的静水深度，ζ(x，y，t)表示未扰动海面以上的水位高度，H＝h+ζ为总水深；g是重力加速度； 

是Coriolis参数，Ω为地球自转角速度， 

为地理纬度；c_d为底摩擦系数，A为水平涡动粘性系数。 

计算网格采用Arakawa C网格，水位ζ位于网格的中心点，节点坐标为(i，j)，u的节点坐标为 

v的节点坐标为 

水深h和水位ζ取在同一点上。 

伴随同化方法以经典的拉格朗日乘子法为理论基础，将有约束最小值问题转化为无约束最小值问题。 

首先构造代价函数： 

$J\left(\mathrm{\&zeta;}\right)=\frac{1}{2}{K}_{\mathrm{\&zeta;}}{\&Integral;}_{\mathrm{\&Sigma;}}{(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_{\mathrm{obs}})}^2\mathrm{d\&sigma;}---\left(28\right)$

其中，ζ_obs为观测值，K_ζ为权重系数，有观测值的点K_ζ取为1，否则取为∑代表整个时间和空间的集合。 

再构造拉格朗日函数： 

$L(\mathrm{\&zeta;},u,v,\mathrm{\&lambda;},\mathrm{\&mu;},\mathrm{\&upsi;})=J\left(\mathrm{\&zeta;}\right)+{\&Integral;}_{\mathrm{\&Sigma;}}[\mathrm{\&lambda;}(\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;t}+\frac{\&PartialD;\left[(h+\mathrm{\&zeta;})u\right]}{\&PartialD;x}+\frac{\&PartialD;\left[(h+\mathrm{\&zeta;})v\right]}{\&PartialD;y})$

$+\mathrm{\&mu;}(\frac{\&PartialD;u}{\&PartialD;t}+u\frac{\&PartialD;u}{\&PartialD;x}+v\frac{\&PartialD;u}{\&PartialD;y}-\mathrm{fv}+g\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;x}+\frac{c_d\sqrt{u^2+v^2}u}{h+\mathrm{\&zeta;}}-A(\frac{{\&PartialD;}^{2}u}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}u}{\&PartialD;y^2}))$

$+\mathrm{\&upsi;}(\frac{\&PartialD;v}{\&PartialD;t}+u\frac{\&PartialD;v}{\&PartialD;x}+v\frac{\&PartialD;v}{\&PartialD;y}+\mathrm{fu}+g\frac{\&PartialD;\mathrm{\&zeta;}}{\&PartialD;y}+\frac{c_d\sqrt{u^2+v^2}v}{h+\mathrm{\&zeta;}}-A(\frac{{\&PartialD;}^{2}v}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}v}{\&PartialD;y^2}))]\mathrm{d\&sigma;}---\left(29\right)$

其中，λ，μ，υ分别为ζ，u，v的伴随变量。 

为了使代价函数达到最小，我们要求： 

$\frac{\&PartialD;L}{\&PartialD;\mathrm{\&lambda;}}=0,$ $\frac{\&PartialD;L}{\&PartialD;\mathrm{\&mu;}}=0,$ $\frac{\&PartialD;L}{\&PartialD;\mathrm{\&upsi;}}=0,---\left(30\right)$

$\frac{\&PartialD;L}{\&PartialD;\mathrm{\&zeta;}}=0,$ $\frac{\&PartialD;L}{\&PartialD;u}=0,$ $\frac{\&PartialD;L}{\&PartialD;v}=0,---\left(31\right)$

$\frac{\&PartialD;L}{{\&PartialD;c}_d}=0,$ $\frac{\&PartialD;L}{\&PartialD;A}=0---\left(32\right)$

经过推导，得到的是原来的模型方程，并进一步得到伴随方程： 

$\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;t}+g\frac{\&PartialD;\mathrm{\&mu;}}{\&PartialD;x}+g\frac{\&PartialD;\mathrm{\&upsi;}}{\&PartialD;y}+u\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;x}+v\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;y}+\frac{c_d\sqrt{u^2+v^2}}{{(h+\mathrm{\&zeta;})}^2}(\mathrm{\&mu;u}+\mathrm{\&upsi;v})=K_{\mathrm{\&zeta;}}(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_{\mathrm{obs}})---\left(33\right)$

$\frac{\&PartialD;\mathrm{\&mu;}}{\&PartialD;t}-\mathrm{\&mu;}\frac{\&PartialD;u}{\&PartialD;x}-\mathrm{\&upsi;}\frac{\&PartialD;v}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&mu;u}}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&mu;v}}{\&PartialD;y}-(f+\frac{c_d\mathrm{uv}}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}})\mathrm{\&upsi;}$

$+(h+\mathrm{\&zeta;})\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;x}-c_d\frac{{2u}^2+v^2}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}}\mathrm{\&mu;}+A(\frac{{\&PartialD;}^2\mathrm{\&mu;}}{\&PartialD;x^2}+\frac{{\&PartialD;}^2\mathrm{\&mu;}}{\&PartialD;y^2})=0---\left(34\right)$

$\frac{\&PartialD;\mathrm{\&upsi;}}{\&PartialD;t}-\mathrm{\&mu;}\frac{\&PartialD;u}{\&PartialD;y}-\mathrm{\&upsi;}\frac{\&PartialD;v}{\&PartialD;y}+\frac{\&PartialD;\mathrm{\&upsi;u}}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&upsi;v}}{\&PartialD;y}+(f-\frac{c_d\mathrm{uv}}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}})\mathrm{\&mu;}$

$+(h+\mathrm{\&zeta;})\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;y}-c_d\frac{u^2+2v^2}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}}\mathrm{\&upsi;}+A(\frac{{\&PartialD;}^2\mathrm{\&upsi;}}{\&PartialD;x^2}+\frac{{\&PartialD;}^2\mathrm{\&upsi;}}{\&PartialD;y^2})=0---\left(35\right)$

代价函数关于底摩擦系数和水平涡动粘性系数的梯度关系式为： 

$\frac{\&PartialD;J}{\&PartialD;c_d}+{\&Integral;}_{\mathrm{\&Sigma;}}[\mathrm{\&mu;}\frac{\sqrt{u^2+v^2}u}{h+\mathrm{\&zeta;}}+\mathrm{\&upsi;}\frac{\sqrt{u^2+v^2}v}{h+\mathrm{\&zeta;}}]\mathrm{d\&sigma;}=0---\left(36\right)$

$\frac{\&PartialD;J}{\&PartialD;A}-{\&Integral;}_{\mathrm{\&Sigma;}}[\mathrm{\&mu;A}(\frac{{\&PartialD;}^{2}u}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}u}{\&PartialD;y^2})+\mathrm{\&upsi;A}(\frac{{\&PartialD;}^{2}v}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}v}{\&PartialD;y^2})]\mathrm{d\&sigma;}=0---\left(37\right)$

用来校正潮波模型中的底摩擦系数和水平涡动粘性系数： 

$c_{\mathrm{dnew}}=c_{\mathrm{dold}}-{\mathrm{\&alpha;}}_{c_d}\frac{\&PartialD;J}{\&PartialD;c_d}---\left(38\right)$

$A_{\mathrm{new}}=A_{\mathrm{old}}-{\mathrm{\&alpha;}}_A\frac{\&PartialD;J}{\&PartialD;A}---\left(39\right)$

这里，α为下降步长。 

在特定海区给定初猜场的情况下，通过伴随反演及迭代过程使得代价函数降低并趋于稳定，最终得到该海区的底摩擦系数和水平涡动粘性系数。 

另外，本方法同时也利用趋近法来同化验潮站资料。此方法在连续方程(1)右端加上一个强迫项F_ζ： 

$F_{\mathrm{\&zeta;}}=\mathrm{\&alpha;}\underset{m}{\mathrm{\&Sigma;}}{w}_m(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_{\mathrm{obs},m}).---\left(40\right)$

这里ζ_obs，m为观测值，实际上是由观测所得调和常数进行预报的nΔt时刻数值： 

${\mathrm{\&zeta;}}_{\mathrm{obs},m}^n=H_m\cos (\mathrm{\&omega;n\&Delta;t}-g_m),---\left(41\right)$

其中m代表观测点；n代表时间步。 

对任一内点(i，j)，对应的经纬度为(λ_i，j，φ_i，j)，w_m值取为： 

其中， 

式中D取1/6°，D_max取1°，a为可调整参数，本研究的α＝0.0014s^-1。 

(三)潮汐预报算方法 

本方法设计了一种较快的潮汐预报算法。借助数值模拟技术构建潮波数值模式和资料同化技术模拟得到的各主要分潮调和常数，实现潮汐预报及水位推算。相对于平均海面的潮高h可表示为 

$h\left(t\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{f}_l\left(t\right)H_l\cos [{\mathrm{\&omega;}}_{l}t+v_{0,l}+u_l\left(t\right)-g_l]---\left(44\right)$

其中，H_l，g_l分别为第l个分潮的振幅和迟角，L是分潮个数；ω_l是分潮角速率；v_o，l为平衡潮第l个分潮在t＝0时刻的位相。这里f，ωt及v₀+u与时间有关，与地点无关；而H和g与地点有关，与时间无关。故上式可写作 

$h\left(t\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}[{\mathrm{\&alpha;}}_l\left(t\right)X_l+{\mathrm{\&beta;}}_l\left(t\right)Y_l],---\left(45\right)$

其中 

$\left\{\begin{array}{c}{\mathrm{\&alpha;}}_l\left(t\right)=f_l\left(t\right)\cos [{\mathrm{\&omega;}}_{l}t+v_{o,l}+u_l\left(t\right)],\\ {\mathrm{\&beta;}}_l\left(t\right)=f_l\left(t\right)\sin [{\mathrm{\&omega;}}_{l}t+v_{o,l}+u_l\left(t\right)];\end{array}\right.---\left(46\right)$

$\left\{\begin{array}{c}{X}_l=H_l\cos g_l,\\ Y_l=H_l\sin g_l.\end{array}\right.---\left(47\right)$

实际计算时不必对每个站点(层次)都计算α和β，可以一次性算出50年逐时的α，β值储存起来。对于每个站点(层次)，只要用该处的X，Y值代入式(45)即可。在本方法中，对潮位我们采用M₂，S₂，N₂，K₂，K₁，O₁，P₁，Q₁，M₄，MS₄，M₆，Sa和Ss_a 13个分潮，即L＝13。 

(四)水位推算 

海面水位变化主要由潮位和余水位来控制，对于目标站位(x，y)在t时刻的瞬时海面高可以表示为： 

h(x，y，t)＝T(x，y，t)+S(x，y，t)      (48) 

其中目标站位在t时刻的瞬时潮位T(x，y，t)可以通过潮汐预报得到，而余水位S(x，y，t)则需要利用验潮站资料推算获得。首先需要分析附近验潮站水位资料，提取非潮汐水位，即余水位，然后根据目标站位与附近验潮站水位的拟合关系，推算目标站位的余水位。基于余水位的水位改正方法可有效补偿非潮汐因素影响，从而有效提高水位推算精度。 

Claims (4)

1.一种稀少验潮站数据控制的高精度水位推算技术，其包括潮波数值模拟、资料同化、潮汐快速预报和余水位推算步骤。

2.如权利要求1所述的潮汐预报及水位推算方法，其特征在于：所述数值模拟为：利用σ坐标下的水动力数值模式构建目标海域潮波数值模式，模式加入引潮势、体潮、负载潮及自吸效应，从而更适用于深水海区。

ζ_g＝ζ+ζ_b+ζ_l.       (4)

其中，ζ_g为地心潮，ζ为海潮，ζ_b为体潮，ζ_l为海洋负载潮

体潮、海潮和负载潮均引起重力势的改变，作用在海水上的重力势Γ方程，

$\mathrm{\&Gamma;}=(1+k_2){\mathrm{\&Phi;}}_2+\underset{n}{\mathrm{\&Sigma;}}(1+{k^{\&prime;}}_n)U_n$

$=(1+k_2){\mathrm{\&Phi;}}_2+\underset{n}{\mathrm{\&Sigma;}}(1+{k^{\&prime;}}_n)g{\mathrm{\&alpha;}}_n{\mathrm{\&zeta;}}_n.---\left(12\right)$

考虑了引潮力的作用后，潮波模式控制方程中的水位ζ应改写为ζ′，

${\mathrm{\&zeta;}}^{\&prime;}={\mathrm{\&zeta;}}_g-\mathrm{\&Gamma;}/g$

$=\mathrm{\&zeta;}-(1+k_2-h_2){\mathrm{\&Phi;}}_2/g-\underset{n}{\mathrm{\&Sigma;}}(1+{k^{,}}_n-{h^{,}}_n){\mathrm{\&alpha;}}_n{\mathrm{\&zeta;}}_n---\left(13\right)$

$=\mathrm{\&zeta;}-\stackrel{\&OverBar;}{\mathrm{\&zeta;}}-{\mathrm{\&zeta;}}_{\mathrm{SAL}}$

其中，Φ_n为引潮势，h_n为Love数。由于主要的半日和全日分潮均属于n＝2的二阶球调和项，h_n′为负载Love数，U_n为ζ_n引起的重力势，

称为修正(即考虑体潮效应)后的平衡潮高，ζ_SAL称为自吸负载(self-attraction/load)项

3.如权利要求1或2所述的潮汐预报及水位推算方法，其特征在于：所述同化模式伴随方程为：

$\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;t}+g\frac{\&PartialD;\mathrm{\&mu;}}{\&PartialD;x}+g\frac{\&PartialD;\mathrm{\&upsi;}}{\&PartialD;y}+u\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;x}+v\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;y}+\frac{c_d\sqrt{u^2+v^2}}{{(h+\mathrm{\&zeta;})}^2}(\mathrm{\&mu;u}+\mathrm{\&upsi;v})=K_{\mathrm{\&zeta;}}(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_{\mathrm{obs}})---\left(33\right)$

$\frac{\&PartialD;\mathrm{\&mu;}}{\&PartialD;t}-\mathrm{\&mu;}\frac{\&PartialD;u}{\&PartialD;x}-\mathrm{\&upsi;}\frac{\&PartialD;v}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&mu;u}}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&mu;v}}{\&PartialD;y}-(f+\frac{c_d\mathrm{uv}}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}})\mathrm{\&upsi;}$

$+(h+\mathrm{\&zeta;})\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;x}-c_d\frac{{2u}^2+v^2}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}}\mathrm{\&mu;}+A(\frac{{\&PartialD;}^2\mathrm{\&mu;}}{\&PartialD;x^2}+\frac{{\&PartialD;}^2\mathrm{\&mu;}}{\&PartialD;y^2})=0---\left(34\right)$

$\frac{\&PartialD;\mathrm{\&upsi;}}{\&PartialD;t}-\mathrm{\&mu;}\frac{\&PartialD;u}{\&PartialD;y}-\mathrm{\&upsi;}\frac{\&PartialD;v}{\&PartialD;y}+\frac{\&PartialD;\mathrm{\&upsi;u}}{\&PartialD;x}+\frac{\&PartialD;\mathrm{\&upsi;v}}{\&PartialD;y}+(f-\frac{c_d\mathrm{uv}}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}})\mathrm{\&mu;}$

$+(h+\mathrm{\&zeta;})\frac{\&PartialD;\mathrm{\&lambda;}}{\&PartialD;y}-c_d\frac{u^2+2v^2}{(h+\mathrm{\&zeta;})\sqrt{u^2+v^2}}\mathrm{\&upsi;}+A(\frac{{\&PartialD;}^2\mathrm{\&upsi;}}{\&PartialD;x^2}+\frac{{\&PartialD;}^2\mathrm{\&upsi;}}{\&PartialD;y^2})=0---\left(35\right)$

其中，t代表时间；x，y是Cartesian坐标，分别取向东和向北为正；u(x，y，t)和v(x，y，t)分别代表x，y方向的流速分量；g是重力加速度；λ，μ，υ分别为ζ，u，v的伴随变量；h(x，y)表示(x，y)处的静水深度，ζ(x，y，t)表示未扰动海面以上的水位高度，H＝h+ζ为总水深；ζ_obs为观测值，K_ζ为权重系数，有观测值的点K_ζ取为1，否则取为0；

是Coriolis参数，Ω为地球自转角速度，

为地理纬度；c_d为底摩擦系数，A为水平涡动粘性系数。

4.如权利要求1或2所述的潮汐预报及水位推算方法，其特征在于：所述的潮汐快速预报为：相对于平均海面的潮高h可表示为

$h\left(t\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{f}_l\left(t\right)H_l\cos [{\mathrm{\&omega;}}_{l}t+v_{0,l}+u_l\left(t\right)-g_l]---\left(44\right)$

其中，H_l，g_l分别为第l个分潮的振幅和迟角，L是分潮个数；ω_l是分潮角速率；v_o，l为平衡潮第l个分潮在t＝0时刻的位相。这里f，ωt及v₀+u与时间有关，与地点无关；而H和g与地点有关，与时间无关。故上式可写作

$h\left(t\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}[{\mathrm{\&alpha;}}_l\left(t\right)X_l+{\mathrm{\&beta;}}_l\left(t\right)Y_l],---\left(45\right)$

其中

$\left\{\begin{array}{c}{\mathrm{\&alpha;}}_l\left(t\right)=f_l\left(t\right)\cos [{\mathrm{\&omega;}}_{l}t+v_{o,l}+u_l\left(t\right)],\\ {\mathrm{\&beta;}}_l\left(t\right)=f_l\left(t\right)\sin [{\mathrm{\&omega;}}_{l}t+v_{o,l}+u_l\left(t\right)];\end{array}\right.---\left(46\right)$

$\left\{\begin{array}{c}{X}_l=H_l\cos g_l,\\ Y_l=H_l\sin g_l.\end{array}\right.---\left(47\right)$

实际计算时不必对每个站点(层次)都计算α和β，可以一次性算出50年逐时的α，β值储存起来。对于每个站点(层次)，只要用该处的X，Y值代入式(45)即可。在本方法中，对潮位我们采用M₂，S₂，N₂，K₂，K₁，O₁，P₁，Q₁，M₄，MS₄，M₆，Sa和Ss_a 13个分潮，即L＝13。

Images