研究生数学建模竞赛优秀论文选《基于卫星调度计海面高度异常资料获取潮汐调和常数方法650.docx
《研究生数学建模竞赛优秀论文选《基于卫星调度计海面高度异常资料获取潮汐调和常数方法650.docx》由会员分享，可在线阅读，更多相关《研究生数学建模竞赛优秀论文选《基于卫星调度计海面高度异常资料获取潮汐调和常数方法650.docx（69页珍藏版）》请在冰豆网上搜索。
研究生数学建模竞赛优秀论文选《基于卫星调度计海面高度异常资料获取潮汐调和常数方法650
“HW杯”第十五届中国研究生数学建模竞赛
学校
参赛队号
1.
队员姓名
2.
3.
“HW杯”第十五届中国研究生数学建模竞赛
题目基于卫星高度计海面高度异常资料获取潮汐调和常数方法及应用
摘要：
利用T/P卫星高度计资料分析海洋潮汐潮流现象，是动力海洋学研究的重要内容。
在描述潮汐潮流特征过程中，潮汐调和常数的获取是重要的一项科学研究，直接影响海洋潮汐同潮图的绘制，为波浪、风暴潮、环流、水团等其他海洋现象的研究进展提供资料，是加强国防建设、促进海洋能源开发、保护环境、建设海港和保护海岸中不可或缺的一部分。
本文基于卫星高度计资料探究潮汐调和常数计算方法及其应用，主要解决以下问题：
问题1针对整个卫星下观测点的海面高度异常值，利用调和分析方法对各个主要分潮（M2、S2、K1、O1）的潮汐调和常数进行了求取，求解过程中，基于达尔文-杜森模型，选取九个分潮（M2、S2、K1、O1、N2、K2、P1、Q1、Sa）进行数据拟合，并结合最小二乘法进行求解，最后通过验潮站数据和上下行轨道交叉点对结果进行验证：
与20个验潮站调和常数结果相比，四个主要分潮绝均差分别为：
M2、S2、K1、O1分潮振幅绝均差为16.71、7.33、15.89和12.21cm；迟角绝均差为6.94º、10.84º、9.15º和9.72º。
20个交叉点的各个分潮振幅和迟角绝均差：
M2、S2、K1、O1振幅差ΔH均方根为0.65、1.09、0.83、0.81cm，矢量差Δ均方根为0.85、1.02、0.79、0.78。
问题2通过绘制四大主要分潮（M2、S2、K1、O1）潮汐调和常数沿卫星轨道
的分布图，发现潮汐调和常数分布在空间上存在由内潮对正压潮调制引起的细节构，这主要表现为局部区域存在较大的数值波动。
鉴于此，本文考虑使用多项式拟合方法，拟合正压潮的调和常数，并在此基础上，通过坐标变换获得内潮的调
和常数，实现正压潮与内潮的分离，并分析了分离后的正压潮的调和常数的分布规律，得到同潮图插值经验规则。
问题3进一步针对南海所有卫星观测轨迹下计算得到的正压潮调和常数，利
用克里金插值法分析各个分潮上的同振幅、同迟角的位置，并结合问题2中的经
验插值规则，从而获得各主要分潮正压潮的同迟角与同幅值分布图。
通过与验潮站数据进行对比分析，结果表明：
8大分潮正压潮幅值、迟角计算结果与验潮站对比的绝均差分别为：
M2、S2、K1、O1、N2、K2、P1和Q1分潮振幅绝均差为9.68、15.22、10.94、12.59、13.64、12.86、11.06和17.41cm；迟角绝均差为8.81º、
8.65º、11.81º、8.85º、9.45º、8.53º、10.35º和7.35º。
问题4针对前两问拟合和插值过程中可能存在的拟合残差较大的问题，本部
分主要从两个方面出发：
1）针对不同次数下的多项式拟合分潮调和常数进行正压潮和内潮的分离分析讨论；2）针对不同插值方法和方式的同潮图绘制的分析
讨论。
得到如下启示：
1）可以分上下行轨道分方法插值，设置权重，增加下行轨道对插值结果的权重比例，降低上行轨道振幅数据对插值结果的影响，从而可以获得较为完备的各个分潮的同潮分布图；2）依据不同分潮类型的实际情况，进行不同的插值方式研究，最后进行叠加，从而达到完整的同潮图绘制。
关键词：
卫星高度计；海面高度异常；潮汐调和常数；最小二乘法；调和分析
一、问题背景及重述
1.1问题背景
海洋潮汐是在天体引潮力作用下形成的长周期波动现象，在水平方向上表现
为潮流的涨落，在铅直方向上则表现为潮位的升降。
潮汐潮流运动是海洋中的基本运动之一，它是动力海洋学研究的重要组成部分，对它的研究直接影响着波浪、风暴潮、环流、水团等其他海洋现象的研究，在大陆架浅海海洋中，对潮汐潮流的研究更具重要性。
海岸附近和河口区域是人类进行生产活动十分频繁的地带，而这个地带的潮汐现象非常显著，它直接或间接地影响着人们的生产和生活。
潮汐潮流工作的开展和研究，可为国防建设、交通航运、海洋资源开发、能源利用、环境保护、海港建设和海岸防护提供资料。
例如，沿海地区的海滩围垦、农田排灌，水产的捕捞和养殖，制盐，海港的选址及建设，以至于潮能发电等活动，无不与潮汐潮流现象有着密切的关系。
区域海洋潮汐的数值模拟需要提供开边界的水位调和常数，而开边界的水位调和常数，或者来源于观测、或者来源于全球海洋潮汐的数值模拟；而全球海洋潮汐的数值模拟，相当耗费资源。
虽然目前有国外学者或研究机构，能够提供区域海洋潮汐的调和常数，但实质上的评价结果难以令人满意。
从区域海洋潮汐的数值模拟的现状来讲，四个主要分潮（M2、S2、K1、O1）的单一分潮的数值模拟与同化可以得到令人满意的结果，但其它分潮（N2、K2、P1、Q1等）的单一分潮的数值模拟与同化，结果却差强人意；这意味着其它分潮的数值模拟，只有与四个主要分潮同时进行数值模拟，才能得到可以接受的结果。
从具体操作来讲，其它分潮由于相对较弱，导致模拟结果的精度难以提高。
长周期分潮（Sa、Ssa、Mm、Mf）的获取，目前已有基于全球长周期分潮数值模拟手段的报道，但其面临的困境，与其它较弱分潮面临的困境没有差别。
从各分潮的调和常数获取的发展史来说，通过对已有观测结果进行插值曾经
是首选，但发展过程中逐渐被数值模拟方法所取代。
高度计资料的出现，引发部分学者开展了插值方法的研究，并取得了一些值得一提的结果，尽管被所谓的主流方式淹没，但也难掩其光芒所在。
鉴于目前已有高度计资料作为支持，其它分潮及长周期分潮的调和常数获取的插值方法研究大有可为。
1.2资料描述
1.2.1地形数据
地形数据来自ETOP5，全球的分辨率为5'⨯5'，图1.1的区域是2~25N，
99~122E。
1.2.2验潮站资料
图1.1南海地形图
中国近海及周边海域770个验潮点的资料，和56个验潮点的资料（是国际
上公开的长期验潮站数据分析得到的调和常数），包括9个分潮（M2、S2、K1、
O1、N2、K2、P1、Q1、Sa）的潮汐调和常数。
图1.2显示了上述资料点所在的位置，从图中可以看出上述验潮点主要分布在近岸或岛屿附近。
图1.2验潮站资料的分布图
1.2.3TOPEX/POSEIDON卫星高度计资料
TOPEX/POSEIDON卫星是由美国国家航空航天局和法国空间局联合于
1992年8月10日发射的，是世界上第一颗专门为研究世界大洋环流而设计的高度计卫星。
其轨道高度达1336km，倾角为66°，覆盖面大，保证了资料的连续性。
轨道的交点周期（绕地球一圈的时间）为6745.8s，轨道运行127圈以后精确重复，轨道重复周期为9.9156天。
相邻最近的轨道之间在赤道上的间隔为360/127=2.835。
卫星在一个周期内的每一圈分为上行轨和下行轨两条轨道，一
个完整的周期内共有254条轨道，沿轨道的两个相邻的星下观测点的距离
5.75km。
高度计系统的定规精度和测高精度较以前有显著提高，其测量精度约为
5cm，是目前观测海面高度精度最高的卫星。
当然，本文只是涉及到TOPEX/POSEIDON卫星高度计资料与潮汐相关的研究，即海面高度异常产品。
1.2.4南海高度计资料
图1.3的给出了2︒~25︒N，99︒~122︒E，TOPEX/POSEION卫星高度计星下
观测点所在的轨道。
一共有超过4000个数据点，每个点都对应一个海面高度异
常的时间序列，从1992年到2017年，时间跨度为25年。
1.3问题提出
图1.3南海TOPEX/POSEIDON高度计资料的星下轨迹
1.3.1问题1：
潮汐调和常数求取与评价
根据沿轨道的星下观测点的海面高度异常值，提取所有星下观测点各主要分
潮（M2、S2、K1、O1）的潮汐调和常数，注意能有效提取那些分潮的潮汐调和常数取决于相应的资料长度；对提取的潮汐调和常数，应利用潮汐验潮点的调和常数给予评价或检验，并给出评价结果的分析或评价。
1.3.2问题2：
正压潮和内潮的分离
得到所有星下观测点各主要分潮（M2、S2、K1、O1）的潮汐调和常数，沿
轨道作图后，可发现潮汐调和常数在沿轨道方向，在空间有细结构，而此细结构是内潮对正压潮的调制；请设法对沿轨道的各分潮的潮汐调和常数进行正压潮和内潮的分离。
1.3.3问题3：
分潮特征同潮图分析
设计数据插值或拟合方法给出南海的各主要分潮的同潮图，并利用潮汐验潮
点的调和常数给予评价或检验，并给出评价结果的分析或评价。
1.3.4问题4：
潮汐调和分析优化过程讨论
如果你们还有时间和兴趣，还可考虑下列：
如果在对沿轨道的潮汐调和常数分离、插值或拟合的过程中，利用了特定的函数进行拟合，是否能够确定出需利用的特定函数的最佳（高）次数？
上述结论是否对第3问有启示或帮助。
二、问题分析
2.1问题1：
潮汐调和分析和评价方法
问题1中主要针对沿轨道的星下观测点的海面高度异常值，需要提取星下测点各个主要分潮的潮汐调和常数。
需要对潮汐调和常数的获取，进行调和分析，利用最小二乘法求解线性线性方程组，最终需要借助可行的评价方法对求得的分潮调和常数进行对比和验证，本文主要利用上下行轨道交叉点调和常数和验潮站的调和常数分别进行验证和对比分析。
2.2问题2：
沿轨道正压潮和内潮分离方法
问题2中主要利用沿轨道做出的各主要分潮的潮汐调和常数（分潮振幅和迟角）分析图，得到潮汐调和常数在轨道方向上，有内潮对正压潮调制引起的细结构这一结论。
鉴于此，需要对正压潮和内潮进行分离。
这方面的研究主要为特定形式的多项式拟合方法，对正压潮和内潮（斜压潮）时间序列进行调和分析，
分离出正压潮和内潮的主要成分即调和常数参量。
从而得到较为平滑的正压潮调
和常数参量，为第三问同潮图的绘制提供基础资料。
2.3问题3：
分潮特征同潮图分析
问题3主要是设计数据插值或拟合方法以得到南海的各主要分潮的同潮图。
利用第二问中分离得到的平滑的正压潮的潮汐调和常数的参量，结合自然邻居插值方法，分别对各个主要分潮（M2、S2、K1、O1）和另外4个分潮（N2、K2、P1、Q1）的调和参数进行插值，从而绘制出其正压潮振幅和迟角的同潮分布图，最后利用验潮站的数据，对插值后得到的靠近验潮站的数据进行评价和检验。
2.4问题4：
潮汐调和分析优化过程讨论
问题4主要针对沿轨道的潮汐调和常数分离、插值或拟合过程中，使用了特定格式的拟合函数，主要从两个方面进行分析和讨论：
1）对于分离过程中，拟合正压潮的多项式，分析了5、7、9、11、13次多项式的拟合结果，优选出（多项式最佳形式）最佳次数下的多项式；2）从空间插值法进行改进，利用部分加权的方法对时间序列差别较大的数据设置权值（例如：
上下行轨道不同权值插值），从而利于各个分潮同潮图的特征展示。
图2.1研究技术路线图
三、模型假设
为了使实际问题的解决具有可行性，同时易于数学模型的建立，在针对性建模之前需要做一定的假设。
（1）对于T/P卫星高度计的测量误差，符合正态分布，不存在测量错误；
（2）在分析实际潮汐的分潮过程中，不考虑非天文分潮；
（3）逐时观测记录完整并且连续，允许潮时误差；
（4）潮高水位平均水位可取近似值；
（5）潮汐调和常数的差比关系用近似代替。
四、名词解释及符号说明
表4.1符号说明表
符号
意义
符号
意义
R
卫星和海表面的距离
N（z）
垂直浮力剖面
∏（z）
水平速度的垂直结构
Φ（z）
内波引起的垂直位移
HKE
水平动能
gj
迟角
hj
振幅
fj
交点因子
μn
Doodson数
TA
混淆周期
ς（ti）
第i时刻的实际潮位
G格
格林尼治时间
σ
分潮的角速率
∆
上升与下降轨在交点处的调和常数的矢量差异
r
噪声
DC
确定性系数
li（t）
拉格朗日立潮位多项式
H2n-1（x）
三次Hermite插值立潮位多项式
S3（x）
三次样条插值立潮位多项式
TsatPeriod
采样间隔
Hk
内潮海表面振幅
gk
内潮海表面迟角
五、模型的建立与求解
5.1问题一：
潮汐调和分析方法
5.1.1卫星测高分潮混叠
执行重复轨迹任务的测高卫星，连续两次通过某上升轨迹或下降轨迹的时间
间隔为其重复周期，也即对海面高的采样间隔TsatPeriod。
根据Nyquist采样定理，当采样时间间隔大于信号的半周期时，会产生频率混叠效应。
测高卫星的时间序列可恢复信号的最高频率为Nyquist频率，为1/（2TsatPeriod），对应的信号周期为2TsatPeriod。
对于T/P卫星而言TsatPeriod、为9.9156天，在此期间卫星绕地球旋转127圈，其结果是直接探测信号的周期在20天以上，远大于潮汐主要半口分潮和口分潮的周期，其相应的折叠频率远低于主要潮汐频率，因而其采样必然产生潮致高频混淆。
折叠频率对应的最小周期为混淆周期，对于周期为TTidePeriod（小时）的分潮而言，分潮的混淆周期TAliasPeriod为：
t=TsatPeriod⨯24/TTidePeriod
（5.1）
t=t-[t]
（5.2）
⎨
[t]表示不大于t的最大整数，且t=⎧tt<0.5。
⎩1-tt>0.5
所以对于分潮TTidePeriod而言，其混淆周期TAliasPeriod为：
TAliasPeriod=TSatPeriod/t
（5.3）
表5.1主要分潮在T/P海面轨迹上采样的混叠时间
分潮
分潮周期（h）
潮汐观测的混叠周期（d）
M2
12.42
62.1
S2
12.00
58.7
N2
12.66
49.5
K2
11.97
86.6
O1
25.82
45.7
P1
24.07
88.9
K1
23.93
173.2
Q1
26.87
69.3
与以往卫星相比，随着设备精度的提高和校正工作的改进，T/P卫星的径向轨道误差有大幅度降低，且其轨道设计使8个主要分潮的混叠周期均低于半年。
对于混淆频率为T1和T2（T1>T2）的两个潮汐信号，实现可靠分离所需要的采样时间Td为：
Td=T1⨯T2/（T1-T2）/4
（5.4）
上式中“4”，表明要区分这两个分潮，采样周期至少要是信号周期的四分之一。
由此计算的在T/P采样情况下8个主要潮汐分潮之间的分辨时间如下表5.2所示。
表5.2主要分潮基本分辨时间/天
M2
S2
N2
K2
O1
P1
K1
Q1
M2
0
271.7
61
54.9
43.3
51.5
24.2
148.8
S2
0
78.8
45.7
51.5
43.3
22.2
96.1
N2
0
28.9
148.8
27.9
17.3
43.3
K2
0
24.2
844.2
43.3
87
O1
0
23.5
15.5
33.5
P1
0
45.7
78.8
K1
0
28.9
Q1
0
从上表可以看出，区分Ma和Sa分潮，需要0.7年（271.7天）的数据；区分Ka和Pi分潮，需要约2.3年（844.2天）的数据，该时间是8个主要分潮中分辨两个分潮所需要的最长时间。
因此，要获取M2、S2、N2、K1、O1、Qi、Ka和Pi稳定的调和常数，至少要100个采样点（约1000天）。
因此，本文采用的T/P卫星高度计时间序列的长度足够对以上主要分潮进行调和分析。
5.1.2利用达尔文-杜森模型（f/u）进行潮汐的调和分析
实际潮汐的分潮从其来源看可分为以下四种：
天文分潮、气象分潮、天文-
气象分潮和浅水分潮。
从分潮的频率分布来看，分潮在频率上的分布是极不均匀的，而是分成族、群和亚群。
在Doodson展开中，按Doodson数μ1区分潮族，按μ2区分群，按μ3区分亚群。
在潮族中一般分为长周期分潮族（μ1=0）、全日分潮族（μ1=1）、半日分潮族（μ1=2）、三分日分潮族（μ1=3）直到十二分日分潮族（μ1=12），共13个潮族。
在每一个潮族中，具有不同数量的群和亚群。
在亚群中的各个分潮的角速度是非常接近的，彼此之间只有微小的差异。
因此，在资料长度有限的情况下，亚群中的各个分潮是无法区分的。
因此，在实际的潮汐分析中，往往将一个亚群合成一个分潮，此时这一分潮的振幅和迟角不再是常数，而是随着升交点的黄经十分缓慢地变化，一般在较短的时间内可近似看作不变。
这样的分潮实质上是准调和的，但习惯上仍叫做调和分潮。
实际水位可以看作是很多个调和分潮迭加的结果，但是在实际分析中只能选取其中有限个较主要的分潮。
假设我们选取了J个分潮，对于任一点的潮位表达式为：
JJ
h=S0+∑fjhjcos（vj+uj-gj）=S0+∑fjhjcos（σjt+v0j+uj-gj）
j=1j=1
（5.1）
其中，S0为余水位，fj为交点因子，uj为交点订正角，hj,gj为分潮的调和常数（振幅和迟角）。
（1）分潮角速度的计算
∙∙∙∙∙∙
ο=μ1τ+μ2s+μ3h'+μ4p+μ5N'+μ6p'
（5.2）
其中：
σ为分潮的角速度，μ1,μ2,μ3,μ4,μ5,μ6为Doodson数，
⎧
∙
⎪τ=14.49205211
⎪∙
⎪s=0.54901653
⎪h
⎪∙'=0.04106864
⎨∙
⎪p=0.00464183
⎪
⎪N
⎪∙'=0.00220641
⎩p
⎪∙'=0.00000196
（2）分潮初相位的计算
（单位：
度/平太阳时）
Y年M月D日t时刻（实际计算中是观测数据的起始时间）的天文初相角：
v0=μ1τ+μ2s+μ3h'+μ4p+μ5N'+μ6p'+μ090
（5.3）
其中：
μ0,μ1,μ2,μ3,μ4,μ5,μ6为Doodson数，
⎧s=277.02+129.3848（Y-1900）+13.1764（n+i+t）
⎪24
⎪t
⎪h'=280.19-0.2387（Y-1900）+0.9857（n+i+）
⎪24
⎪t
⎪p=334.39+40.6625（Y-1900）+0.1114（n+i+）
⎨24
⎪t
⎪N'=100.84+19.3282（Y-1900）+0.0530（n+i+）
⎪24
⎪p'=281.22+0.0172（Y-1900）+0.00005（n+i+t）
⎪24
⎪=15t-s+h'
⎪⎩τ
（5.4）
式中i为1900年至Y年的闰年数，i=int（Y-1901）；n为从Y年1月1日开
4
始计算的累积日期序数，1月1日的日期序数为0，t为时间（单位：
小时）。
以上各式中的单位是度。
表5.3部分分潮的Doodson数、分潮角速度和交点因子与订正角
分潮符号
Doodson数
分潮角速度
单位：
度/平太阳时
交点因子与订正角
μ1
μ2
μ3
μ4
μ5
μ6
μ0
f
u
Sa
0
0
1
0
0
0
0
0.0410686
1
0
Ssa
0
0
2
0
0
0
0
0.0821373
1
0
Mm
0
1
0
-1
0
0
0
0.5443747
Mm
Mm
MSf
0
2
-2
0
0
0
0
1.0158958
M2
-M2
Mf
0
2
0
0
0
0
0
1.0980331
Mf
Mf
Q1
1
-2
0
1
0
0
-1
13.3986609
O1
O1
分潮符号
Doodson数
μ1μ2μ3μ4μ5μ6μ0
分潮角速度
单位：
度/平太阳时
交点因子与订正角
fu
O1
1
-1
0
0
0
0
-1
13.9430356
O1
O1
M1
1
0
0
0
0
0
1
14.4920521
M1
M1
P1
1
1
-2
0
0
0
-1
14.9589314
P1
P1
S1
1
1
-1
0
0
0
2
15.0000000
1
0
K1
1
1
0
0
0
0
1
15.0410686
K1
K1
J1
1
2
0
-1
0
0
1
15.5854434
J1
J1
OO1
1
3
0
0
0
0
1
16.1391017
OO1
OO1
N2
2
-1
0
1
0
0
0
28.4397295
M2
M2
M2
2
0
0
0
0
0
0
28.9841042
M2
M2
L2
2
1
0
-1
0
0
2
29.5284789
L2
L2
S2
2
2
-2
0
0
0
0
30.0000000
1
0
k2
2
2
0
0
0
0
0
30.0821373
k2
k2
M4
4
0
0
0
0
0
0
57.9682085
M2
2
2M2
MS4
4
2
-2
0
0
0
0
58.9841043
M2
M2
M6
6
0
0
0
0
0
0
86.9523127
M3
2
3M2
表中交点因子及交点订正角的含义说明：
例如，表中M6分潮的交点因子是M2分潮的交点因子的三次方，M6分潮的交点订正角是M2分潮的交点订正角的三倍。
（3）fj和uj的计算
由于fj和uj随时间变化非常缓慢，一般情况下取资料序列的中间时刻计算。
各分潮的
升级会员 