用MATLAB的T_TIDE程序_省略_行包含误差估计的经典潮汐调和分析_Ric

【实例简介】
用MATLAB的T_TIDE程序_省略_行包含误差估计的经典潮汐调和分析_Ric
直接按照矩阵的形式而写,容易理解和修改。最后,为了区分頁实确定性(线性)频率和 广谱变化,所估计的潮汐参数的置信区间是在几种用户可选择的算法中选取一种进行计算 这个程序包由很多文件组成,每个文件里有一个或多个函数。用户可调用的函数一般有“t” 前缀,以避免命名冲突 本文由4部分组成。第一部分是对平衡位的形式进行描述;第二部分措述的是估计振 幅/相位的数学基础;第三部分是说明置信区间怎样产生:;最后一部分是举例进行讨论。 2.潮汐位 月亮和太阳产生的万有引力矢量F可表小为一个标量位V的梯度,F=-VV。在任何 时刻,地球表面的潮汐位大小都取决于地球、太阳、月亮的相对位置。在 Doodson的研究 结果中( Doodson,1954,见 Godin,1972),这个位被表示为月亮时τ(定义为从“月亮午 夜”起始)和下列天文变量(都是时间的函数)的函数 S是月亮的平经度,h是太阳的平经度,p是近地点的经度,N是升交点黄经的负数, p'是近日点经度,都以周期为单位。如果凵知儒略目( Julian date),这些变量叮以用 t astron 函数进行计算,公式见 Scidclmann(1992)。对于一个特殊的成分{,’k',r,m,n},它们 的结果和 Doodson数集结合成个天文参数V=+j+h+p+mW+np't是波类码 (l=0,1,2分别对应长周期波,日波,半日波),ljk'为分波码。成份频率σ定义为 σ=2ml/dt。潮汐位被写作: y=∑G(0)∑ ArrKIm'n, coS((2m)+G1(0)∑ Brim'n,sin(2m)(1) kr'm'n 对于一个给定的 Doodson数集,A"或B'是非零的,但不会同时非零。这些常数被制作成 表格并存储在一个数据结构里,可以用 t gctconsts程序加载。大地测量函数G和G1,随着 和纬度θ变化,同时也取决于地球的半径,地球、月亮和太阳的质量以及三者问的距离。 每个成份的平衡振幅定义为GA/g或者是G"B'/g。式中g是重力加速度,对于特定的纬 度可以用 t equilib程序得出 3.相位/振幅估计 这里采用的估计相位/振幅的算法是基于 Godin(1972)、 Foreman(1977和 Foreman(1978) 等人用 Fortran语言编写的算法。和那些作者不同,我们是直接使用复数而不是分别用正弦、 余弦拟合。这样做的优点是能统一对标量(如压力)和矢量(如水平流)时间序列的处理, 用复数u+ⅳ表示。要注意,潮流的复薮表示形式是基于一个潮流旋转矢量的物理模型,只 对线性或者接近线性的潮汐波满足。在某些情况下,将沿着和垂直通过通道的潮流看作两 个独立的标量时间序列比较好。 考虑一个时间序列观测值y(t),t=12t2,t为一列矢量。观测时间间隔An(默认为 1小时),M是一个奇数(必要时可以去掉一个端点)。时间轴的原点(或中央标准时间) 为1(M+D2缺失的观测值可以在输入矢量中用“缺失数据”标记处理(在 MATLAB中用 19 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct NaN表示,即IEEE算法里的Not-a- Number)。这种常规的吋间间隔约束不是由最小二乘拟 合引起,而是由自动选择成份的算法引起,在置信区闫算法估算频谱时也是需要的。这个 时间序列可以是实数也可以是复数。它和其他(多半是任选的)参数一起进入ttde分析 程序 潮汐响应被拟合为 x()=b+b1t+∑a +aie k=1,V 式中采用了N个成份(每种成份有特定的 Doodson数集)。随着潮汐位的发展,每个成份 的频率a是已知的,而复振幅ak未知,尽管时间序列y(4)是实的时间序列,ak和ak是 复共轭。表达式的前两项表示可能的偏移和(任选的〕线性漂移。传统的方法是用实的正 弦形式表示: x()=6o+b,t+ 2A cos(or t)+ Bk sin(o t k=l.N 通过A=a1+ak、B1=i(a-aA)和方栏(2)相关联。实数表达式在下面将提到的线性 误差分析里会更方便。 从45个天文成份和101个浅水成份里挑选所需要的成份。数据结构包含这些成份的名 字和其他信息,可用 t getconsts进行加载。选择成份可用很多种算法,一般釆用自动挑选 算法( Foreman,1977),运作如下。将所有主要的天文成份和24个最重要的浅水成份收集 在一起ε根据预先按平衠振幅确定的重要性次序,将这些成份一一列岀。除了比较重要的 频率成分外,频率小于 Rayleigh分辨极限a(M)的次要成份(默认为a=1)要去掉。必 要时可以补充规定一些“浅水成份”。如果两种成份的相对相位/振幅是通过其他资料知道 的,可以釆用推理程序。在其他情况下,成分列表规定待很清晰。 最小乘法拟合是将下式最小化的系数 E=(n)-y(n)2=a-)2 式中y=[y(41),y(2)…y(t),a=[b,b1,a1a1…,aN],T是Mx(2N+2)矩阵,是观 测时间得到的线性正弦基函数。解是用 MATLAB矩阵除法运算符待出的 在拟合之后可以通过tw进行各利校正。第一,各种成份响应的相位记作“格林威治 相位g,也就是参照经度为0°(格林崴治子午线)时的平衡响应的相位。这也可以破认为 是在0经度处当平衡加压处于最人正值时的响应的相位。最简单的是找出记录的中央标准 时间(t=0)的拟合相位;此时对于在相应的儒略日计算的给定成分,平衡相位v刚好等 于V,但可能婁作14、1/2或34周期的调整,视A或B是否非零及其正负号而定 第二,如果纬度给定,按照如下方式进行节点和卫星校正。考虑在下标为矿的卫星对 下标为k的成份的主峰值的影响。不同卫星的作用是在不同阶段缓慢地调制主峰值的相位 和振幅,通常超过8年。在一些时间段内,我们的拟合响应ak可以被写为主成分“真实” 响应ak的词制,意味着由于有卫星,振幅和相位分别通过一个因子f和角度lk发生了改 变。f和u分别称作节点校正振幅和相位,即 20 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct =faeo=aeo+∑ (5) 消去公共项,得: k(06:-k) 1+∑ 只要(k-ak)一直很“小”(即N6t<8年),最终近似就会成立。一般而言,卫星的真实 相位和振幅是不知道的。然而,山于它们的频率和峰值的频率非常类似,我们假定真实 振幅之比与平衡响应中振幅之比是一样的,真实相位差等于平衡响应相位差。于是节点校 正便可由平衡响应方程(1)得到。由人地测量函数可得知纬度。在赤道处,G1为0,为 了避免日校正过大,采用了不是很严格的限制条件。用纬度-相关的平衡响应去预测一部分 海洋盆地的动态特征的有效性还不清楚。如果记录长度超过1年,可以将下一年的进行过 和没有进行过节点校正的分析结果加以对比,以测试程序的有效性。注意,如果时间序列 长度超过18.6年,“真实”卫星振幅/相位可以直接估计( Foreman and neufeld,1991),但 是现在 T TIDE还不能完成 上述分析结果是一对复数(ak,a},对每种成价k可能进行了节点调制。可以转换 为标准参数: (7) 6=m8(ak)+0(ak) mod1 80 (9) 8k=vk -ang(a)+Ok (10) 对于水平沏流,这些参数描述的是山速度矢量末端勾绘出的一个椭圆的特点:椭圆半长轴 和半短轴的长度(分别为L),北半短轴逆时针向正东的倾角a以及格林威治相位gk 如果lk>0(<0),椭圆轨迹是逆时针(顺时针)方向。对标量时间序列,L是振幅,,64=0 (桁圆沿正轴退化成一条直线)。注意,把定乂限制为冋北轴倾斜(通过模运算符),意味 着对椭圆成东西向的成分的分析,由于噪声的缘故,轨迹会有接近0到将近180°的倾斜。 这种明显的巨大跳跃是限制造成的,并不表示有同样巨大的物理特征的变化。 旦椭圆参数确定,就可以用于进一步的分析。可以用 t predic程序在其他时间进行 预测。 t predic程序是在时间序列的中点计算节点校正,和ttde正好相反。 4.置信区间 经典调和分析的一个缺点是,相对宽带非潮汐作用的能量而言,它不能桷定一个已知 成份代表真实潮汐能的程度。从两方面来说,这是一个非常有用的信息:首先,它可以让 我们对潮汐特征傚一个更好的估计;其次,对于不同的分析我们可以做一个定量的比较 通过两步生成置信区间。第一,我们必须估计非湖汐成份的特征或者影响a(或者是 21 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct A4,Bk)的残余“噪声”;第二,我们必须将估计值转换为置信区间,以便通过非线性映射 获取标准参数。下面我们讨论实数时间序列的情形。 4.1残余噪声(实数) 在对N点实数时间序列ν(t)进行凋和分析之后,我们检查残余序列的构成。在最简单 的情况下,残差在统计上服从晑斯分布并且在时间上不相关。如果满足该情况,则总残差 功率丹=a2=PM,式中P是双边谱密度。在一个频率间隔4=(M△)内,拟合的正 弦和余弦项的振幅(分别为A,B)将会受到未分解噪声成分引起的误差的污染。因此, o=GB=PN=2/N。一个地球物理序列的谱线不可能是平坦的,在 t tidc里有更复杂 的方法米寻求与相邻成分适合的P的局部值,它是通过对时间序列残差做谱估计(即去除 所有拟合成分之后),并在仼一成份的频率周围的窗口对频点的功率求平均值,拟合成分所 在的频点忽略不计。这里我们选择宽度为0.4cpd,中心点在1,2,3,cpd的一个窗口序列 适合于半日成份的P值可由第二个频点得到。 4.2向标准参数(实数)的转化 用线性分析可以将正弦/余弦振幅的误差转化到标准参数(振幅和相位)的误差。考虑 成份k。令5=F(AB)这些参数的非线性函数,可以是振幅或格林威治相位。如果{A.B2} 是独立随机变量,我们可以找出ξ的标准差的线性化估计值,表示为正弦振幅的标准差 OF OF )2a2+()2 (l1) aB 式中的偏导数可由方稈(7)—(10)精确推导。 非线性映射可直接采用“参数辅助程序 bootstrap’进行处理( Efron and tibshirani,193)。 在这种情况下,代码用残余变量的估计值来模拟分析的一系列实现和复制,取正弦曲线振 幅的佔计值)加入有适合变量的高斯噪声。所有的实现用方程(11)—(10)非线性地转 化为标准参数,标准误差的估计值直接由这个复制数据集获得。 旦标准误差确定,95%的置信区间可由标准方法估计得到。作为一和选择,一个信 噪比(SNR)可以基于振嗝和振幅误差比的平方计算得到。用 t synth程序进行模拟,把用 不冋噪声实现的固定数据集所作分析的变化与所估计的置信区间相比较,发现只要 SNR>10,线性程序对于实数时间序列(如潮高)是足够的,SNR低为2或3的时候,结 果可能也不坏。SNR较扃吋,非线性程序和线性程序的结果很类似;SNR较低吋,结果更 精确 4.3残差信号(复数〕 一个复残差时间序列n+i可被模拟为元自噪声,并且计算方差2,O2和协方差 σ灬。如果我们进步假定两个分量的噪声都是不相关的(≈0),那么可以采用个二元 有色噪声模型,方差分被赋值到基于局部谱密度待到的成份振幅的实部和虚部。如果怀 疑σ,≠0,时间序列应该进行坐标转换(如转換到沿通道/跨通道轴或主轴)。 4.4标准参数的转换(复数 由于A=A+i4和B=B+iB,有实部和虚部,已线性化的分析包括4个变量的函数。 偏导数解析表达式 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct OF CF F )2+()2a+(2) (12) aA aB 变得很大。假定4个变量是独立的,可能会使解形式简化 铺助程序 bootstrap也能运用丁复系数ak。但是有一个小问题,除非噪声是坏形的 σ2=σ,σn=0),否则a和aA的误差会彼此相关。因此辅助稈序需要相关噪声的 重复试验。 5.实例 Foreman(1977)捉供的一个数据集的实例小于图1和表1。这个例子被收集在数据文 1975年 Tuktoyuktuk潮高 200 210 240 250 时间/年.天 潮高和残差 210 230 240 时间/年.天 95%置信区间的分析 10 频率/cph 主要成分的相位 归270 妥堡长姿 频率/cph 图1 Tuktoyuktuk分析例子。(A)原始时间序列;(B)上面的曲线是去掉潮汐信号后的时间序 列,下面的曲线是采用显著成份合成的潮汐序列;(C)所有分析成份的振幅,显著水平95%,显 著成份用实心圆圈表示;(D)显著成份的相位,95%置信区间 23 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct 表1 Tuktoyuktuk站点观测数据分析结果 File rame: paperout. txt Date: 17-Aug-2001 Nobs= 1584, ngood=1510, record length(days =66.00 Start time:06-Ju1-197501:00:00 Rayleigh cr Greenwich phase computed with nodal corrections applied to amplitude and phaee relative tc center time x0=1.98. x trend=0 Varx)=0.2196var(xF)-0.21224yar(xre)=0.60972 Percent var predicted=25.8% Tidal amplitude and phase with 95 CI estimates Tide err Pharr C.015 61.4 0.18 C.00282 0.1561 0.083 C.034 0.12 2Q1 C.03571 82,69 105.21 0.31 C.03722 C.03873 0.0764 74.23 43.35 I01 C,04027 0.0290 0.035 238,14 74,68 0.69 C.04155 0.0465 7188 70.96 C.0178 0.1405 6481 23.9 5.7 C.0432 0.025 0.050 7.32 129.76 C.04483 0.0531 0.059 235.75 72.96 0.81 C.04634 0.0298 0. 0 EPS2 C.07618 0.0211 0.030 184.59 104.65 0.51 M2 C.07769 0.0419 C.07900 0.0833 44.52 25.54 5.9 0.4904 0.035 4.51 19e+02 08202 0.0213 35.22 113.22 0.33 C.0833 0.059 0.043 149,12 45.60 ETA C.08507 0.0071 0.033 246.05 207.25 043 MO3 C.11924 0.0148 0.014 C.12077 0.0123 0.014 261.57 0.81 C.12229 0.012 331 144.92 0.18 SK3 C.1251 0.0023 0.010 237.69 219.86 0.054 M C,16102 0.0126 0.011 291.78 1,4 C.16384 0.0010 339.35 243 0.015 0.16667 0.0047 0.010 142.32 0.23 C.20280 0.005 310.10 0.067 2SK5 C.20845 0.0045 0,006 104.00 99.71 0.64 C.24002 0.0035 0.07 271.24 133.30 C,24153 0.0017 0.006 158.83 197.43 0.093 C.24436 0.54 2SM6 C,24718 0.007 298.92 175.13 0.11 212.25 2.1 C.32205 0.0030 0.004 75.29 0.55 M10 0.0009 0.003 0.089 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct 件 t example, mat内。这个数据集包含66天的海拔高度的小时值,其中有3天的数据是缺 失的。潮汐变化明显叠加在潮下的变化上面。时间序列可用演示程序τdemo进行加载和分 析。在脚本文件里对程序进行了描述。在这个例子里采用成份自动选择算法,选择了35个 成份。此外,一个浅水成份(M1o)被人为地添加进去,并对两个其他成份作了推理分析 f成份由K1进行推珥,K2由S2进行推理。进行了节点校。线性趋势没有包含在分析里。 彩色 bootstrap分析被用于决定显著性水平和置信区间。表I给出程序的输出结果。第一列 是给出成份的名字。显著成份(最后一列中SNR>1)用“*”标记。SNR是振幅(第3列) 和振幅误差(第4列)之比的平方。在重复分析中,山于 bootstrap程序的随机性,误差因 了(SNR随之)会发生很小的改变,但是振幅和相位本身不会变。频率(第1列)的单位 是cph,格林威治相位/相位误差(第5和第6列)的单位是“°”有1种成份(SNR截 止为2时,只有6个成份)被认为是显著的。图1B表示的是残差序列和由显著成分合成 的高程序列(注意,它将缺失的3天数据弥合)。分析的结果显小在图1C。有几个高频成 份勉强可以算为显著,但大多数显著成分是1日和半日带宽(大约分别为≈0,04和≈0.08 cph)。尽管大量能量是在2星期带宽(约0.002cph),拟合的成分却并不显著。显著成份的 相位示于图1D。显著成份一般有较小的相位误差 6.总结 在任何海洋时间序列的分析里,区分潮汐和非潮汐成份是一项重要工作。这里,我们 讨论经典调和分析的理论基础和所用 MATLAB程序包的使用细节。采用有关残差噪声结构 的3组不同假设中的一个,这个稈序包还能计算潮汐参数的置信区间。本文还用实例显示 了典型的结果。代码可访问htt:/ww. orgy. ubc.ca/^rich或IAMG服务器。 致谢 这个用于潮汐调和分析的 MATLAB工只箱的基本原理是近年来发展起来的,是作者们 对遡汐现象研究的一部分。我们感谢 Mikc forman的重要贡献,他在1977—1978年间向公 众提供了 FORTRAN代码和文件,成为我们的工作起点。 Richard Signer首先推导出本文所 用的线性化误差分析的公式, Julio candera和 Jim irish还提供了有用的信息。资助此项工作 的有加拿大自然科学和工程研究理事会的OGPO194270(RP)项目以及美国国家科学基金会 和海车研究署(SL和RB)。 *E: Classical tidal harmonic analysis including error estimates in MATLAB using T TIDE, Computers Geosciences 28(2002),929937 翻译:李进武 校对:张怀素 21994-2015ChinaacAdemicJOurnalElcctronicPublishingHousc.Allrightsrcscrved.http://www.cnki.nct 【实例截图】
【核心代码】