如何用调和常数推算潮汐

航海学（下）易错点总结7潮汐与潮流7．1潮汐7.1.1潮汐不等现象周日不等：在同一太阳日所发生的两次高潮或两次低潮的潮高以及相邻的高、低潮的时间间隔不相等。

成因：月赤纬≠0°且地理纬度ψ≠0°。

赤纬越大周日不等越明显。

分点潮无周日不等，回归潮周日不等最显著。

现象：一天一次高潮与一次低潮的条件ψ≥90°-Dec半月不等成因：月引潮力与太阳引潮力合力的变化；日、月与地球相互位置关系不同；月相不同。

现象：大潮和小潮潮汐半月变化规律：潮差的变化是以半个太阴月为周期（约14.5天）。

太阳的赤纬不等于0时，也会发生潮汐的周日不等现象。

视差不等：由地球和月球距离变化（注意：不是相对位置的变化）而产生的潮汐不等的现象。

周期：一个恒星月（约27.3天）太阳潮中也存在视差不等现象。

周期：一个回归年（约365.24日）简言之，视差不等是由于日、月、地三者空间距离的变化。

7.1.2潮汐类型半日潮型：一个太阴日内出现两次高潮和两次低潮，前一次高潮和低潮的潮差与后一次高潮和低潮的潮差大致相同，涨潮过程和落潮过程的时间也几乎相等（6小时12.5分）。

我国渤海、东海、黄海的多数地点为半日潮型，如大沽、青岛、厦门等。

全日潮型：一个太阴日内只有一次高潮和一次低潮。

如南海汕头、渤海秦皇岛等。

南海的北部湾是世界上典型的全日潮海区。

混合潮型：一月内有些日子出现两次高潮和两次低潮，但两次高潮和低潮的潮差相差较大，涨潮过程和落潮过程的时间也不等；而另一些日子则出现一次高潮和一次低潮。

我国南海多数地点属混合潮型。

如榆林港，十五天出现全日潮，其余日子为不规则的半日潮，潮差较大。

从各地的潮汐观测曲线可以看出，无论是涨、落潮时，还是潮高、潮差都呈现出周期性的变化，根据潮汐涨落的周期和潮差的情况，可以把潮汐大体分为如下的4种类型：正规半日潮：在一个太阴日(约24时50分)内，有两次高潮和两次低潮，从高潮到低潮和从低潮到高潮的潮差几乎相等，这类潮汐就叫做正规半日潮。

潮汐调和常数计算方式及其应用梁国亭李文学张晨霞（黄委员黄河水利科学研究院郑州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为地方迟角。

海洋要素计算作业之二——潮汐（威海2013年五月份）一．本次潮汐调和分析共选取了十三个分潮：MSf，Q1，O1，K1，P1，K2，N2，M2，S2，MK3，M4，MS4，M6为使您查看方便，将本次大作业的放在本文件夹各文件内，具体参考如下：1.原数据为：qd.dat;2.Fortran编程见该文件夹内：tide.f90文件；3.求各分潮调和常数H、g的值及其中间过程得到的各值见：qd_tide.dat文件；二．对比回报值和实测值：1. 回报1968年一月份的水位值见：huibao.dat；2. 用matlab绘制的潮汐过程曲线见：潮汐过程曲线.bmp3. 用给定的六个分潮求得的高潮和低潮发生的时刻及潮位值见—：gaodichao.dat;运行tide.f90后求得威海地区2013年5月份的平均潮差。

由图可知：由于只计算了一个月的潮汐数据，所以回报值和实测值相符的不是很好，如果计算一年的数据，应该会取得比较良好的结果。

三．程序%% 潮汐过程曲线图clear,clc%%huibao=load('G:\chaoxi\huibao.dat');% huibao=fread(fhuibao);shice=load('G:\chaoxi\qd.dat');% shice=fread(fshice);%huibao_y=zeros(1,12*62);%shice_y=zeros(1,12*62);huibao=double(huibao');huibao_y=double(huibao(:));%shice_y=reshape(shice',1,[])%for i=1:12;% for j=1:62% huibao_y(i)=huibao(i,j)% shice_y(i)=shice(i,j)%end%endshice=double(shice');shice_y=double(shice(:));x=linspace(1,31,length(huibao_y));plot(x,huibao_y,'r-')hold onplot(x,shice_y,'b-')title('威海(37°31′N ,122°08′E)2013年五月潮汐调和分析图') legend('回报值','实测值')xlabel('时间（2013年五月份）')ylabel('水位（m）')。

潮 汐计 算一 求任意时间的潮高和任意潮高的潮时（1）公式法1）求任意时间的潮高A ）以高潮为基准高、低潮时分别为T HW 、T LW ，高、低潮高分别为H HW 、H LW ，潮差R=H HW -H LW 。

θ为相位角。

△h 称为潮高改正数任意时潮高H t ： h H H HW t ∆-=而 )cos 1(2cos 22θθ-=-=∆RR R h︒⨯--=180HWLW HWT T T T θ所以，任意时潮高H t ： )cos 1(2θ--=RH H HW t )]180cos(1[2︒⨯----=HWLW HWHW T T T TRHT H T TH HT .D .BR △h△h′R （） /2/2R H 求任意时潮高高潮面任意时水面平均海面低潮面潮高基准面B ）以低潮为基准（对应图中的θ'与h '∆），任意时刻T 的潮高H t 可由下式求得：h H H LW t'∆+=而 )cos 1(2θ'-='∆R h 所以，任意时潮高H t ： )cos 1(2θ'--=R H H LW t )]180cos(1[2︒⨯---+=LWHW LW LW T T T T R H 例：求1992年2月8日铜沙T=1200的潮高H t 。

已知铜沙该日潮汐为0428 108；0959 418；1737 101；2219 350。

解：因为 6.471805909371759090012180︒=︒⨯--=︒⨯--=mh m h mh m h HW LW HW T T T T θ )(317101418cm H H R LW H W =-=-=)cos 1(2θ-=∆R h =)(5.51)6.47cos 1()2/317(cm =︒-⨯所以 )(36652418cm h H H H W t =-=∆-=以低潮为基准计算H t ：4.1321805909371700123717180︒=︒⨯--=︒⨯--=mh m h mh m h HW LW LW T T T T θ )cos 1(2θ'-=∆R h )(4.265)4.132cos 1(2317cm =︒-= )(366265101cm H H LW t =+==两种方法的计算结果完全一致。

!利用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。

江苏省燕尾港天文潮的推算戴荣，梁忠民河海大学水文水资源学院，南京（210098）E-mail ：dairong@摘 要：本文首先对燕尾港潮位站的实测资料进行潮汐调和分析，以求得潮汐各分潮的调和常数。

接着根据所得到的调和常数，进行潮汐推算。

计算结果表明该法具有精度高、速度快、简便实用等优点，很适合于燕尾港天文潮的推算。

关键词：潮汐调和分析，调和常数，潮汐推算，天文潮1 前言燕尾港位于灌云县的灌河口。

东临黄海，南濒灌河，与响水县相望。

由于来自太平洋通过东海，进入黄海的前进波在台风等气象因素的作用下，该海区经常受到风暴潮袭击，造成海堤损坏、海岸坍淤、部分地段的港槽摆动不定，不利于沿海地区的进一步开发，因而风暴潮灾害十分严重。

因为天文潮是风暴潮的主部，所以准确地预报出天文潮是提高风暴潮预报精度的关键之一。

虽说，国家每年都刊登下一年各潮位站的潮汐表，但由于人工将潮汐表中的数据输入到计算机中，耗费时间长，且易出现输入错误等缺点，因此本文利用燕尾港潮位站的资料，开发了自己的燕尾港天文潮推算程序。

2 天文潮的分析及预报计算天文潮预报方法通常可分为调和法和非调和法两类。

本文拟采用调和法进行江苏省燕尾港的天文潮预报。

调和法以潮汐静力学理论为基础，经过潮汐动力学理论的补充与完善。

根据调和法的计算过程，将其分为以下两步：前一步为调和分析，后一步为潮汐预报。

2.1 潮汐调和分析2.1.1 潮高公式海洋中的潮汐是在地球、月球和太阳三个天体相对运动中，由月球和太阳对地球的引潮力引起的[1]，将其按球面三角函数处理，并进一步用级数公式展开，可得出如下潮高公式： ])(cos[)(100∑=−+++=Q j j j j j j p g u t H fx t h νσ (2.1)式中，)(t h p 为t 时刻的潮位；0x 为平均潮位；H ,g 为调和常数；f ，u 为交点因子和订正角；σ为分潮的角速率；0v 为0=t 时刻的天文相角；Q 为分潮数。

潮汐调和分析的方法和应用研究‘.分类号学号:密级:玉??河海大哮硕士学位论文潮汐调和分析的方法和应用研究童章龙垂翅亟塑援.指导教师姓名??一盟渔太堂盔塞盛盗滥墨巫型王型国基重盛塞堕窒.申请学位级别专业名称堡堂亟±塑理洹注堂论文提交日期年月日论文答辩日期年月日学位授予单位塑塑盔堂学位授予日期 .生旦答辩委员会主席韭盈坌熬援论文评阅人毖盈丝塑堡选塑塑塾撞年月中国?南京摘要海岸附近和河口区域,是人类进行生产活动十分频繁的地带,往往也是人口最密集、经济最发达、开发程度最高的地区,而在这个地带潮汐现象显著,它直接或间接地影响着人们的生活。

研究潮汐,对人类的社会的发展有重要意义。

目前对潮汐的研究方法主要分为:潮汐调和分析和潮汐动力学两个方面。

本文主要对潮汐调和分析的一些理论和方法进行探讨。

本文分以下方面进行研究:第一、归纳潮汐的基本理论和基本概念,着重对使用一年资料进行调和分析的方法进行阐述,包括潮汐预报精度的评价标准。

第二、针对恶劣天气或仪器等意外情况导致潮位数据缺测或者具有重大误差等情况发生时,根据多次调和分析方法的思想,给出了基于连续函数最乘法的潮汐迭代调和分析方法,给出了方法收敛的条件。

与通常的多次调和分析法相比,该迭代方法不仅能够大大减少计算量,而且不用事先采取某种方法补全或替换原始资料。

然后将新建立的方法应用于多种实际情况中。

结果表明本方法是有效的。

第三、原始的天文相角是通过杜德森数表示的,而为了简化计算,现有的调和分析方法是通过角速度来表示天文相角的,省略了时间的二次项和三次项,这必然会引起误差。

本文直接用杜德森数表示天文相角,并用它直接进行调和分析,建立了基于杜德森数的调和分析方法,并用这种方法进行了实际的调和分析,并证明了通过角速度表示天文相角的可行性。

关键词:迭代法;调和分析:连续函数最小二乘法;多次分析:杜德森数., , ,,.,’’ .,:。

, . . : ., ,.,’, ,.,,,. .’ .,矾, , , .,. , .,.:; ;;学位论文独创性声明:本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。

一个简化的潮汐预报准调和分析方法王如云1，2,李慧娟1，2，蒋风芝2（1 河海大学水文水资源与水利工程科学国家重点实验室，南京210098；2 河海大学海洋学院，南京210098）摘要：在用现有的浅水港日潮汐准调和分析预报方法进行潮汐分析预报时，发现最小二乘法的法方程组的系数矩阵条件数很大，数量级在108，因此矩阵是坏条件的(或为病态的)，算法不稳定。

根据潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建立了一个简化的浅水准调和分析模型。

利用连云港的多年实测数据检验，简化的准调和分析模型相对于原准调和分析模型来讲，最小二乘法的法方程组系数矩阵条件数小很多，因此简化后的模型计算更为稳定。

在实测数据时间较长的情况下，简化前后的模型预报精度相当。

但当实测数据较短时，简化前的原模型却没有传统的调和分析模型的预报结果精度高，而简化后的模型却能保持比传统的调和分析模型的预报结果有一定的改善。

特别是简化后比简化前的模型计算时间减少了68%。

关键词：浅水潮汐；准调和分析；潮汐预报1引言在潮汐预报方面，一般采用调和分析方法，在深水区域此方法可以获得很好的预报效果，但在浅水区域尤其是河口区域，由于浅水潮汐的复杂性，采用此方法往往不能获得满意的效果。

例如杜德森提出的60个分潮[1]，其结果不能令人满意。

为此，杜德森后来又提出了一个直接对高低潮进行浅水改正的方法[2]，该方法虽然使高低潮的预报精度有了提高，但把它应用到逐时潮位预报上则有许多困难和不便之处。

在浅水区域由于非线性效应的加大，潮波往往产生畸变。

此时，高频振动的作用必须予以充分考虑。

为了提高浅水区域潮汐预报的精度，从调和分析方法来讲就必须增加高频的浅水分潮。

在水深不太浅的区域，浅水分潮的振幅会随着阶数的增高而迅速减小，所以在一般港口采用较少数目的主要浅水分潮即可满足潮汐预报的要求。

但在浅水区，常常需要考虑到六阶甚至更高阶的相互作用，才能满足潮汐预报的要求。