海洋占据了地球表面71%, 与我们人类的生存发展密切相关。潮汐作为海洋中最普遍的运动过程之一, 是海水在天体引潮力作用下形成的某些特定频率上的周期性波动现象, 在垂直方向上表现为潮位的周期性升降, 而在水平方向上则表现为潮流的周期性运动^[1]。潮汐的主要成因是地球表面上各点离月球和太阳的相对位置不同, 各点所受到的引潮力有所差异, 从而导致地球上的海水的相对运动^[2]。潮汐的原动力——引潮力, 可分解为一系列频率(周期)的分力, 而每一个频率对应一个简谐振动, 称之为分潮^[3]。对实际观测的水位资料进行经典潮汐调和分析可以求出各个主要分潮的振幅和迟角(称之为调和常数), 结合天文参数(交点因子和交点订正角)可以推算未来时刻的潮汐, 进行潮汐预报。潮汐与国民经济建设密切相关, 潮汐、潮流的准确预报在国防建设、航运交通、水产养殖、海洋资源开发、能源利用、环境保护、近海和远海工程建设以及海岸防护等各方面起着至关重要的作用^[1]。

潮汐可以直接影响海洋内波、海洋环流、风暴潮、海气相互作用、海洋生态系统等。潮致混合可以将次表层冷水带到海洋表面, 进而使海表面温度降低^[4]。一般认为潮致混合的强度与潮流流速的平方成正比^[5-6], 而潮流流速存在周期约半月的大小潮循环和18.61 a循环, 这导致潮汐混合强度也存在半月和18.61 a循环, 最终海表面温度也存在相应的周期^[7-11]。海面温度直接影响海表大气温度以及海洋生态系统^{[6, 9, 12]}。Royer^[13]指出在美国阿拉斯加锡特卡气温的30%的低频变化与18.61 a循环有关。Souza等^[14]指出南加利福利亚海域海表面温度和硅藻丰度变化的半月周期与潮致混合有关。Tanaka等^[15]发现在海气耦合模式里考虑潮汐混合的18.61 a循环可以提高年代际气候预报的准确性。

在层结的海洋中, 当正压潮流流经陆架、陆坡、海山、海洋中脊、岛弧等变化的海底地形时会产生与天文潮频率一致的海洋内波, 称为内潮^[16-18]。地形变化越剧烈的海域, 往往会产生越显著的海洋内潮。海洋内潮可以引起海洋内部等密度面上百米的起伏^[19-21], 在水平方向上可以传播1 000 km以上^[22-24], 同时也会在深海产生强烈的剪切与混合^{[19, 21, 25-26]}。海洋内潮作为从天文潮到混合过程的一个重要中间环节^[19], 是连接不同尺度海洋动力过程的重要桥梁^[27]。海洋内潮也是全球深海混合过程的主要能量来源。根据Munk等^[28]的估算, 维持全球海洋的层结需要约2.1 TW的能量, 而其中约一半的能量来自深海潮汐。海洋内潮及其引起的海水混合会影响局地环流流场, 进而影响海域内的物质输运和生态系统^{[27, 29-30]}。同时海洋内潮也是海洋内孤立波产生的一个重要因素^[31], 而且内潮还能与其他频率的海洋波动(比如近惯性内波)和中尺度涡等发生相互作用^{[27, 32-36]}。另外, 深海潮汐也是海洋物质输运的动力源泉, 特别是潮汐在海洋近底层可产生较强的剪切应力, 导致海底沉积物再悬浮和物质输运^[37-41]。总之, 潮汐学的研究不仅具有重要的实际应用价值, 还具有重要的学术价值^{[1, 42-44]}。

1 潮汐数据分析方法

经典潮汐调和分析是使用最广泛的潮汐数据分析方法。在经典潮汐调和分析方法中, 各个分潮的振幅和迟角都是常数, 这种假设对于平稳潮是合理的。但是对于河流潮汐、内潮等非平稳潮, 这种假设是不合理的。河流潮汐是天文潮与地形、底摩擦和河流径流非线性相互作用的结果。潮波在沿河传播的过程中会逐渐被扭曲和衰减^[45-51]。对于这些高度非平稳潮汐时间序列, 经典潮汐调和分析方法只能得到各主要分潮的振幅和迟角的时间平均值。为了能够更好地理解潮与非潮过程的相互作用, 国内外学者开展了大量的相关研究。Jay等^[52]首次将小波分析引入河流潮汐水位资料的分析, 小波分析相对于傅里叶分析和经典调和分析的优势在于它能得到某些频率振幅的时间变化, 因此被广泛用于包括河流潮汐在内的非平稳潮的分析^[52-57]。然而, 小波分析虽然能反映信号时域和频域振幅的变化, 但是却丢失了潮频段分潮的分辨能力^[58]。

Pan等^[59]首次将经验模态分解(empirical mode decomposition, 简称EMD)引入河流潮汐的分析。经验模态分解^[60]是依据数据自身的时间尺度特征来进行信号分解, 无须预先设定基函数。这一点与建立在先验性的谐波基函数(小波基函数)上的傅里叶分析(小波分析)方法具有实质性的差别。正是由于这样的特点, EMD方法在理论上可以应用于任何类型的信号的分解, 因而在处理非平稳及非线性数据上, 具有很高的信噪比, 有明显的优势。因此EMD方法一经提出就在不同的工程领域得到了广泛的应用, 例如在海洋、大气、天体观测资料与地震记录分析、机械故障诊断、密频动力系统的阻尼识别以及大型土木工程结构的模态参数识别方面^{[59, 61-65]}。EMD方法能使复杂信号分解为有限个本征模函数(intrinsic mode function, 简称IMF), 所分解出来的各IMF分量包含了原信号的不同时间尺度的局部特征信号。Flandrin等^[66]研究表明EMD是一个二分滤波器(每个EMD模态的周期都近似是前面的2倍)。对于河流潮汐数据进行EMD分解得到的第一模态代表半日潮, 第二模态代表全日潮, 剩余模态的叠加是平均海平面(mean water level, MWL)。EMD成功捕捉到了全日分潮, 半日分潮以及平均海平面的由径流引起的非平稳变化特征。和小波分析一样, EMD同样丢失了潮频段分潮的分辨能力^[58-59]。张亮等^[67]使用了改进的集合经验模态分解(MEEMD)对海洋潮汐进行分析, 该方法解决了EMD的模态混淆问题, 但是仍然不具备潮频段分潮的分辨能力。

Matte等^[58]首次尝试改进经典潮汐调和分析模型, 使之能应用于河流潮汐。Matte等^[58]基于Kukulka等^[68-69]的潮汐-径流相互作用的理论框架, 将非平稳的外力强迫直接引入到调和分析基础函数中, 开发了专门用于河流潮汐分析的工具包NS_TIDE。Matte等^[58]将NS_TIDE用于分析美国哥伦比亚河的Vancouver验潮站的水位数据, Gan等^[70]和Yu等^[71]将NS_TIDE应用于中国长江口水位数据分析, 他们的结果均表明NS_TIDE在河流潮汐数据分析中的表现远远优于经典调和分析工具包T_TIDE^[72]。

Jin等^[73]提出了非平稳潮汐调和分析方法(enhanced harmonic analysis, EHA), 并首次应用于内潮流速资料的分析, 理念上将经典调和分析理论中的振幅和迟角由常数改进为随时间变化的函数, 将独立点方案^{[44, 74-76]}和三次样条插值方法引入调和分析方法中, 利用最小二乘方法对随时间变化的振幅和迟角进行求解。Pan等^[77]基于广泛使用的T_TIDE工具包为EHA开发了S_TIDE工具包, 并首次使用S_TIDE研究了美国哥伦比亚河口的潮汐-径流相互作用。相比NS_TIDE只能应用于河流潮汐, S_TIDE理论上可以应用于所有类型的潮汐。

NS_TIDE和S_TIDE作为潮汐调和分析模型的最新进展, 引起了大家的兴趣和关注。本文首先回顾经典潮汐调和分析模型的原理, 然后重点介绍NS_ TIDE和S_TIDE模型的原理, 指出它们的优缺点, 并探讨未来潮汐调和分析模型的发展方向。

2 经典潮汐调和分析模型

在经典潮汐调和分析模型中[如式(1)所示], 水位被认为是一系列三角函数线性叠加的结果^[78-81], 每一个三角函数都代表了一个分潮:

$ Z(t)=S_{0}+\sum\limits_{j=1}^{J}\left[H_{j} \cos \left(\sigma_{j} t-g_{j}\right)\right], $

(1)

Z(t)是t时刻观测到的水位, σ_j、H_j以及g_j是第j个分潮对应的角速度、振幅和迟角。J是分辨的分潮的个数。S₀是平均海平面。式(1)可以使用变量a_j和b_j线性化表达并改写为式(2):

$ Z(t)=S_{0}+\sum\limits_{j=1}^{J}\left(a_{j} \cos \sigma_{j} t+b_{j} \sin \sigma_{j} t\right), $

(2)

式(1)中

$ H_{j}=\sqrt{a_{j}^{2}+b_{j}^{2}}, \quad g_{j}=\arctan \left(b_{j} / a_{j}\right). $

(3)

在做线性回归时, 根据式(2)用矩阵的形式表示, 即式(4):

$ \boldsymbol{Z}=\boldsymbol{A} \boldsymbol{X}, $

(4)

其中, Z是观测水位矩阵, A是已知的系数矩阵, X是待求解的未知参数矩阵, 具体表达如式(5)和(6)所示:

$ \boldsymbol{Z}=\left(\begin{array}{c} Z\left(t_{1}\right) \\ Z\left(t_{2}\right) \\ \cdots \\ Z\left(t_{N}\right) \end{array}\right), \quad \boldsymbol{X}=\left(\begin{array}{c} S_{0} \\ a_{1} \\ \cdots \\ a_{J} \\ b_{1} \\ \cdots \\ b_{J} \end{array}\right), $

(5)

$ \boldsymbol{A}=\left(\begin{array}{ccccccc} 1 & \cos \sigma_{1} t_{1} & \cdots & \cos \sigma_{J} t_{1} & \sin \sigma_{1} t_{1} & \cdots & \sin \sigma_{J} t_{1} \\ 1 & \cos \sigma_{1} t_{2} & \cdots & \cos \sigma_{J} t_{2} & \sin \sigma_{1} t_{2} & \cdots & \sin \sigma_{J} t_{2} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 4 & \cos \sigma_{1} t_{N} & \cdots & \cos \sigma_{J} t_{N} & \sin \sigma_{1} t_{N} & \cdots & \sin \sigma_{J} t_{N} \end{array}\right). $

(6)

根据最小二乘法, X的解为:

$ {\mathit{\boldsymbol{X}}_{\bf{s}}} = {\left( {{\mathit{\boldsymbol{A}}^{\bf{T}}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\bf{T}}}\mathit{\boldsymbol{Z}}. $

(7)

为了减小潮汐调和分析所求得的分潮振幅和迟角的置信区间, Godin^[82]提出了瑞利准则来限制调和分析里使用的分潮的个数[式(8)]:

$ \left|\sigma_{1}-\sigma_{2}\right|>\frac{1}{m \Delta t} $

(8)

其中, $ {\sigma _1} $和$ {\sigma _2} $是两个相邻分潮的频率, m是采样点的个数, $ \Delta t $是采样间隔, $ m\Delta t $是数据长度(length of record, LOR)。根据瑞利准则, 调和分析里能分辨的分潮的最小频率差只与数据长度有关。数据越长, 能分辨的分潮越多。

对美国哥伦比亚河5个验潮站观测的水位数据应用经典潮汐调和分析(借助于T_TIDE工具包), 结果如表 1所示。Astoria是离海洋最近的站点, 河流径流影响较小, 经典潮汐调和分析回报水位解释了原信号96.7%的变化, 均方根误差0.15 m, 最大绝对误差1.18 m。潮波在沿河道向内陆传播的过程中, 受河流径流的影响越来越大, 而经典潮汐调和分析并没有考虑潮汐-径流的非线性相互作用, 导致离海洋越远的站点, T_TIDE的表现就越差(表 1)。Vancouver是离海洋最远的站, 也是T_TIDE表现最差的, 其回报的水位只能解释原信号的52.4%, 均方根误差是1.60 m, 最大绝对误差是4.16 m。

3 非平稳潮汐调和分析模型NS_TIDE 3.1 NS_TIDE模型介绍

经典潮汐调和分析模型认为每一个分潮都是完美的三角函数, 即分潮的振幅和迟角都是不随时间变化的常数。在感潮河段, 由于潮汐-径流非线性相互作用, 导致主要分潮的振幅和迟角不再是常数, 而是河流径流的非线性函数。由于不合理的平稳性假设, 经典潮汐调和分析方法在河口地区表现不尽人意。基于Kukulka等^[68-69]的潮汐-径流相互作用的理论模型, Matte等^{[58, 83]}在T_TIDE的基础上开发了专门用于河流潮汐数据分析的工具包NS_TIDE。NS_TIDE认为分潮的振幅和迟角都是随时间变化的[式(9)], 并将非平稳的外力强迫(河流径流和外海潮汐)直接引入到调和分析基础函数中[式(10)]。

$ Z(t) = {S_{0, 0}}(t) + \sum\limits_{j = 1}^J {\left( {{S_{1, j}}(t)\cos {\sigma _j}t + {S_{2, j}}(t)\sin {\sigma _j}t} \right)} , $

(9)

$ {S_{l, j}}(t) = {d_{0, l, j}} + {d_{1, l, j}}Q{(t)^p} + {d_{2, l, j}}\frac{{R{{(t)}^q}}}{{Q{{(t)}^r}}} , $

(10)

其中, Q代表河流径流, R代表外海参考站的全日潮差。J代表分析的分潮总个数。l是模型系数(l=0代表亚潮模型系数, l=1, 2代表潮汐-径流模型系数)。p、q以及r是待定的常数, 其理论值由Kukulka等^[68-69]给出。如果知道Z、Q、R, 通过最小二乘的方法就能求出未知的模型系数d, 进而量化非平稳外力对于潮汐振幅和迟角的影响。在NS_TIDE中, 为了减少离群值对于调和分析结果的影响, NS_TIDE使用了迭代权重最小二乘法(the iteratively reweighted least squares, IRLS)^[84-85]。

对于河口潮等非平稳潮, 瑞利准则[式(8)]是不适用的。Munk等^[86]指出潮汐谱线周围形成的尖点宽度与非线性强迫的低频谱有关。从这一事实出发, NS_TIDE定义了适用于非平稳潮的瑞利准则, 它主要考虑到主要分潮附近的小分潮对于主要分潮能谱的扭曲:

$ \left| {{\sigma _1} - {\sigma _2}} \right| > \max \left( {\frac{1}{{m\Delta t}}, \Delta \sigma } \right) , $

(11)

$ \Delta \sigma $由下式给出:

$ \frac{{\int_0^{\Delta \sigma } {P(\sigma )d\sigma } }}{{\int_0^\infty {P(\sigma )d\sigma } }} = 1 - \eta , $

(12)

其中, $ P(\sigma ) $是$ {Q^p} $或者是$ 1/{Q^r} $的标准化能量谱, $ \eta $是用户自己定义的参数, 代表占全部谱能量的比例。$ \eta $越小, $ \Delta \sigma $就越大, 能分辨的分潮就越少, 计算得到的分潮振幅里来自邻近的未分辨的分潮的贡献就越大, 反之亦然^[58]。

在哥伦比亚河上游, NS_TIDE的表现要明显优于T_TIDE。以离海洋最远的Vancouver验潮站为例: NS_ TIDE的回报水位解释了观测水位的96.0%的方差变化, 均方根误差为0.155 m, 最大绝对误差为1.018 m。NS_ TIDE得到的Vancouver处的M₂和K₁的时间变化的振幅如图 1所示。根据T_TIDE计算结果, 在Vancouver, M₂分潮平均振幅为0.21 m, 平均迟角为247.66°; 而K₁分潮平均振幅为0.12 m, 平均迟角为15.85°。可以看到NS_TIDE得到的振幅和迟角的变化都在T_TIDE得到的调和常数附近波动, 随着河流径流的变化而变化。当河流径流显著增加时, 潮波传播受到的阻碍也显著增加, 潮波能量被消耗, 分潮振幅会显著减小, 迟角会有快速变化, 反之亦然。K₁和M₂振幅和迟角的变化也存在由于大小潮产生的准半月波动。

3.2 NS_TIDE模型优缺点

NS_TIDE不仅保持了潮频段分潮的分辨能力, 而且还可以量化非平稳外力(河流径流)对于潮汐非平稳变化的贡献。因此, NS_TID在美国哥伦比亚河^{[58-59, 77]}、加拿大圣劳仑斯河^{[58, 87-88]}、中国长江^{[70-71, 89]}和珠江^[90]等全球各大河口都得到了广泛应用。但是NS_TIDE也存在一些不足^{[59, 77]}。第一, 河流高频径流数据的缺乏会限制NS_TIDE的使用。第二, NS_TIDE基于潮汐-径流相互作用模型, 这也就导致了NS_TIDE不能用于其他的非平稳潮(比如内潮)。第三, NS_TIDE需要河流参考站的逐时水位数据。而且, NS_TIDE只考虑了潮汐-径流的相互作用, 忽略了其他的重要物理过程比如沿岸上升流^[91]对于水位变化的影响。第四, NS_TIDE的待定指数需要进行迭代计算和优化, 导致计算效率降低。

4 非平稳潮汐调和分析模型S_TIDE 4.1 S_TIDE模型介绍

与NS_TIDE类似, 如式(13)所示, S_TIDE同样认为分潮的振幅和迟角是随时间变化的:

$ Z(t) = S(t) + \sum\limits_{j = 1}^J {\left( {{a_j}(t)\cos {\sigma _j}t + {b_j}(t)\sin {\sigma _j}t} \right)} . $

(13)

但是不同于NS_TIDE, 在S_TIDE中, 式(13)的求解使用了独立点方案^{[44, 73, 92-95]}。独立点方案的基本思想是基于样条函数方法。如图 2所示, 假定调和变量(比如振幅)是随时间变化的(图 2中黑色实线), X={x₁, x₂, x₃, …, x₉} 是不同的时刻。首先, 选择$ \widetilde X $= {x₁, x₃, x₅, x₇, x₉}作为独立点(independent points, IPs)代表整个变量空间。然后, 独立点处的调和变量的值(图 2中的红色圆圈)通过特定的算法可以求解(细节稍后给出)。最后, 调和变量在非独立点处(x₂, x₄, x₆, x₈)的值可以通过独立点处的插值得到。

S_TIDE的求解过程可以表示为如下步骤。对于平均海平面S和第j个分潮潮汐调和变量(a_j和b_j), 在时间域上均匀选取一系列点作为独立点(标记为S_i, a_{i, j}和b_{i, j}), 而非独立点处的值可以通过独立点处的插值得到。因此, 随时间变化的MWL和调和变量可以用独立点来表达S_i、a_{i, j}和b_{i, j}:

$ S(t) = \sum\limits_{i = 1}^{{M_s}} {{f_{t, i}}{S_i}} , \;{a_j}(t) = \sum\limits_{i = 1}^M {{f_{t, i}}{a_{i, j}}} , \;{b_j}(t) = \sum\limits_{i = 1}^M {{f_{t, i}}{b_{i, j}}} , $

(14)

其中, f_{t, i}是第i个独立点在t时刻的插值权重, 只与选用的插值方法有关。M_s和M分别是MWL和分潮使用的独立点个数。要注意的是, 目前S_i、a_{i, j}和b_{i, j}还是待定的, 但是基于独立点方案, 随时间变化的S(t)、a_j(t)和b_j(t)被表达成了独立点(S_i, a_{i, j}, b_{i, j})处值的线性叠加[式(14)]。将式(13)和(14)组合得到式(15):

$ Z(t) = \sum\limits_{i = 1}^{{M_s}} {{f_{t, i}}} {S_i} + \sum\limits_{j = 1}^J {\left( {\sum\limits_{i = 1}^M {{f_{t, i}}{a_{i, j}}} \cos {\sigma _j}t + \sum\limits_{i = 1}^M {{f_{t, i}}{b_{i, j}}} \sin {\sigma _j}t} \right)} . $

(15)

三次样条插值由于其稳定、收敛和光滑的性质在数据处理中得到了广泛使用^{[44, 60, 96-97]}。因此, S_TIDE选择了三次样条插值连接独立点和非独立点。关于三次样条插值权重的具体计算公式可以参考Pan等^[77], 这里不再赘述。

如果有N个时刻的水位观测值, 即Z = Z₁, Z₂, …, Z_N, 在t = t₁, t₂, …, t_N, 时刻, 可以得到下面的方程组:

$ \left\{ \begin{array}{l} \sum\limits_{i = 1}^{{M_s}} {{f_{{t_1}, i}}} {S_i} + \sum\limits_{j = 1}^J {\left( {\sum\limits_{i = 1}^M {{f_{{t_1}, i}}{a_{i, j}}} \cos {\sigma _j}{t_1} + \sum\limits_{i = 1}^M {{f_{{t_1}, i}}{b_{i, j}}} \sin {\sigma _j}{t_1}} \right)} = Z({t_1}) \hfill \\ \sum\limits_{i = 1}^{{M_s}} {{f_{{t_2}, i}}} {S_i} + \sum\limits_{j = 1}^J {\left( {\sum\limits_{i = 1}^M {{f_{{t_2}, i}}{a_{i, j}}} \cos {\sigma _j}{t_2} + \sum\limits_{i = 1}^M {{f_{{t_2}, i}}{b_{i, j}}} \sin {\sigma _j}{t_2}} \right)} = Z({t_2}) \hfill \\ \cdots \cdots \hfill \\ \sum\limits_{i = 1}^{{M_s}} {{f_{{t_N}, i}}} {S_i} + \sum\limits_{j = 1}^J {\left( {\sum\limits_{i = 1}^M {{f_{{t_N}, i}}{a_{i, j}}} \cos {\sigma _j}{t_N} + \sum\limits_{i = 1}^M {{f_{{t_N}, i}}{b_{i, j}}} \sin {\sigma _j}{t_N}} \right)} = Z({t_N}). \hfill \\ \end{array} \right. $

(16)

方程组(16)里有(2M·J+M_s)个未知参数, 当观测数N远远大于(2M·J+M_s)时可以通过最小二乘法来求解这些未知参数(S_i, a_{i, j}和b_{i, j})。一旦独立点处的值确定, 就可以三次样条插值得到时间变化的MWL、分潮振幅和迟角。

在做线性回归时, 与经典潮汐调和分析类似, 式(16)也可以写成矩阵的形式, 即式(17):

$ \boldsymbol{Z}=\boldsymbol{K} \boldsymbol{U} \text {, } $

(17)

其中Z是观测水位矩阵, K是已知的系数矩阵, U是待求解的未知参数矩阵, 具体表达如式(18)、(19)和(20)所示:

$ Z = {\left[ {Z\left( {{t_1}} \right), Z\left( {{t_2}} \right), \cdots , Z\left( {{t_N}} \right)} \right]^T} , $

(18)

$ U = {\left( {{S_1}\;\;\; \cdots \;\;\;{S_{{M_s}}}\;\;\;{a_{1, 1}}\;\;\; \cdots \;\;\;{a_{M, 1}}\;\;\; \cdots \;\;\;{a_{1, J}}\;\;\; \cdots \;\;\;{a_{M, J}}\;\;\;{b_{1, 1}}\;\;\; \cdots \;\;\;{b_{M, 1}}\;\;\; \cdots \;\;\;{b_{1, J}}\;\;\; \cdots \;\;\;{b_{M, J}}} \right)^T}, $

(19)

$ K = \left( {\begin{array}{*{20}{c}} {{f_{{t_1}, 1}}}& \cdots &{{f_{{t_1}, {M_s}}}}&{{f_{{t_1}, 1}}\cos {\sigma _1}{t_1}}& \cdots & \cdots &{{f_{{t_1}, M}}\sin {\sigma _J}{t_1}} \\ {{f_{{t_2}, 1}}}& \cdots &{{f_{{t_2}, {M_s}}}}&{{f_{{t_2}, 1}}\cos {\sigma _1}{t_2}}& \cdots & \cdots &{{f_{{t_2}, M}}\sin {\sigma _J}{t_2}} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {{f_{{t_N}, 1}}}& \cdots &{{f_{{t_N}, {M_s}}}}&{{f_{{t_N}, 1}}\cos {\sigma _1}{t_N}}& \cdots & \cdots &{{f_{{t_N}, M}}\sin {\sigma _J}{t_N}} \end{array}} \right) . $

(20)

根据最小二乘法, U的最优解为:

${\mathit{\boldsymbol{U}}_{\bf{s}}} = {\left( {{\mathit{\boldsymbol{K}}^{\bf{T}}}\mathit{\boldsymbol{K}}} \right)^{ - 1}}{\mathit{\boldsymbol{K}}^{\bf{T}}}\mathit{\boldsymbol{Z}}. $

(21)

独立点个数和使用的分潮是S_TIDE最核心的两个输入参数, 如果这两个参数取值不合适, S_TIDE输出结果可能是没有物理意义的。根据目前的研究^{[77, 98]}, 独立点个数越多, S_TIDE得到的分潮振幅和迟角的变化就越复杂; 独立点越少, 得到的振幅和迟角的变化就越简单。当独立点个数为2时, S_TIDE得到的分潮振幅和迟角的变化是线性的; 当独立点个数为1时, S_TIDE得到的分潮振幅和迟角都是不随时间变化的, 与经典调和分析结果是完全一致的。从某种意义上来说, 经典潮汐调和分析模型可以看成是S_TIDE在独立点为1时的一个特例。独立点个数和分潮的个数是相互制约的。使用的独立点越多, S_TIDE得到的分潮振幅和迟角的变化就越复杂, 能分辨的分潮就越少, 计算得到的分潮振幅里来自邻近的未分辨的分潮的贡献就越大, 反之亦然。

以下以M₂分潮振幅的年循环为例, 介绍为什么在S_TIDE模型里独立点个数和分潮的个数是相互制约的。如式(22)所示, σ_a、h_a以及g_a是M₂分潮振幅年循环的频率、振幅和迟角, σ_M2、h_M2以及g_M2是M₂分潮振幅的频率、振幅和迟角, t是时间。对式(22)做积化和差展开, 得到式(23)。虽然式(22)与(23)在数学上是完全等价的, 但是它们却代表了对同一信号的两种完全不同的解读方式。式(22)认为信号里只存在一个频率, 那就是M₂分潮的频率, 但是M₂分潮振幅是存在年变化的。式(23)中出现3个分潮的频率, 分别是M₂、H₁和H₂这3个分潮的频率, 而且这3个分潮振幅都是不随时间变化的常数。式(22)代表了S_TIDE理解信号的方式, 即分潮的振幅和迟角是时间变化的; 而式(3)—(10)则代表了T_TIDE理解信号的方式, 即分潮的振幅和迟角都是常数。如果要用S_TIDE提取M₂分潮振幅的年变化, 那么在S_TIDE模型里, H₁和H₂分潮不能区分, 因为H₁和H₂分潮是由于M₂分潮的年循环产生的。

$ \left[ {{h_{M2}} + {h_a}\cos \left( {{\sigma _a}t + {g_a}} \right)} \right]\cos \left( {{\sigma _{M2}}t + {g_{M2}}} \right) , $

(22)

$ {h_{M2}}\left( {\cos {\sigma _{M2}}t + {g_{M2}}} \right) + 0.5{h_a}\cos \left[ {\left( {{\sigma _{M2}} + {\sigma _a}} \right)t + {g_{M2}} + {g_a}} \right] + 0.5{h_a}\cos \left[ {\left( {{\sigma _{M2}} - {\sigma _a}} \right)t + {g_{M2}} - {g_a}} \right] . $

(23)

S_TIDE在哥伦比亚河的表现要略优于NS_ TIDE。以Vancouver验潮站为例: S_TIDE回报的水位解释了原始观测99.1%的方差变化, 均方根误差为0.075 m, 最大绝对误差为0.48 m。图 3展示了S_TIDE回报水位, 几乎与观测水位一致, 回报的残差非常小。相比之下, T_TIDE的回报误差非常大。

4.2 S_TIDE模型优点

1) 相比小波分析和NS_TIDE, S_TIDE更加简单易懂, 更容易实现。由S_TIDE得到的结果与经验模态分解, 小波分析还有NS_TIDE高度一致, 也符合河口潮解析理论。但是相比NS_TIDE, S_TIDE回报更加准确。

2) 与经典潮汐调和分析类似, S_TIDE只假定预先知道潮汐频率, 无需知道其他的动力过程。也就是说S_TIDE不像NS_TIDE那样依赖于河流径流数据或者其他外界强迫数据, 这使得S_TIDE可以广泛应用于其他非平稳潮汐过程(比如内潮)。不仅如此, 尽管S_TIDE是为非平稳潮设计的, 但是其实它也可以用于平稳潮, 可以用来分析振幅和水位的长期变化和季节变化等。

3) 可以使用S_TIDE来准确分离水位里的潮汐和亚潮波动。也可以通过调整独立点的个数来提取MWL、分潮振幅和迟角在不同时间尺度上的波动。

4.3 S_TIDE模型的不足

1) S_TIDE是处理非平稳非线性数据的有力工具。NS_TIDE可以用于预报; 而S_TIDE目前还未实现, 尚需深化研究。

2) 当独立点个数较大时, 相比其他的方法, S_ TIDE需要更多的计算资源和时间。另外, 太多的独立点可能会导致计算机内存的溢出和不符合实际的结果。例如, 对于哥伦比亚河Vancouver验潮站7年的水位数据, 当分潮的独立点超过了500个时, 反演得到的振幅的大小潮变化由于过度拟合被明显夸大了。因此, 在使用S_TIDE时要非常谨慎选择独立点个数和使用的分潮。

5 结论

经典潮汐调和分析作为使用最广泛的潮汐数据分析方法, 从1921年至今, 已历经100年的发展历程。Doodson^[78-80]发展并完善了经典的潮汐调和分析理论。方国洪于1960年发展了Doodson短期观测的分析方法, 对引潮势进行了准调和展开, 提出了潮汐分析和预报的准调和分析方法^[99-101], 为我国的潮汐、潮流分析预报工作做出了重要贡献。当代我国的潮汐潮流研究, 陈宗镛^[102]的《潮汐学》, 方国洪等^[1]的《潮汐和潮流的分析和预报》, 黄祖珂等^[103]的《潮汐原理与计算》等专著发挥了重要作用。

经典调和分析方法理论上没有考虑季节变化, 而主要分潮的季节变化是客观存在的, 亟待准确的刻画和提取。尽管在经典调和分析的框架下, 能够获取主要分潮的季节变化, 但在具体求解过程中必须假定某个时间段内各分潮的调和常数是不变的, 而这个时间段长度如何确定是个无法解决的问题: 这个时间段长度过短则无法准确分离不同分潮, 过长则无法准确提取出各个分潮随时间变化; 需要特别指出的是上述的求解过程不可避免地会弱化振幅(迟角)的变化幅度。根据调和分析理论和方法的发展趋势和现实需求, 潮汐调和分析模型研究未来可能的发展方向之一是改进经典的调和分析方法, 把经典调和分析方法中主要分潮的调和常数(振幅和迟角)推广位变化的振幅(迟角), 开展变化振幅(迟角)获取的理论和方法研究。

致谢: 谨以此文致敬在潮汐调和分析模型发展及应用历程中有所建树的先贤和前辈们, 特别感谢引领我入门物理海洋学研究的授业恩师方国洪先生。本文使用的S_TIDE工具包可以从如下网址获取: https://www.researchgate.net/project/A-non-stationary-tidal-analysis-toolbox-S-TIDE。