南海是一个半封闭海, 其与外界水体交换的主要通道是吕宋海峡^[1], 南海北部具有宽广的大陆架和大陆坡, 陆架区潮汐现象显著^[2]。因此, 针对南海潮汐的研究一直在国内外颇受重视。TOPEX/ POSEIDON(T/P)卫星是由美国国家航空航天局和法国空间局联合于1992年8月10日发射的世界上第一颗专门为研究世界大洋环流而设计的高度计卫星。T/P卫星在2002年9月变轨, 原轨道的任务由2001年12月发射的Jason-1卫星执行。Jason-1卫星于2009年2月变轨到T/P卫星变轨后的轨道上继续运行。原有轨道测量任务随后交给2008年6月发射的Jason-2卫星。T/P、Jason-1和Jason-2三颗卫星高度计系统的定轨精度和测高精度较以前有显著提高, 其测量精度约为5 cm, 为潮汐的研究提供了重要的数据源。鉴于T/P和Jason卫星高度计观测海面高度精度高的特性, 国内外很多学者将其应用到潮汐研究中。Yanagi等^[3]基于3年的卫星高度计资料及沿岸验潮站的调和分析结果, 以指数函数为插值系数插值得到了中国东海和黄海八个主要分潮的调和常数。Fang等^[4]利用卫星高度计资料及验潮站数据进行二次多项式拟合, 得到了渤、黄、东海网格点的调和常数, 进而得到了渤、黄、东海五个主要分潮的调和常数。范丽丽^[5]基于物理量在空间上的连续性, 提出一种通过轨道上的调和常数拟合轨道间任意点调和常数的方法, 并将该方法应用于夏威夷海域。

本文数据采用1992~2016年T/P和Jason卫星高度计海面高度资料, 对南海海域(2°~25°N, 99°~ 122°E, 图 1)的潮汐信号进行了提取。然后, 采用以格林函数为基础的双调和样条插值法对所提取的M₂、S₂、K₁和O₁四个主要分潮调和常数进行了空间二维插值, 得到了南海四个主要分潮的同潮图, 并利用全球潮汐模型TPXO7.2、验潮站数据和前人研究中无潮点位置对插值结果进行了评估。

1 数据和方法 1.1 卫星高度计数据

变轨前的T/P卫星轨道高度达1 336 km, 倾角为66°, 覆盖面大, 轨道的交点周期(绕地球一圈的时间)为6 745.8 s, 轨道运行127圈以后精确重复, 重复周期为9.915 6 d。本文高度计数据来自南海海域(2°~ 25°N, 99°~122°E)变轨前T/P卫星轨道上1992~2016年多卫星数据, 包含: 1992~2002年未变轨的T/P卫星高度计数据、2002~2008年未变轨的Jason-1卫星高度计数据和2008~2016年Jason-2卫星高度计数据(下载地址: ftp://avisoftp.cnes.fr/AVISO/pub/), 共有4 300多个星下观测点, 每个星下观测点都对应一个海面高度的时间序列。

1.2 验潮站数据

本文中的141个验潮站(验潮站分布情况见图 1)数据来自于中国近海及周边海域验潮站的资料和国际上公开的验潮站资料, 包括M₂、S₂、K₁和O₁四个主要分潮的潮汐调和常数。

1.3 调和分析方法

实际水位可以看作是很多个分潮所致水位迭加的结果, 但是在实际分析中只能选取有限个主要分潮。假设我们选取J个分潮, 任意一点的潮位h的表达式为:

$ h = {S_{\rm{0}}} + \sum\limits_{j = 1}^J {{f_j}{H_j}\cos ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j} - {g_j})}, $

(1)

其中, S₀为余水位, f_j为第j个分潮的交点因子, u_j为第j个分潮的交点订正角, H_j、g_j分别为第j个分潮的振幅和迟角, 统称为调和常数, σ_j为第j个分潮的分潮角速度, V_0j为第j个分潮的初相位, t为时间。

某一分潮所引起的水位变化为:

$ {f_j}{H_j}\cos ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j} - {g_j}). $

(2)

式(2)可化为式(3):

$ \begin{array}{l} {f_j}\cos ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j}){H_j}\cos {g_j} + \\ {f_j}\sin ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j}){H_j}\sin {g_j} \end{array}, $

(3)

其中, 时间t和分潮角速度σ_j均为已知, 交点因子交点f_j、交点订正角u_j和分潮初相位V_0j可根据Doodson数和时间t计算得出^[6]。对应J个分潮, 则有:

$ \begin{array}{l} h = {S_{\rm{0}}} + \sum\limits_{j = 1}^J {{f_j}\cos ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j})} {H_j}\cos {g_j} + \\ \sum\limits_{j = 1}^J {{f_j}\sin ({\sigma _j}t + {V_{{\rm{0}}j}} + {u_j})} {H_j}\sin {g_j} \end{array}. $

(4)

如果在n个时刻$t = {t_{\rm{1}}}, \;{t_{\rm{2}}}{\rm{, }}\; \cdots , \;{t_n}$, 有n个潮位观测值$h = {h_{\rm{1}}}, \;{h_{\rm{2}}}, \; \cdots , \;{h_n}$, 那么, 就可以建立如下由n个方程构成的方程组:

$ \left\{ {\begin{array}{*{20}{c}} {{S_{\rm{0}}} + \sum\limits_{j = 1}^J {{x_j}{f_j}\cos ({\sigma _j}{t_1} + {V_{\rm{0}}} + {u_j}) + \sum\limits_{j = 1}^J {{y_j}{f_j}\sin ({\sigma _j}{t_1} + {V_{\rm{0}}} + {u_j}) = {h_{\rm{1}}}} } }\\ {{S_{\rm{0}}} + \sum\limits_{j = 1}^J {{x_j}{f_j}\cos ({\sigma _j}{t_2} + {V_{\rm{0}}} + {u_j}) + \sum\limits_{j = 1}^J {{y_j}{f_j}\sin ({\sigma _j}{t_2} + {V_{\rm{0}}} + {u_j}) = {h_{\rm{2}}}} } }\\ \vdots \\ {{S_{\rm{0}}} + \sum\limits_{j = 1}^J {{x_j}{f_j}\cos ({\sigma _j}{t_n} + {V_{\rm{0}}} + {u_j}) + \sum\limits_{j = 1}^J {{y_j}{f_j}\sin ({\sigma _j}{t_n} + {V_{\rm{0}}} + {u_j}) = {h_n}} } } \end{array}} \right.. $

(5)

方程(5)中x_j对应${H_j}\cos {g_j}$, y_j对应${H_j}\sin {g_j}$, 它们和S₀共同构成了方程组中的全部未知量。潮位观测个数n往往大于未知量个数2J+1, 因此需要利用最小二乘法对上述方程组进行求解。求得x_j和y_j后, 即可得到各个分潮的调和常数H和g。本文中J取4, 即只考虑M₂、S₂、K₁和O₁四个主要分潮。

1.4 双调和样条插值

双调和样条插值法是以格林函数为基础的内插法, 通过插值计算, 找到通过空间数据点的曲率最小的曲线或曲面。在二维空间中, 满足平面上N个点的双调和方程为:

$ \begin{array}{*{20}{c}} {{\nabla ^4}\zeta (x) = \sum\limits_{j = 1}^N {{\alpha _j}\delta (x - {x_j})} }\\ {\zeta ({x_j}) = {\zeta _j}, } \end{array} $

(6)

其中, ${\nabla ^4}$为双调和算子, x为二维平面中的一个坐标点, ζ(x)为内插结果在x处的位移。α_j为在x_j处的加权系数, δ(x–x_j)为狄拉克函数, 即

$ \delta (x - {x_j}) = \left\{ \begin{array}{l} 0, x \ne {x_j}\\ \infty , x = {x_j} \end{array} \right.且\int_{\; - \infty }^{\; + \infty } {\delta (x){\rm{d}}x} = 1.$

(7)

双调和方程的通解为

$ \zeta (x) = \sum\limits_{j = 1}^N {{\alpha _j}\Phi (x - {x_j})} . $

(8)

加权系数α_j可以通过求解线性方程${\zeta _j} = \sum\limits_{i = 1}^N {{\alpha _j}\mathit{\Phi }({x_i} - {x_j})} $得到, 其中$\mathit{\Phi }(x) = |x{|^2}({\rm{ln}}|x| - 1)$为格林函数。

双调和样条插值综合利用数据值和斜率确定插值曲面, 比样条插值更为灵活, 此外, 即使是带噪声的数据, 也能在最小二乘意义下实现有效插值, 特别适合对不规则的卫星高度计数据进行插值^[7-8]。

1.5 评估方法

本文采用振幅差ΔH和迟角差Δg分别评估振幅和迟角的误差, 采用矢量差Δ综合评估振幅和迟角的误差^[4]。矢量差Δ表达式如下:

$ \Delta {\rm{ = [(}}{H_{\rm{a}}}\cos {g_{\rm{a}}} - {H_{\rm{b}}}{\rm{cos}}{g_{\rm{b}}}{{\rm{)}}^2}{\rm{ + (}}{H_{\rm{a}}}{\rm{sin}}{g_{\rm{a}}} - {H_{\rm{b}}}{\rm{sin}}{g_{\rm{b}}}{{\rm{)}}^2}{{\rm{]}}^{1/2}}, $

(9)

其中, H_a和g_a为卫星高度计数据或者插值所得的振幅和迟角, H_b和g_b为验潮站资料提取的振幅和迟角或者全球潮汐TPXO7.2模式的振幅和迟角。矢量差Δ、振幅差ΔH和迟角差Δg的平均值分别用$\bar \Delta $、$\Delta \bar H$和$\Delta \bar g$表示。

2 结果分析 2.1 星下观测点处调和常数

参照Fang等^[4]和暴景阳^[9]中所述方法, 计算了M₂、S₂、K₁和O₁四个主要分潮的混淆周期, 结果见表 1。然后根据Rayleigh准则, 计算得到了分离四个主要分潮所需的最短资料长度为2.97 a(表 2)。剔除观测时间小于2.97 a的星下点数据, 得到了可以分离出四个主要分潮调和常数的4309个星下点(图 1)。

利用调和分析方法得到了4 309个星下观测点处M₂、S₂、K₁和O₁四个主要分潮的调和常数(注:下文的迟角都为东8区迟角)。以15 km的空间窗口对4 309个星下观测点与南海141个验潮站的位置在空间上进行了匹配, 其中有40个验潮站符合要求。

对比星下观测点处调和分析的结果和验潮站处观测结果(表 3)可知, 利用卫星高度计数据所得星下点处的四个主要分潮调和常数结果与验潮站处观测结果误差较小, 半日分潮的振幅绝均差整体大于全日分潮的振幅绝均差。其中, M₂分潮的振幅绝均差和矢量均差均最大, 分别为3.16 cm和4.99 cm; K₁分潮与验潮站处观测结果的振幅绝均差最小, 为1.75 cm; S₂分潮与验潮站处观测结果的矢量均差最小, 为3.43 cm。

通过与前人结果进行对比(表 4)可知:相对于李培良等^[10], 本文调和分析所得星下观测点处K₁分潮的振幅和迟角与验潮站处观测结果的误差均相对较小; M₂、S₂和O₁分潮振幅和迟角拟合效果与前人结果相当, 其中S₂分潮迟角绝均差相对较大, 但仍与暴景阳等^[9]、王延强等^[11]和Fang等^[12]在南海的误差结果相近。

总的来看, 调和分析所得星下观测点处K₁分潮的误差相对最小, M_2、S₂和O₁分潮结果总体偏差与前人研究相当, 上述结果表明:利用1992~2016年的卫星高度计数据进行调和分析所得结果是可信的。

2.2 南海同潮图 2.2.1 双调和样条插值所得同潮图

以美国国家海洋大气局发布的水深数据ETOPO5为参考, 制作了分辨率为5′×5′的插值参考点。将调和分析所得M₂、S₂、K₁和O₁四个主要分潮的振幅H和迟角g转化为Hsing和Hcosg, 然后采用双调和样条插值法对其进行插值, 再将插值所得网格点处的Hsing和Hcosg转化为振幅H和迟角g, 进而绘制了四个主要分潮的同潮图(图 2)。

2.2.2 与TPXO7.2模式同潮图的对比

TPXO7.2是由俄勒冈大学开发的全球海洋潮汐模型, 该模式以拉普拉斯潮汐方程组为基础, 同化了T/P、Jason卫星高度计和验潮站数据^[14]。高秀敏等^[14]采用南海60个验潮站的调和常数资料对TPXO7.2模式结果进行了评估, M₂、S₂、K₁和O₁四个主要分潮的拟合程度均大于83%。此外, TPXO7.2模式结果与本文141个验潮站处M₂、S₂、K₁和O₁四个主要分潮的矢量均差分别为: 12.68、6.32、9.32和7.23 cm。上述结果表明: TPXO7.2模式结果在南海有较高的准确度, 可用于评估双调和样条插值方法所得南海同潮图。

利用TPXO7.2模式结果绘制了南海四个主要分潮的同潮图(图 3), 并将其作为基准, 与插值所得同潮图进行了对比。从同潮图整体分布来看, 双调和样条插值所得四个分潮的同潮图在南海中部、南海南部和泰国湾等绝大部分研究区域与TPXO7.2模式所得同潮图的等振幅线和等迟角线的变化趋势较为接近; 在雷州半岛附近海域、台湾海峡和泰国湾湾顶等少数浅水海域与TPXO7.2模式的同潮图存在一定的误差。

为了定量计算插值结果和TPXO7.2模式结果之间的误差, 将南海48 500个网格点处双调和样条插值所得调和常数和TPXO7.2模式结果进行了比较, 四个主要分潮的矢量差在南海的分布见图 4, 整体误差统计结果见表 5。可以看出, 四个分潮的矢量差在远离陆地的开阔海域都小于4 cm, 矢量差大于10 cm的网格点大多分布在雷州半岛附近海域和台湾海峡等少数近岸浅水海域, 这可能与TPXO7.2模式在我国近海海域精度较差有关^[15]。车贝雪夫多项式广泛应用于气象要素插值^[16], 车贝雪夫多项式插值所得的O₁分潮的同潮图与TPXO7.2模式结果误差最小, 振幅绝均差和迟角绝均差分别为: 4.49 cm和17.58°, 而双调和样条插值所得O₁分潮的同潮图与TPXO7.2模式结果振幅绝均差和迟角绝均差分别为: 1.43 cm和6.06°, 双调和样条插值所得结果明显优于车贝雪夫多项式插值结果。此外, 双调和样条插值所得M₂分潮的同潮图与TPXO7.2模式结果误差最大, 振幅绝均差和迟角绝均差分别为: 3.06 cm和7.37°, 但也在合理范围之内。以上结果表明:双调和样条插值结果和TPXO7.2模式结果之间整体误差较小。

2.2.3 与验潮站数据的对比

为了进一步验证利用双调和样条插值法所得同潮图的精度, 比较了141个验潮站处双调和样条插值所得四个主要分潮调和常数和验潮站观测结果, 四个主要分潮振幅插值结果与验潮站观测结果的对比见图 5, 振幅绝均差、迟角绝均差和矢量均差结果见表 6。可以看出, M₂、S₂、K₁和O₁四个主要分潮的矢量均差分别为22.59、10.26、10.24和8.51 cm, 平均相对误差均小于36%。车贝雪夫多项式插值所得M₂、S₂、K₁和O₁分潮的矢量均差分别为46.15、14.93、16.60和14.83 cm, 均大于双调和样条插值方法所得误差, 进一步说明双调和样条插值所得结果优于车贝雪夫多项式插值结果。以上结果说明利用双调和样条插值所得结果与利用验潮站观测资料提取的调和常数误差较小, 能够较好地刻画四个主要分潮的特征。

李培良等^[10]利用正交潮响应法对南海T/P高度计资料进行了潮波分析, 并按照潮汐在海洋内部的变化向岸边插值, 与53个验潮站处的主要分潮数据进行了对比, 结果表明: K₁分潮振幅和迟角的平均绝对误差分别是4.73 cm和11.6°, M₂分潮振幅和迟角的平均绝对误差分别是11.91 cm和28.4°。高秀敏等^[17]利用一个考虑内潮耗散的二维潮波模型模拟了南海M₂、S₂、K₁和O₁四个主要分潮, 并与南海63个验潮站和24个T/P卫星高度计轨道交叉点处的调和常数进行了对比, M₂、S₂、K₁和O₁四个主要分潮的振幅绝均差分别为: 13.57、3.09、12.23和3.81 cm, 迟角绝均差分别为: 17.41°、13.92°、24.97°和19.27°。可以看出, 利用双样条调和插值方法所得南海四个主要分潮的误差与前人利用正交潮响应法或者数值模拟方法得到的误差相近, 但本文所采用的双调和样条插值方法计算量更小、更易实施, 是一种有效获取南海同潮图的方法。

2.2.4 无潮点位置的对比

本文得到的M₂、S₂、K₁和O₁四个主要分潮的无潮点均分布在泰国湾(具体位置见表 7), M₂分潮在泰国湾有两个旋转潮波系统, 湾口的无潮点呈顺时针旋转, 潮波系统特征和无潮点位置与Fang等^[18]在泰国湾的潮汐数值拟合结果相近, 湾顶的逆时针无潮点离陆地很近, 位置与曹德明等^[19]的结果最接近(见表 7), S₂、K₁和O₁分潮均在泰国湾有一个无潮点, S₂分潮的潮波系统呈顺时针旋转, K₁和O₁分潮的潮波系统都是逆时针方向旋转, 三者的无潮点位置均与前人的研究结果基本一致。M₂、S₂、K₁和O₁四个主要分潮无潮点位置的比较结果表明:双调和样条插值所得结果不仅能够在整体上刻画南海四个主要分潮的特征, 还能较为准确地再现无潮点位置这一局部特征。

3 结论

基于1992年至2016年卫星高度计海面高度资料, 运用调和分析方法得到了南海区域卫星星下观测点处M₂、S₂、K₁和O₁四个主要分潮的调和常数, 并通过和验潮站观测数据以及前人研究结果的对比验证了调和分析结果的准确性。

利用双调和样条插值方法对星下观测点处四个主要分潮的调和常数进行插值, 进而绘制了南海M₂、S₂、K₁和O₁四个主要分潮的同潮图。插值结果在在雷州半岛附近海域和台湾海峡等少数近岸浅水区域与TPXO7.2模式结果误差较大, 在其他海域均与TPXO7.2模式结果误差较小; 与141个验潮站处观测结果的矢量均差分别为22.59、10.26、10.24和8.51 cm; 所得四个主要分潮的无潮点数量和位置与前人研究结果基本一致。此外, 利用双样条调和插值方法所得南海四个主要分潮的误差与前人利用正交潮响应法或者数值模拟方法得到的误差相近。上述结果表明:双调和样条插值方法是一种有效的插值方法, 能够通过对卫星高度计所得调和常数数据进行插值获取较为准确的南海同潮图。