下载文档原格式
潮汐调和分析的算法的报告，600字
潮汐调和分析（Tidal Harmonic Analysis）是一种分析水潮变
化的方法，用于对不同时间段内水潮高度变化进行分析。
潮汐调和分析技术可以帮助渔业、航海、灾害预测等行业理解水潮的演变。
潮汐调和分析的基本原理是，根据一定时间段内的水潮数据，从中提取出水潮变化的周期性变化。
通过统计学方法，把水潮变化的周期性变化表示为一系列的正弦函数相加的方式。
计算这些正弦函数的振幅、相位和频率就可以得到水潮变化的谐和非谐和分析结果。
潮汐调和分析一般采用Armadillo或Matlab来进行计算。
首先，将水潮数据转换为TXT文件格式，然后输入到Armadillo和Matlab软件中，使用潮汐调和分析的相关功能，对水潮数据
进行允许的处理。
主要的计算步骤是首先找到数据的频率，然后计算频率对应的振幅、相位和幅度分解值及其相关值，最后进行谐和非谐和分析，从而得出水潮变化的周期性特征。
潮汐调和分析有助于在一定时间段内对水潮变化情况进行分析，可以有效应用于渔业、航海和灾害预测等行业，有助于提高企业和行业的生产经营效率和绩效。
但是，由于潮汐调和分析的数据处理技术较为复杂，容易出现误差，因此在实际应用中应当备份数据、加强数据处理能力，以确保最终结果的准确性和可靠性。
潮汐调和常数计算方式及其应用梁国亭李文学张晨霞（黄委员黄河水利科学研究院郑州450003）（黄委会黄河水利技术学院475001）摘要在潮位预测中最大体的工作，第一就是计算潮汐调和常数。
本文在前人研究的基础上，开发了短时刻观测资料和30天观测资料的潮汐调和常数计算及潮汐预测模型，为深切研究黄河口的治理计划提供重要工具。
关键词潮汐分潮调和常数潮位预测1前言潮汐是河口最重要的海洋动力之一。
由于它周而复始的作用于河口、海岸的冲淤转变和入海泥沙等物质的扩散运移，直接影响着河口防洪安全、河口航运事业进展和人们的日常生活。
关于定点潮汐现象的预测研究，主如果按如实测资料通过调和分析，把复杂的潮汐曲线分解成许多调和项，即许多分潮，然后再按照调和常数和天文要素的转变推算潮汐。
所以调和常数是潮汐推算和进行潮波散布数值计算不可缺少的数据[1][2]。
河口泥沙数学模型除具有一般泥沙数学模型的特点外，还必需解决数学模型的下边界条件，即河口开边界的潮位进程线。
因此，研究潮汐调和常数的计算方式是研制河口泥沙数学模型的一项基础性研究工作。
2潮汐调和常数潮汐转变取决于地球、月球和太阳相对位置的转变。
按照万有引力定律，潮高的表达式，通过度解可取得月球平衡潮如下形式：ζ月=34(ME) (aD)3 a [(12-32Sin2ϕ) (23-2Sin2δ)+Sin2ϕSin2δCosT1+Cos2ϕCos2δCos2T1 ](1) 式中，M、E别离为月球和地球的质量，a～地球平均半径，D～地、月中心距，ϕ～地理纬度，δ～月球赤纬，T1～月球时角。
对(1)式中的变量，赤纬用经度、月球时角用太阳时替换，并引进辅助春分点，展开后略去4次方项，取得许多主要的调和项，即称分潮。
由于系数决定着潮差，相角决定着分潮周期。
在实际的海洋中，由于水流运动存在惯性、摩擦等缘故，天体在天顶时刻潮位并非发生最高，往往要掉队一段时刻才出现高潮，因此，通过对式(1)进行简化，可得：ζ=fHCos(σt+v0+u-K) (2) 式中的H为平均振幅，K为地方迟角。
潮汐简便计算法人们通过长期的实践、观察，发现海水有规律的涨落，而涨落的时间和高度又有着周期性的变化，由此人们把这种海水涨落的现象叫潮汐。
而随着海水的涨落、水位的升降，出现了海水的水平流动，这种海水流动的现象叫潮流。
海水有周期性涨落规律，如在每日里出现两次大潮和两次小潮。
通过长期实践、观察、发现每日的高潮大多出现在月亮的上、下中天（即过当地子午线时1前后。
低潮时间则在月出月落前后，并且每日的高（低）潮时间逐日后程约48分钟，即每天晚48分钟（0.8小时）。
每月的两次大潮是农历初一、十五附近几天，两次小潮是在农历的初七、八和甘二、廿三附近几天。
人们还发现，潮汐现象同月亮、太阳、地球的相对运动有密切的关系。
地球在一定轨道上绕太阳运转，月亮又在一定轨道上绕地球运转，它们之间有一定的吸引力和离心力，这种力就是产生潮汐现象的基本因素。
但实际潮汐涨落的主要成因却是月球对地球（表层）的吸引力，其次是太阳对地球的吸引力，太阳的乍用较小，约为月球的2/5，因月球离地球较近，故此月球的乍用较大。
据科学推测是：月球绕地球转，每一个月（29.5天多一点）转一圈，当月、日、地三者成一直线时，潮涨落的最大，这时是新月和望月（初一、十五）的时候，当日、月、地三者成直角三角形时潮涨落的最小，这是月上弦（初七、八）和下弦（廿二、廿三）的时候。
但在实际上形成大潮和小潮的时间，并不正好是上述时间，因为地球形状很复杂，所以各地发生最大潮和最小潮的时间要比理论上拖后几天。
如：山东半岛沿海每月的初三和十八潮的涨落最大，而初十和廿五前后潮的涨落又最小。
由于地球本身的自转，使地球上某点与月球的相对位置随时发生变化，这种变化每天（太阳约24时48分）为一周期。
每24时48分，发生两次高潮和两次低潮。
由高潮到低潮约经过6时12分，由第一个高潮到第二个高潮约经过12时24分。
潮汐的时间，在理论上应该与月球的上中天或下中天的时刻相符合，但实际上常常推迟。
!利用1996年7月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报program workimplicit nonecharacter*80::a1character(len=5),dimension(62,16)::aainteger::bb(62,12),c(62,2),caita(-371:371),i,i1,i2,j,t1real::N0,n(13,6),a(0:13,0:13),b(1:13,1:13),s,s0,s1,s2,s3,sa,hh !n代表Doodson代码；a,b为系数矩阵real::xiaoa(0:13),xiaob(13),gg1,gg2,pjchaocha,t,ma,mi!计算法方程所需的参数real,dimension(1:13)::w,u,f,V0,f1(0:13),f2!f1和f2为法方程右边系数real,dimension(13)::sita,h,g,r,h0(13),g0(13),h1(13),g1(13)!调和常数参数real,dimension(-371:371)::caita1,caita3,caita4,caita8,caita9,caita5,caita11 !主分潮、浅水分潮的潮高数值real,dimension(:),allocatable::hightide,lowtide,chaocha!高低潮数值integer,dimension(:),allocatable::hightrq,lowtrq,hight,hightt,lowt,lowtt!读取数据，把潮位数据赋值给bb，把年月份数据赋值给copen(unit=2,file='XM_July1996.dat')read(2,'(a)')a1print*,'数据文件的第一行信息:',a1do i=1,62read(2,'(16a5)')aa(i,:)end dodo i=1,62read(aa(i,5:16),*)bb(i,:)read(aa(i,3:4),*)c(i,:)end dodo i=1,62c(i,2)=int(real(c(i,2))/10.0)end doclose(2)!计算分潮角速率ww=(/0.002822,0.037219,0.038731,0.041781,0.163845,0.241534,0.078999 ,&&0.080511,0.083333,0.122292,0.161023,0.041553,0.083561/)w=360*wprint*print*,'角速率w:',w!计算N0 (middle time:1996-7-16 ; data sum:744, middle number:372 ) N0=259.157-19.32818*(1996-1900)-0.05295*(31*3+30*2+29+15+int((9 5.0)/4.0)) !初始升交点平均黄经N0=-(0.00220641*3+N0)print*!转换成格林威治时间print*,'N0:',N0!数字序号对应选取的分潮，但将5、6(P1、K2)分别与12、13(MS4、M6)对调，其中P1、K2为随从分潮!计算交点订正角uu(3)=10.8*sind(N0)-1.34*sind(2*N0)+0.19*sind(3*N0)u(4)=-8.86*sind(N0)+0.68*sind(2*N0)-0.07*sind(3*N0)u(8)=-2.14*sind(N0)u(13)=-17.74*sind(N0)+0.68*sind(2*N0)-0.04*sind(3*N0)u(1)=-u(8)u(2)=u(3)u(7)=u(8)u(9)=0u(10)=u(8)+u(4)u(11)=2*u(8)u(5)=u(8)u(6)=3*u(8)u(12)=0 !print*print*,'交点订正角u:',u!计算交点因子ff(3)=1.0089+0.1871*cosd(N0)-0.147*cosd(2*N0)+0.0014*cosd(3*N0) f(4)=1.006+0.115*cosd(N0)-0.0088*cosd(2*N0)+0.0006*cosd(3*N0) f(8)=1.0004-0.0373*cosd(N0)+0.0003*cosd(2*N0)f(13)=1.0241+0.2863*cosd(N0)+0.0083*cosd(2*N0)-0.0015*cosd(3*N0 )f(1)=f(8)f(2)=f(3)f(7)=f(8)f(9)=1f(10)=f(8)*f(4)f(11)=f(8)**2f(5)=f(8)**2f(6)=f(8)**3f(12)=1 !print*print*,'交点因子f:',f!查表得到的Doodson代码n(1,:)=(/0,2,-2,0,0,0/)n(2,:)=(/1,-2,0,1,0,0/)n(3,:)=(/1,-1,0,0,0,0/)n(4,:)=(/1,1,0,0,0,0/)n(5,:)=(/4,2,-2,0,0,0/)n(6,:)=(/6,0,0,0,0,0/)n(7,:)=(/2,-1,0,1,0,0/)n(8,:)=(/2,0,0,0,0,0/)n(9,:)=(/2,2,-2,0,0,0/)n(10,:)=(/3,1,0,0,0,0/)n(11,:)=(/4,0,0,0,0,0/)n(12,:)=(/1,1,-2,0,0,0/)n(13,:)=(/2,2,0,0,0,0/)!计算V0do i=1,13V0(i)=(14.49205212*3+180)*n(i,1)+(0.54901653*3+277.025+129.3848 *96+13.1764*(220))*n(i,2)+&&(0.04106864*3+280.190-0.23872*96+0.98565*(220))*n(i,3)+(0.00464 183*3+334.385+40.66249*96+&&0.11140*(220))*n(i,4)-(0.00220641*3+259.157-19.32818*96-0.05295* (220))*n(i,5)+&&(0.00000196*3+281.221+0.01718*96+0.000047*(220))*n(i,6) end doprint*print*,'初始幅角V0:',V0!设caita为潮高数据do i=1,61j=-371+(i-1)*12caita(j:j+11)=bb(i,:)end docaita(361:371)=bb(62,1:11)!计算法方程等式右边的数据，相邻数据时间间隔为1小时f1(0)=sum(caita)!f1为A阵中除第一行外的等式右边一维数据f1(1:13)=0do i=1,13do j=-371,371f1(i)=f1(i)+caita(j)*cosd(j*w(i))end doend do!f2为B阵中等式右边的一维数据f2=0do i=1,13do j=-371,371f2(i)=f2(i)+caita(j)*sind(j*w(i))end doend do!计算A阵中的系数矩阵Aa(0,0)=743do j=1,13a(0,j)=sind(743.0/2*w(j))/sind(0.5*w(j))a(j,0)=a(0,j)end dodo j=1,13a(j,j)=0.5*(743+sind(743.0*w(j))/sind(w(j)))end dodo i=1,13do j=i+1,13a(i,j)=0.5*(sind(743.0/2*(w(i)-w(j)))/sind(0.5*(w(i)-w(j)))+& & sind(743.0/2*(w(i)+w(j)))/sind(0.5*(w(i)+w(j)))) a(j,i)=a(i,j)end doend doprint*print*,'系数矩阵A:',a!计算B阵中的系数矩阵Bdo j=1,13b(j,j)=0.5*(743-sind(743.0*w(j))/sind(w(j)))end dodo i=1,13do j=i+1,13b(i,j)=0.5*((sind(743.0/2*(w(i)-w(j)))/sind(0.5*(w(i)-w(j)))-&& sind(743.0/2*(w(i)+w(j)))/sind(0.5*(w(i)+w(j))))) b(j,i)=b(i,j)end doend doprint*print*,'系数矩阵B:',b!Guass-Seidel迭代法求解方程组h=0;g=0;i1=0doh0=hg0=g!A阵do i=0,11s1=0do j=0,11s1=s1+xiaoa(j)*a(i,j)end doxiaoa(i)=-s1/a(i,i)+f1(i)/a(i,i)+xiaoa(i)-xiaoa(12)*a(i,12)/a(i,i)-xiaoa(13)* a(i,13)/a(i,i)end do!B阵do i=1,11s1=0do j=1,11s1=s1+xiaob(j)*b(i,j)end doxiaob(i)=-s1/b(i,i)+f2(i)/b(i,i)+xiaob(i)-xiaob(12)*b(i,12)/b(i,i)-xiaob(13) *b(i,13)/b(i,i)end do!计算调和常数h,gdo j=1,11sita(j)=atand(xiaob(j)/xiaoa(j))+180r(j)=sqrt(xiaoa(j)**2+xiaob(j)**2)g(j)=V0(j)+u(j)+sita(j)h(j)=r(j)/f(j)end dodo i=1,11do while(g(i)>360.or.g(i)<0)if(g(i)>360)thendog(i)=g(i)-360if(g(i)>0.and.g(i)<360)exitend doelseend ifif(g(i)<0)thendog(i)=g(i)+360if(g(i)>0.and.g(i)<360)exitend doelseend ifend doend dogg1=g(4)-g(3)gg2=g(9)-g(8)do while(gg1>230.or.gg1<-130)if(gg1>230)thendogg1=gg1-360if(gg1<230.and.gg1>-130)exitend doelseend ifif(gg1<-130)thendogg1=gg1+360if(gg1>-130.and.gg1<230)exitend doelseend ifend dodo while(gg2>230.or.gg2<-130)if(gg2>230)thendogg2=gg2-360if(gg2<230.and.gg2>-130)exitend doelseend ifif(gg2<-130)thendogg2=gg2+360if(gg2>-130.and.gg2<230)exitend doelseend ifend dog(12)=g(4)-0.075*(g(4)-g(3))g(13)=g(9)+0.081*(g(9)-g(8)) sita(12)=-(u(12)+V0(12)-g(12)) sita(13)=-(u(13)+V0(13)-g(13)) h(12)=h(4)*0.324h(13)=h(9)*0.282do j=12,13r(j)=h(j)*f(j)xiaoa(j)=r(j)*cosd(sita(j))xiaob(j)=r(j)*sind(sita(j)) end dog1=g-g0h1=h-h0i1=i1+1if(all(abs(h1)<10.0).and.all(abs(g1)<2.0))exit !退出循环条件end doprint*print*,'系数a:',xiaoaprint*print*,'系数b:',xiaobprint*print*,'迭代循环次数:',i1print*print*,'调和常数h:',hprint*print*,'调和常数g:',g!计算平均水位s0=xiaoa(0)print*print*,'平均水位s0:',s0!向文件中输入数据!将各分潮的调和常数写入'hg.txt'open(unit=2,file='hg.txt')do i=1,13write(2,*)h(i),g(i)end doclose(2)!将所有潮汐数据写入'tides.txt'open(unit=2,file='tides.txt')do i=1,62write(2,'(12i5)')bb(i,:)end doclose(2)!将所有主要分潮(3,4,8,9)、浅水分潮(5,11)数据按随时间的变化情况写入向量中do j=-371,371caita3(j)=s0+f(3)*h(3)*cosd(w(3)*j-sita(3))caita4(j)=s0+f(4)*h(4)*cosd(w(4)*j-sita(4))caita8(j)=s0+f(8)*h(8)*cosd(w(8)*j-sita(8))caita9(j)=s0+f(9)*h(9)*cosd(w(9)*j-sita(9))caita5(j)=s0+f(5)*h(5)*cosd(w(5)*j-sita(5))caita11(j)=s0+f(11)*h(11)*cosd(w(11)*j-sita(11))caita1(j)=caita3(j)+caita4(j)+caita8(j)+caita9(j)+caita5(j)+caita11(j)-5*s0 end do!计算高低潮个数i1=0;i2=0do j=-370,370if(caita(j)>caita(j-1).and.caita(j)>caita(j+1))theni1=i1+1 !高潮个数end ifif(caita(j)<caita(j-1).and.caita(j)<caita(j+1))theni2=i2+1 !低潮个数end ifend doprint*print*,'高潮个数:',i1print*,'低潮个数:',i2if(allocated(hightide))deallocate(hightide)if(allocated(lowtide))deallocate(lowtide)if(allocated(hight))deallocate(hight)allocate(hightide(1:i1),lowtide(1:i2))!将观测时的高潮写入文件open(unit=2,file='realh.txt')do i=1,size(hightide)write(2,*)hightide(i)end doclose(2)!将观测时的低潮写入文件open(unit=2,file='reall.txt')do i=1,size(lowtide)write(2,*)lowtide(i)end doclose(2)!潮汐预报部分!计算高低潮对应的潮位及时刻if(allocated(hightide))deallocate(hightide) if(allocated(lowtide))deallocate(lowtide) if(allocated(hight))deallocate(hight)if(allocated(lowtrq))deallocate(lowtrq)if(allocated(hightrq))deallocate(hightrq)i=min(i1,i2)allocate(hightide(1:i1),lowtide(1:i2),hightrq(i1),lowtrq(i2),hight(i1),hightt (i1),lowtt(i2),lowt(i2),chaocha(i))i1=0;i2=0do j=-370,370!高潮潮位及时刻if(caita1(j)>caita1(j-1).and.caita1(j)>caita1(j+1))theni1=i1+1do t=j-1,j+1,0.01s=s0+f(3)*h(3)*cosd(w(3)*(t)+V0(3)+u(3)-g(3))+&f(4)*h(4)*cosd(w(4)*(t)+V0(4)+u(4)-g(4))+&f(5)*h(5)*cosd(w(5)*(t)+V0(5)+u(5)-g(5))+&f(8)*h(8)*cosd(w(8)*(t)+V0(8)+u(8)-g(8))+&f(9)*h(9)*cosd(w(9)*(t)+V0(9)+u(9)-g(9))+&f(11)*h(11)*cosd(w(11)*(t)+V0(11)+u(11)-g(11)) sa=s0+f(3)*h(3)*cosd(w(3)*((t+1/60.0))+V0(3)+u(3)-g(3))+& f(4)*h(4)*cosd(w(4)*((t+1/60.0))+V0(4)+u(4)-g(4))+&f(5)*h(5)*cosd(w(5)*((t+1/60.0))+V0(5)+u(5)-g(5))+&f(8)*h(8)*cosd(w(8)*((t+1/60.0))+V0(8)+u(8)-g(8))+&f(9)*h(9)*cosd(w(9)*((t+1/60.0))+V0(9)+u(9)-g(9))+&f(11)*h(11)*cosd(w(11)*((t+1/60.0))+V0(11)+u(11)-g(11)) if(s<sa)thenma=sahightt(i1)=int((t-floor(t))*60)hight(i1)=floor(t)+371end ifend dohightide(i1)=maelseend if!低潮潮位及时刻if(caita1(j)<caita1(j-1).and.caita1(j)<caita1(j+1))theni2=i2+1do t=j-1,j+1,0.01s=s0+f(3)*h(3)*cosd(w(3)*(t)+V0(3)+u(3)-g(3))+&f(4)*h(4)*cosd(w(4)*(t)+V0(4)+u(4)-g(4))+&f(5)*h(5)*cosd(w(5)*(t)+V0(5)+u(5)-g(5))+&f(8)*h(8)*cosd(w(8)*(t)+V0(8)+u(8)-g(8))+&f(9)*h(9)*cosd(w(9)*(t)+V0(9)+u(9)-g(9))+&f(11)*h(11)*cosd(w(11)*(t)+V0(11)+u(11)-g(11)) sa=s0+f(3)*h(3)*cosd(w(3)*((t+1/60.0))+V0(3)+u(3)-g(3))+& f(4)*h(4)*cosd(w(4)*((t+1/60.0))+V0(4)+u(4)-g(4))+&f(5)*h(5)*cosd(w(5)*((t+1/60.0))+V0(5)+u(5)-g(5))+&f(8)*h(8)*cosd(w(8)*((t+1/60.0))+V0(8)+u(8)-g(8))+&f(9)*h(9)*cosd(w(9)*((t+1/60.0))+V0(9)+u(9)-g(9))+&f(11)*h(11)*cosd(w(11)*((t+1/60.0))+V0(11)+u(11)-g(11)) if(s>sa)thenmi=salowtt(i2)=int((t-floor(t))*60)lowt(i2)=floor(t)+371end ifend dolowtide(i2)=mielseend ifend do!将高潮位写入'hightide.txt'，第1列为潮位，第2、3列为时刻open(unit=2,file='hightide.txt')do i=1,i1write(2,*)hightide(i),hight(i),hightt(i)end doclose(2)do i=1,i1j=1if(hight(i)<=23)thenhightrq(i)=1end ifdo while(hight(i)>23)hight(i)=hight(i)-24j=j+1hightrq(i)=jend doend doopen(unit=2,file='hightidexiu.xls')do i=1,i1hightide(i)=int(hightide(i)+0.5)*1.0write(2,*)hightide(i),hightrq(i),hight(i),hightt(i) !将高潮位及时刻写入'hightidexiu.xls'，第1列为潮位，第2、3、4列为天数、小时、分钟end doclose(2)!将低潮位写入'lowtide.txt'，第1列为潮位，第2、3列为时刻open(unit=2,file='lowtide.txt')do i=1,i2write(2,*)lowtide(i),lowt(i),lowtt(i)end doclose(2)do i=1,i2j=1if(lowt(i)<=23)thenlowtrq(i)=1endifdo while(lowt(i)>23)lowt(i)=lowt(i)-24j=j+1lowtrq(i)=jend doend doopen(unit=2,file='lowtidexiu.xls')do i=1,i2write(2,*)lowtide(i),lowtrq(i),lowt(i),lowtt(i) !将低潮位及时刻写入'lowtidexiu.xls'，第1列为潮位，第2、3、4列为天数、小时、分钟end do!将四大主分潮潮位数据写入'zhuyaotides.txt'中，第1列为O1，第2列为K1，第3列为M2，第4列为S2open(unit=2,file='zhuyaotides.txt')do i=-371,371write(2,'(4f7.1)')caita3(i),caita4(i),caita8(i),caita9(i)end doclose(2)!将浅水分潮潮位数据写入'qianshuitides.txt'中,第1列为MS4，第2列为M4open(unit=2,file='qianshuitides.txt')do i=-371,371write(2,'(2f7.1)')caita5(i),caita11(i)end doclose(2)!计算平均潮差s=0do i=1,size(chaocha)s=s+(hightide(i)-lowtide(i))pjchaocha=s/size(chaocha)print*print*,'平均潮差:',pjchaochaprint*deallocate(hightide,lowtide,hight,lowt,chaocha,hightrq,lowtrq) !释放掉潮位数据!判断XM站潮汐类型hh=(h(3)+h(4))/h(8)if(hh<0.5)thenprint*,'潮型数：',hh,'所以XM站潮汐类型为半日潮'else if(hh>=0.5.and.hh<2.0)thenprint*,'潮型数：',hh,'所以XM站潮汐类型为不规则半日潮混合潮' else if(hh>=2.0.and.hh<=4.0)thenprint*,'潮型数：',hh,'所以XM站潮汐类型为不规则日潮混合潮' else if(hh>=4.0)thenprint*,'潮型数：',hh,'所以XM站潮汐类型为全日潮'elseprint*,'error'end ifend program work。
潮汐调和分析的方法和应用研究‘.分类号学号:密级:玉??河海大哮硕士学位论文潮汐调和分析的方法和应用研究童章龙垂翅亟塑援.指导教师姓名??一盟渔太堂盔塞盛盗滥墨巫型王型国基重盛塞堕窒.申请学位级别专业名称堡堂亟±塑理洹注堂论文提交日期年月日论文答辩日期年月日学位授予单位塑塑盔堂学位授予日期 .生旦答辩委员会主席韭盈坌熬援论文评阅人毖盈丝塑堡选塑塑塾撞年月中国?南京摘要海岸附近和河口区域,是人类进行生产活动十分频繁的地带,往往也是人口最密集、经济最发达、开发程度最高的地区,而在这个地带潮汐现象显著,它直接或间接地影响着人们的生活。
研究潮汐,对人类的社会的发展有重要意义。
目前对潮汐的研究方法主要分为:潮汐调和分析和潮汐动力学两个方面。
本文主要对潮汐调和分析的一些理论和方法进行探讨。
本文分以下方面进行研究:第一、归纳潮汐的基本理论和基本概念,着重对使用一年资料进行调和分析的方法进行阐述,包括潮汐预报精度的评价标准。
第二、针对恶劣天气或仪器等意外情况导致潮位数据缺测或者具有重大误差等情况发生时,根据多次调和分析方法的思想,给出了基于连续函数最乘法的潮汐迭代调和分析方法,给出了方法收敛的条件。
与通常的多次调和分析法相比,该迭代方法不仅能够大大减少计算量,而且不用事先采取某种方法补全或替换原始资料。
然后将新建立的方法应用于多种实际情况中。
结果表明本方法是有效的。
第三、原始的天文相角是通过杜德森数表示的,而为了简化计算,现有的调和分析方法是通过角速度来表示天文相角的,省略了时间的二次项和三次项,这必然会引起误差。
本文直接用杜德森数表示天文相角,并用它直接进行调和分析,建立了基于杜德森数的调和分析方法,并用这种方法进行了实际的调和分析,并证明了通过角速度表示天文相角的可行性。
关键词:迭代法;调和分析:连续函数最小二乘法;多次分析:杜德森数., , ,,.,’’ .,:。
, . . : ., ,.,’, ,.,,,. .’ .,矾, , , .,. , .,.:; ;;学位论文独创性声明:本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。
中期水位观测资料的最小二乘分析报告摘要：本次实验中采用了KM站（28.05N,121.17E）1997年8月的逐时潮位数据，运用中期水位观测资料的最小二乘分析方法，通过奇异值剔除、调和分析、逐时潮位回报、高低潮及余水位计算等工作，对此验潮站的数据进行了初步分析，并了解了中期水位资料分析的基本流程。
报告主要内容：（1）11个主要分潮（MSf Q1 O1 K1 N2 M2 S2 MK3 M4 MS4）及两个随从分潮（P1 K2）的调和常数H和g（2）图像和数据文件的基本信息（3）平均潮差和潮汐类型（4）余水位特征分析（5）误差分析（6）程序的相关说明（1）潮汐的调和常数：利用最小二乘原理，通过引入差比关系的方法，我们可以成功得到11个主要分潮和2个随从分潮的调和常数如下：分潮名称调和常数H 调和常数gMSf 121.2222 -32.38747Q1 62.95736 233.5120O1 225.5294 238.7111K1 266.1612 113.2537N2 420.5689 167.2492M2 1922.772 174.8581S2 679.3940 197.3759MK3 33.19594 252.1002M4 32.43390 121.7806MS4 33.60584 198.8826M6 3.762754 94.29744P1 73.46050 109.5160K2 192.9479 201.4156程序运行结果如图：其中H关系到分潮的振幅，g关系到分潮的相位。
从表中可以看出，M2分潮的振幅最强，对当地潮位的贡献最大，这与实际情况相符，但K1分潮的调和常数H仅有266.1412，结果偏小。
（2）图像和数据文件的基本信息：本次报告中包含以下数据文件：1. KM9708new.dat数据原始文件。
2. KM9708new_02.dat经过奇异值订正的数据文件，为方便画图时读取，没有输入数据质量信息。
怎么分析潮汐趋势的方法潮汐是地球上因引力而形成的周期性海水涨落现象。
随着时间的推移，潮汐会出现趋势性变化，例如海平面逐渐升高或降低。
这些趋势可能与气候变化、海洋和大气动力学过程以及地球物理过程等因素有关。
因此，分析潮汐趋势具有重要的科学和管理意义。
以下是几种常用的方法。
一、线性回归分析法线性回归分析法是一种常用的趋势分析方法，可以用来研究潮汐现象的长期趋势。
具体来说，可以使用时间序列数据拟合一条直线，根据回归方程的斜率确定趋势的方向和大小。
如果斜率为正，说明趋势向上增加；如果斜率为负，说明趋势向下减少。
线性回归分析要求样本量充足、数据稳定，并且假设回归残差服从正态分布。
由于潮汐随时间的变化通常具有季节性变化，因此可以将每年的同一月份的潮汐数据放在一起进行分析。
二、滑动窗口分析法滑动窗口分析法是一种非常有效的方法，可以用来确定潮汐短期趋势的变化。
该方法将时间序列数据分为连续的多个时间片段，在每个时间片段内计算平均值或其他统计量。
通过比较相邻时间片段的平均值，可以确定趋势的方向和大小。
如果平均值在相邻时间段内增加，则趋势向上增加；如果平均值在相邻时间段内减少，则趋势向下减少。
滑动窗口分析法的优点在于可以捕捉到潮汐短期趋势的变化，而不受长期趋势的影响。
三、小波分析法小波分析法是一种时频分析方法，可以用来研究潮汐的时间和频率特征。
该方法将时间序列数据经过小波变换，将其变换到时频域中。
在时频域中，可以看到不同时间尺度和频率的潮汐成分。
通过对不同成分的分析，可以确定潮汐趋势的变化。
例如，在低频成分中，可以确定潮汐的长期趋势；在高频成分中，可以确定潮汐的短期波动。
四、离群点检测法离群点检测法是一种异常检测方法，可以用来识别异常或伪异常的潮汐观测数据。
在潮汐数据中，可能存在因仪器误差、天气因素及人为操作等因素引起的异常或伪异常现象。
这些异常或伪异常数据可能会对潮汐趋势的分析和预测造成干扰。
因此，离群点检测方法可以帮助识别和过滤掉这些异常数据，从而提高潮汐趋势的准确性。
一个简化的潮汐预报准调和分析方法王如云1，2,李慧娟1，2，蒋风芝2（1 河海大学水文水资源与水利工程科学国家重点实验室，南京210098；2 河海大学海洋学院，南京210098）摘要：在用现有的浅水港日潮汐准调和分析预报方法进行潮汐分析预报时，发现最小二乘法的法方程组的系数矩阵条件数很大，数量级在108，因此矩阵是坏条件的(或为病态的)，算法不稳定。
根据潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建立了一个简化的浅水准调和分析模型。
利用连云港的多年实测数据检验，简化的准调和分析模型相对于原准调和分析模型来讲，最小二乘法的法方程组系数矩阵条件数小很多，因此简化后的模型计算更为稳定。
在实测数据时间较长的情况下，简化前后的模型预报精度相当。
但当实测数据较短时，简化前的原模型却没有传统的调和分析模型的预报结果精度高，而简化后的模型却能保持比传统的调和分析模型的预报结果有一定的改善。
特别是简化后比简化前的模型计算时间减少了68%。
关键词：浅水潮汐；准调和分析；潮汐预报1引言在潮汐预报方面，一般采用调和分析方法，在深水区域此方法可以获得很好的预报效果，但在浅水区域尤其是河口区域，由于浅水潮汐的复杂性，采用此方法往往不能获得满意的效果。
例如杜德森提出的60个分潮[1]，其结果不能令人满意。
为此，杜德森后来又提出了一个直接对高低潮进行浅水改正的方法[2]，该方法虽然使高低潮的预报精度有了提高，但把它应用到逐时潮位预报上则有许多困难和不便之处。
在浅水区域由于非线性效应的加大，潮波往往产生畸变。
此时，高频振动的作用必须予以充分考虑。
为了提高浅水区域潮汐预报的精度，从调和分析方法来讲就必须增加高频的浅水分潮。
在水深不太浅的区域，浅水分潮的振幅会随着阶数的增高而迅速减小，所以在一般港口采用较少数目的主要浅水分潮即可满足潮汐预报的要求。
但在浅水区，常常需要考虑到六阶甚至更高阶的相互作用，才能满足潮汐预报的要求。
上个世纪六十年代，一些潮汐学者试图通过扩充高频分潮的数目以使预报结果获得改进，如Zelter and Cumimngs[3]以及Rossiter and Lennon[4]曾将分潮的数目扩充到110多个，但效果并不理想。
方国洪等人认为，不理想的原因在于随着频率的增加，高频分潮的数目极速的增加，不可能从中挑选出少数分潮近似代替所有分潮，难以用有限数目的浅水分潮来体现总的浅水效应。
可以认为通过增加浅水分潮以改进潮汐预报，其效果可能是比较有限的。
基于以上分析，方国洪等提出了一个浅水潮汐预报的准调和分析方法[5]，可以用来推算任意时刻的潮高，也可以用来推算高、低潮，效果比传统的调和分析法有了显著的改进。
但我们使用此方法对连云港的多年潮位实测数据进行分析预报时，发现最小二乘法的法方程组系数矩阵条件数很大，算法不稳定。
为此，我们对浅水准调和分析模型进行了简化，简化后的模型计算更为稳定，计算时间大为减少。
2 准调和分析方法介绍方国洪等人[5]提出的浅水预报准调和方法思路是把潮高分做两部分，一部分为低频部分，由基金项目：水文水资源与水利工程科学国家重点实验室开放研究基金(2005407411)；中国教育部科学技术研究重点项目（104104）；中国江苏省普通高等学校高新技术产业发展项目(JH03-010)作者简介：王如云(1963-)，教授，男，安徽芜湖人，从事计算物理学研究，E-mail:wangry@潮族0，1，2组成，主要是天文源潮波，另一部分为高频部分，属于浅水分潮。
分析分两步进行，首先对于低频部分，即潮族0，1，2，用如下表1列出的分潮计算。
表1中包含了Doodson [5]所用的所有属于潮族0，1，2的分潮。
对实测潮汐()t ξ进行调和分析，求出各分潮的调和常数。
然后用实测水位减去平均水位和0，1，2族的潮位，依据剩余值再作进一步分析。
表1 调和分潮假如利用一年潮位资料()t ξ，对其进行调和分析，计算各分潮的潮汐调和常数，再依据下式将长周期、全日、半日潮族分别加以组合，得出三个基本准调和分潮，然后把三分日及以上的高频潮表示为这些基本准调和分潮的函数。
随着频率的增加，只增加少量的准调和项，每一个准调和项可以近似看作一群频率相近的分潮之和。
计算每小时0，1，2族的振幅()R t 和()r t 。
各族的振幅和位相的变化已不再是常数，而具有缓慢的变化。
称这些量是准调和分量。
5001cos cos[()]i i i i i i R r f H t v u g σ==++-∑5001sin sin[()]i i i i i i R r f H t v u g σ==++-∑32116cos cos[()]i i i i i i R r f H t v u g σ==++-∑ （1）32116sin cos[()]i i i i i i R r f H t v u g σ==++-∑582233cos cos[()]i i i i i i R r f H t v u g σ==++-∑582233sin cos[()]i i i i i i R r f H t v u g σ==++-∑从动力学的原因来看，浅水分潮由两种非线性产生。
一种是运动方程中的平流项，如u ux∂∂，...，和连续方程中的非线性项，如()u xς∂∂，...，能够产生高级摄动项。
另一种是由运动方程中的摩擦项所产生，如u u ρ（这里u 是流速，ς是水位，ρ是摩擦系数）。
将上式定义的准调和项视为单一潮波，并用cos()kA a θ-表示由非线性产生的属于高频部分的各阶摄动项，则潮高可表达为：cos()j j j j h k A a ξθ=+-∑ （2） 式中5801cos[()]i i i i i i h A f H t v u g σ==+++-∑ （3）式(2)中cos()KA a θ-为浅水分潮部分，用34个准调和项表示浅水效应，式中,K θ为准调和常数，,A a 为已知参数，它们的表达式见表2。
预报时，首先由表1中第1到第58个分潮的调和常数按(3)式计算h 。
再由(1)式计算各族的R 和r ，根据R 和r 按表2计算,A a 。
利用文[1]中的潮汐预报准调和分析方法对连云港多年的实测潮位资料进行了分析预报，得出的结论是改善效果明显，但浅水系数矩阵的条件数极大（见表4），这意味着方程组病态很严重，该算法不稳定。
在进行浅水准调和分析时，选择的浅水分潮并非越多越好，选择的分潮达到一定个数后，浅水系数矩阵奇异严重，预报精度不仅没有改善，反而会由于矩阵奇异给预报可靠度带来负面作用。
实际上在观测记录数据资料的过程中，会有恶劣天气或仪器磨损等意外情况导致潮位数据缺测或者具有重大误差等情况发生及观测数据长度等的不同，算法若不稳定，这些原始误差就会在预报过程中被扩大，导致预报结果的可靠性下降，这对预报是极为不利的。
鉴于这些，我们对原准调和分析方法进行了简化。
根据潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建立了一个简化的浅水准调和分析模型。
表2 浅水准调和项中,A a 表达式jj A j ajj Aj a 12 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1721R 21R R 221R R 21R R 221R R 321R R 31R 21R R 221R R 231R R 22R 32R 220R R 221R R 21R R 221R R 321R R13r 21r r + 21r r + 122r r - 122r r -122r r - 14r 212r r + 212r r + 212r r + 22r 22r 22r 213r r + 212r r +212r r +212r r + 1819 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34221R R 321R R 221R R 22R 32R 42R 221R R 321R R321R R 421R R 32R 42R 421R R 42R 52R 52R 62R123r r - 123r r - 2122r r + 23r 23r 23r 213r r + 213r r + 124r r - 124r r - 24r 24r 214r r + 25r 25r 26r 26r3 简化的准调和分析模型在忽略了摩檫力作用的前提下，研究潮波在一维等深等宽半无限长沟渠中传播的情况，简化后的方程组00u u u g t x xu u h u tx x x ζζζζ∂∂∂⎧++=⎪⎪∂∂∂⎨∂∂∂∂⎪+++=⎪∂∂∂∂⎩中，第一个方程的第二项，第二个方程的第三和第四项与相应方程的主要项相比，其量级即使在浅水地区也小于1，作为零级近似，把它们略去了。
因此倍潮波的振幅与相应源潮波振幅的整数次方成正比例，且幂次与频率的倍数相等；复合潮波的振幅则正比例于源潮波振幅的乘积，这个结果对由无摩檫引起的线性浅水分潮大致上是成立的，但是受摩檫力非线性作用引起的浅水潮波，它们的振幅与源潮波的振幅则不一定遵从以上关系，实际的振幅变化幅度要小于理论的变幅[6]。
因为摩檫力比例于潮汐振幅的平方，而不是线性关系[1]，不能认为引潮力增大若干倍，潮汐振幅也相应增大相同的倍数。
由此我们假设浅水分潮的一般表达式为cos()K a θ-（表达式中a 的含义同原来的模型），同时由于简化后的计算模型避免了对表2中浅水准调和项表达式A 的计算，可节省计算时间。
简化模型求解浅水调和常数时方法同原来的相似。
最后进行预报并与原模型进行分析比较。
对于表2中a相同的项，我们只须保留一项即可，否则浅水系数矩阵的行列式由于存在完全相同的两行或多行，其行列式为零，从而系数矩阵奇异，方程组无解。
简化后的准调和项为17项，a的表达式见表3。
表4为简化前后模型的预报结果比较，仍以1967年实测潮位作为分析资料，来预报其它年份水位。
为书写方便，我们称简化前模型为模型1，简化后模型为模型2。
表3 准调和项中a的表达式表4 模型1和模型2均方差比较表年份调和分析均方差（㎝）模型1均方差（㎝）模型2均方差（㎝）1962 30.5 28.6 29.51963 24.2 21.3 22.61964 24.5 21.7 231965 24.8 22.1 23.31966 24.6 21.7 231967 20.8 17.3 18.91968 24.2 21.6 22.91969 28.3 26.1 27.11970 26 23.5 24.71971 24.8 22.1 23.51972 27.6 25.3 26.41973 27.1 24.8 25.91974 28.5 26.3 27.31975 26.8 24.6 25.51982 24.4 21.2 22.6浅水系数矩阵条件数模型1：35912244 模型2：14分析表4，预报精度虽没有所提高，但两者均方差相差仅为0.9cm到1.6cm之间，同时模型2的浅水系数矩阵条件数相比原来的算法显然小很多，说明简化后的算法稳定。
