潮汐调和分析方法的探讨-张凤烨

68海洋科学/2011年/第35卷/第6期潮汐调和分析方法的探讨张凤烨1,2,魏泽勋1,2,王新怡1,2,王永刚1,2,方国洪1,2(1.海洋环境科学与数值模拟国家海洋局重点实验室,山东青岛266061;2.国家海洋局第一海洋研究所,山东青岛266061)摘要:为实现对不等时距潮汐资料的分析,基于Matlab内部函数功能,提出了一种调和分析方法。基于这种方法,分别对大连、北海两个站位1985年的全年等时间间距取样的资料和非等时间间距取样的资料进行了调和分析,结果显示,由等时间距资料和非等时间距资料计算的调和常数基本吻合。对大连、北海两个站位的全年资料进行多个不同时间间距取样分析,发现当分潮频率大于取样频率的二分之一时,分潮发生频率混淆。若分潮周期明显大于样品长度,该分潮的分析结果产生很大误差。最后得出结论为:此调和分析方法,适合对非等时间间距、非连续潮汐潮流资料进行调和分析,并且能够获得与传统方法精度相当的结果。关键词:Matlab;潮汐;调和分析;取样中图分类号:P731.23文献标识码:A文章编号:1000-3096(2011)06-0068-08在实际的海洋调查中,尤其是长时间的观测,由于受到仪器故障、恶劣天气、地理位置制约,观测方式等因素的影响,很难得到从观测初始时刻到结束时刻这段时间内完整的高质量数据资料,也可以说获得的海洋资料是不完整的、不连续的。对于潮汐资料,调和分析方法作为一种主要的分析方法,发展至今,理论基础已经很成熟,但是传统的调和分析大多是对观测时间间隔为等距的资料进行调和分析,这是由于非等时序列数据的法方程系数矩阵很大,所需的计算量也非常大,通常要求观测数据等时间间隔以减少计算量,而且常用的观测资料也是以观测时间间隔为等距离为主。若资料中出现数据缺测漏测,或者因为异常天气和其他原因使得部分实测数据相对正常潮汐存在大的误差,就首先要对数据进行前期处理[1],然后应用最小二乘法对潮汐进行调和分析,这种方法在计算时,首先将原始矛盾方程进行求导处理,从而得到矛盾方程的法方程,这样大大增加了计算的难度[2]。如果我们从线性代数这门学科的角度来考虑,矛盾方程组只是一个线性代数中的超定方程组,只需要直接解这个超定方程组即可,基于Matlab内部函数功能求解线性方程组能很好地解决这一问题,并且方程组所得解即为最小二乘解[3-4],从而简化了调和分析的算法。有多少样本,就有多少个方程,更能体现最小二乘法的灵活性,完全不受等时或非等时取样的限制。本文采用基于Matlab内部函数功能的最小二乘法实现潮汐资料的调和分析,并应用于大连、北海两站位的水位观测资料分析,通过对两站位全年资料的调和分析,从非理论推导的角度,分析基于Matlab内部函数功能的调和分析方法是否可应用于等时间间隔和不等时间间隔的观测资料的调和分析,以及观测资料的取样间隔对于调和分析精度的影响,然后通过对两站位7月份资料的调和分析,来分析采样长度对潮汐调和分析精度的影响。1改进的潮汐调和分析方法介绍调和分析是建立在分潮的概念上[5],将潮汐看成是以不同频率传播的各种潮波叠加产生的现象,调和分析的目的就是求出各个分潮的振幅和迟角,迟角是指:某时刻、某一地点实际的分潮的相角与理论上该时刻的分潮相角的差值,迟角分为3种:有区时迟角、地方迟角、格林威治迟角。现在国际上通用的是格林威治迟角,用字母g表示,本文所用的迟角即为格林威治迟角[6]。收稿日期:2011-01-24;修回日期:2011-04-21基金项目:国家自然科学基金项目(40976016);国家海洋局第一海洋研究所基本科研业务费专项资金项目(GY02-2010G07);国家高技术研究发展计划(863计划)重点项目(2009AA121402)作者简介:张凤烨(1985-),男,河北唐山人,硕士研究生,主要从事物理海洋学研究,E-mail:[e-mail];魏泽勋,通信作者,研究员,E-mail:[e-mail]/Vol.35,No.6/201169潮位表示式如下:00Greenwichicos[()]iiiiiAfHtVugςω=+++−(1)式中A0为平均水位高度,i为分潮序列号,ζ表示分潮潮高,ω为分潮的角速度,f,u分别表示由于月球轨道18.61a的周期变化引进的对平均振幅H和相角的订正值,即交点因子和交点订正角,V0为格林威治初相角,H和g为分潮的调和常数。本文利用潮位公式,建立一种基于Matlab内部函数功能的潮汐调和分析方法,其区别于传统的调和分析方法,其特殊之处在于,这种基于Matlab内部函数功能的调和分析方法可以对不等时间间隔的潮汐、潮流资料进行调和分析,而且计算易于实现。下面将介绍这种基于Matlab内部函数功能的潮汐调和分析的方法和原理。设有观测资料:ζ(j),j=1,2,3…N,N代表数据个数,该资料每∆t小时为一次记录,此处的∆t既可以是整数,也可以是小数或者是一个不定值。取数据的开始时刻为时间原点。设潮汐为M个分潮叠加而成:00101cos()(cossin)MkkkkMkkkkkARtAAtBtωθωω===+−=++∑∑ς(2)其中00cos,sinkkkkkkARBRθθ==。由于每个分潮的ω是一个定值,所以每个时刻的costω,sintω是一个常数,这里A0是平均水位,我们可以把它看做是ω为0的一个分潮。这样就形成了一个有2M+1个未知量的方程个数为N的线性方程组,然后应用基于Matlab的内部函数功能的最小二乘法,直接对线性方程组进行求解,从而得最小二乘解,即每个对应分潮的A,B。每个分潮的ω根据下式求得:123456shpNpωµτµµµµµ′′′=+++++(3)其中:123456,,,,,µµµµµµ为杜德森数,,,',,,'shpNpτ为只与时间有关的基本天文要素。通过Matlab可以求出一个平均水位、N/2个A值和N/2个B值。根据22kkkRab=+,k=θarctg()jjba得出各分潮的振幅和相位。最后只要依据天文要素计算出初始时刻的f,u和V0,此处f,u和V0为格林威治时刻值。根据kkkHRf=,0kkkkgVuθ=++,即可求出调和常数,其中f,u的计算方法参照引潮力的第四展开式[1,8]。2站点水位的潮汐调和分析下面将采用基于Matlab内部函数功能的潮汐调和分析方法对大连、北海两个站点等间隔资料、非等间隔资料以及不同时距的等时间隔进行调和分析,讨论本方法的适用性以及适用条件。2.11a资料的调和分析利用基于Matlab内部函数功能的调和分析方法对两个站点的1a的水位资料进行调和分析,从潮族中选取具有代表意义的分潮,选取了30个分潮,分别为K1,O1,Q1,P1,MP1,J1、SO1,OO1,M1,M2,S2,N2,K2,OQ2,L2,O3,MO3,K3,MK3,M4,MS4,Mk4,SN4,2MK5,MSK5,M6,MSk6,2Mk6,Sa,Ssa。分别对大连、北海1a(1985年)的观测资料进行调和分析,这里原始观测数据取样的时间间隔是1h,而对于非等时的观测资料的分析,首先对原始数据进行处理,去掉3月份的所有的观测资料,然后随机删除一些水位观测资料,使其达到一种非等时距采样的效果,然后对其进行调和分析,从图1可以看到,等时调和分析预报的结果和不等时调和分析预报的结果与原始的观测资料基本上是重合的,从表1、表2可以看出,等时间序列调和分析得到的调和常数与非等时间序列的调和常数的值非常接近,表中提到的原始值代表的是利用传统的调和分析方法对潮汐观测数据进行调和分析所得的调和常数,由于本文中提到的原始调和常数是对26a数据进行调和分析由135个分潮所得,所以一些分潮的调和常数与本文计算的调和常数有一定的差别,如因为M1分潮的杜德森数自身的问题,所以其调和常数不予考虑,但是主要分潮的调和常数与本文计算所得分潮调和常数相当接近。这说明基于Matlab内部函数功能的调和分析方法适用于潮汐资料的调和分析。下面的实验主要是探讨数据的取样间隔对该方法所得调和常数的影响,为了从结果上讨论取样频率是否对基于Matlab内部函数功能的调和分析方法有影响,本文分别对大连、北海两站位时长为1a、取样时间间隔分别为1,2,3,4,5,6h的观测数据,应用基于Matlab内部函数功能的调和分析方法进行调和分析,然后对调和常数进行对比,结果见表3、表4、表5、表6。70海洋科学/2011年/第35卷/第6期图11985年1月大连、北海潮位观测数据、等时和非等时后报值图Fig.1Thechartoftidevalues,showingtheobservationalandthehindcastvaluesofharmonicanalysisofdatasampledatequaltimeintervalsandnon-equaltimeintervalsinDalianandBeihaiinJanuary1985表1大连原始、等时与不等时资料调和常数比较Tab.1Thecomparisonofharmonicanalysisofdatasampledatequaltimeintervalsandnon-equaltimeintervalsinDalianH(cm)g(°)分潮原始值(传统方法)等时(本方法)不等时(本方法)原始值(传统方法)等时(本方法)不等时(本方法)M296.8095.7595.89287.38287.0286.97S230.4330.6030.67341.04339.94340.13K125.6025.3725.561.520.310.0O118.1418.4418.61321.55322.63322.21M40.810.710.69111.80113.36117.49MS40.250.290.28192.43168.37163.09Q13.513.903.93296.14300.21298.98P17.497.537.35354.07356.02355.82N218.2218.9318.38257.53259.80259.06K28.438.708.62339.29336.43335.92M61.040.991.0581.8881.8983.21Sa24.6826.9826.40126.77125.47122.89Ssa2.641.523.13339.31290.97314.14MP10.570.540.68169.34132.83128.79J11.401.121.1742.5344.5344.81SO10.930.840.84185.03208.54202.78OO10.620.570.5980.0994.9498.07M10.890.690.67342.84153.35147.22OQ20.540.560.34232.19248.80259.63L24.977.907.72326.83323.44326.45O30.090.040.0570.49175.0232.96MO31.191.261.22111.04117.12116.43K30.070.190.21219.1060.0262.83MK31.491.561.61157.54154.06155.58Mk40.080.100.13240.84258.56267.60SN40.040.110.1551.5686.0776.042MK50.150.100.0929.1745.1817.48MSK50.110.120.1663.2536.06343.81MSk60.190.220.20187.23204.36190.502Mk60.390.360.33134.59137.73131.87MarineSciences/Vol.35,No.6/201171表2北海原始、等时与不等时资料调和常数比较Tab.2Thecomparisonofharmonicanalysisofdatasampledateq