潮汐观测数据调和分析及预报成图

首先，需要说明两点。
第一，本案例使用的是单月潮汐观测数据，处理方法则是基于长期观测资料的调和分析来进行处理。中期观测资料的分析需要分别求主要分潮、随从分潮，短期观测资料分析则还需要计算不同观测序列的权重，但是核心算法与长期观测资料分析是一致的，都是建立矛盾方程，然后使用最小二乘法建立法方程，求出法方程系数，再求出矩阵X、Y。
第二，本文主要介绍方法步骤，所用代码大多为关键步骤实现，仅供参考。如需完整代码，请关注博主的另一篇资源。
下面开始介绍本案例的处理，从理论上讲，首先，我们需要选取分潮，确立我们所要分析的天文分潮，本案例用8个主要分潮——M2、S2、N2、K2、K1、O1、P1、Q1，四个半日潮，四个全日潮（其实去看潮汐相关研究的文献就会发现，基本都是以这八个分潮为主的）。相邻数据时间间隔为一小时，以所有数据的中间数对应时刻作为时间原点，然后对观测记录数据进行排序。
从实际出发，在matlab中，首先清理空间、准备环境（这是一个良好习惯），然后需要导入数据，对分潮进行排序。

%%导入数据
clear all;close all;clc;
data=importdata('C:\Users\STAR\Desktop\TideData_01.txt');%%分潮排序（八个）
M2=1;
S2=2;
N2=3;
K2=4;
K1=5;
O1=6;
P1=7;
Q1=8;

常量的准备，输入杜德森数，以及天文元素随时间的变化速度。

%%杜德森数（八个分潮各有七个杜德森数）
miu{M2}=[2,0,0,0,0,0,0];
miu{S2}=[2,2,-2,0,0,0,0];
miu{N2}=[2,-1,0,1,0,0,0];
miu{K2}=[2,2,0,0,0,0,0];
miu{K1}=[1,1,0,0,0,0,1];
miu{O1}=[1,-1,0,0,0,0,-1];
miu{P1}=[1,1,-2,0,0,0,-1];
miu{Q1}=[1,-2,0,1,0,0,-1];%%天文元素随时间的变化速度
rateOfChange=[14.49205211,0.54901653,0.04106864,0.00464183,0.00220641,0.00000196];

然后需要计算时间原点。个人习惯将代码模块化，本处使用自编函数TimeCalculation。

year=2003;
month=3;
day=1;
hour=0;
[year,month,day,hour]=TimeCalculation(2003,3,1,0,360);

计算分潮角速度

for i=M2:Q1sigma(i)=AngularVelocity(miu{i},rateOfChange);
end
sigma=deg2rad(sigma);   %注意，这里涉及到一个角度转弧度的操作

然后是计算各天文元素

[tao,ss,hhh,pp,NNN,ppp]=AstronomicalElements(year,month,day,hour);
tao=deg2rad(tao);
ss=deg2rad(ss);
hhh=deg2rad(hhh);
pp=deg2rad(pp);
NNN=deg2rad(NNN);
ppp=deg2rad(ppp);
astronomicalElements=[tao,ss,hhh,pp,NNN,ppp];

计算分潮初始位相，同样使用自编函数

for i=M2:Q1v0(i)=InitialPhase(miu{i},astronomicalElements);
end
v0(2) = 6.2832;%由于S2分潮的特殊性，直接赋值

计算交点因子及订正角

deltaMiu4{M2}=[0,0,0,2,2];
deltaMiu5{M2}=[-2,-1,0,0,1];
rho{M2}=[0.0005,-0.0373,1,0.0006,0.0002];
deltaMiu4{K2}=[0,0,0,0];
deltaMiu5{K2}=[-1,0,1,2];
rho{K2}=[-0.0128,1,0.2980,0.0324];
deltaMiu4{K1}=[-2,0,0,0,0,0];
deltaMiu5{K1}=[-1,-2,-1,0,1,2];
rho{K1}=[0.0002,0.0001,-0.0198,1,0.1356,-0.0029];
deltaMiu4{P1}=[0,0,0,2,2];
deltaMiu5{P1}=[-2,-1,0,0,1];
rho{P1}=[0.0008,-0.0112,1,-0.0015,-0.0003];
deltaMiu4{O1}=[0,0,0,2,2,2];
deltaMiu5{O1}=[-2,-1,0,-1,0,1];
rho{O1}=[-0.0058,0.1885,1,0.0002,-0.0064,-0.0010];
for i=[M2,K2,K1,P1,O1][f(i),u(i)]=IntersectionFactorAndCorrectionAngle(deltaMiu4{i},deltaMiu5{i},rho{i},pp,NNN);
end
f(Q1)=f(O1);
u(Q1)=u(O1);
f(N2)=f(M2);
u(N2)=u(M2);
f(S2)=1;
u(S2)=0;

建立法方程并求x和y

deltaT=1;
N=data(end,1);
NN=(N-1)/2;
A(0+1,0+1)=N;
for i=M2:Q1A(0+1,i+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2));A(i+1,0+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2));A(i+1,i+1)=(N+(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2;B(i,i)=(N-(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2;
end
for i=M2:Q1for j=M2:Q1if ~(i==j)A(i+1,j+1)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))+(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2 ;B(i,j)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))-(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2 ;end end
end
F1(0+1)=sum(data(:,2));
for i=M2:Q1F1(i+1)=sum(data(:,2).*cos((data(:,1)-361)*sigma(i)*deltaT));F2(i)=sum(data(:,2).*sin((data(:,1)-361)*sigma(i)*deltaT));
endX=F1/A(:,1:end);
Y=F2/B(:,1:end);

计算准调和振幅R和位相theta

R=(X(2:end).^2+Y.^2).^0.5;
for i=M2:Q1theta(i)=CalculatedPhase(R(i),X(i+1),Y(i));
end

计算分潮的调和常数

H=R./f;
g=theta+v0+u;
for i=M2:Q1g(i)=rad2pi(g(i));
end

潮汐预报、计算自报余差

S0=X(0+1);
for i=1:721h(i)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g));
end
r=data(1:end,2)'-h(1:end);
delta=sum(r.^2)^0.5/data(end,1);

计算预报潮位

starttime=(31-1+30+1)*24+1;
endtime=starttime+31*24+1;
forecastTime=1:31*24+1+1;
for i=starttime:endtimeforecastTide(i-starttime+1)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g));
end

绘制预报图

figure(2);
plot(forecastTime,forecastTide);
title("Forecast(May 1 to June 1)");
xlabel("Serial number");
ylabel("Stage");
legend("Forecast",'Location','Best');

写在文末，潮汐体现的是海洋动力及水文要素的变化规律和控制机制，从古至今，关于潮汐的研究从未停止，远有沈括著作《梦溪笔谈》，近有牛顿平衡潮理论之广泛应用，希望本文能为有需要的人提供一点思路及方法。

潮汐观测数据调和分析及预报成图相关推荐

自定义View之仿小米MIUI天气24小时预报折线图
本篇文章已授权微信公众号 hongyangAndroid(鸿洋)独家发布. 效果图本控件是仿MIUI8天气24小时预报折线图,用小米手机的可以打开天气软件看一下.本文是对自定义View的练手作品,要 ...
华测数据导出方法_输电线路塔基测量成图方法研究与应用
1.前言塔基测量成图软件是EATool第一个电力勘测模块,也是笔者个人刚参加工作时开发的第一个与工作直接相关的软件.该模块开发开始于2010年,于2012基本成型.记得第一次参加线路工程终勘定位,对 ...
在计算机中描述景物结构形状与外貌,在计算机中通过描述景物的结构、形状与外貌，然后将它绘制成图在屏幕上显示出来，此类图像称为_____。...
并实打印现共享文件和机等共享功能,计中通屏完全地相信平等互通,网工需的则所作模局域式是,网络的各没有台计一个内部主次之分若要算机. 算机述景包括条件主要区域规划发展. 过描的咨项目询包括(准备建设阶段 ...
无人机航测流程详解：航线规划、像控点布设、CC刺点建模及CASS成图
无人机航测是传统航空摄影测量手段的有力补充,具有机动灵活.高效快速.精细准确.作业成本低.适用范围广.生产周期短等特点,在小区域和飞行困难地区高分辨率影像快速获取方面具有明显优势,随着无人机与数码相机 ...
【终结扩散模型】Consistency Models.OpenAI开源新模型代码，一步成图，1秒18张
[终结扩散模型]Consistency Models.OpenAI开源新模型代码,一步成图,1秒18张 0.前言 Abstract 1. Introduction 2. Diffusion Model ...
Python 基础之词云（词的频率统计大小成图）的简单实现（包括图片词云，词云颜色，词的过滤）
Python 基础之词云(词的频率统计大小成图)的简单实现(包括图片词云,词云颜色,词的过滤) 目录
３０Ｍ　DEM数据处理成图
最近写论文要用到DEM数据,顺便整理了以下. 数据来源:NASA SRTM1 v3.0 由于数据是hgt文件类型,需要进行栅格数据转换步骤为: 打开ARCGIS 点击目录下Data manageme ...
全国大学生“高教杯“成图大赛：关于齿轮的快速建模研究
前言成图大赛全名为:"高教杯"全国大学生先进成图技术与产品信息建模创新大赛,是机械类全国大学生比赛制图类的最高等级比赛.目前已历经十四届.而机械类计算机建模的赛题,目前已相对固定 ...
全国大学生“高教杯”成图大赛——如何绘制锥度斜度
前言: 成图大赛全名为:"高教杯"全国大学生先进成图技术与产品信息建模创新大赛,是机械类全国大学生比赛制图类的最高等级比赛.目前已历经十四届.而机械类计算机建模的赛题,目前已相对固 ...

潮汐观测数据调和分析及预报成图

潮汐观测数据调和分析及预报成图相关推荐

最新文章

热门文章