[发明专利]一种潮汐预报方法

权利要求书

1.潮汐预测的方法，其特征主要包括如下步骤：

一、步骤(1)：安装SVM工具箱到MATLAB 7.8文件夹中

①将下载的SVM工具箱的文件夹放在/MATALB 7.8/toolbox/下；

②打开MATLAB-＞File-＞Set Path中添加SVM工具箱的文件夹；

③测试一下安装成功与否：在MATLAB中输入which svmoutput回车，如果可以正确显示路径，就证明添加成功了，例如：C:/Program Files/MATLAB71/toolbox/svm/svcoutput.m。

二、步骤(2)：利用MATLAB 7.8中SVM工具箱和历史潮汐数据构建潮汐预测网络模型，具体方法为：

①根据潮汐记录的时间(年、月、日、时、分、秒)，可计算

1)儒略世纪数：

其中Y表示年份， 的整数部分，D1为计算年首至计算日零时之间的天数；

2)儒略日：JD＝T×36525+2451545

3)月球平均纬度：s＝218.31643+481267.88128×T-0.00161×T²+0.000005×T³

4)太阳平均纬度：h＝280.46607+36000.7698×T+0.0003×T²

5)月球近地点平均纬度：p＝83.35345+4069.01388×T-0.01030×T²-0.00001×T³

6)太阳近地点平均纬度ps＝282.93835+1.71946×T+0.00046×T²+0.000003×T³

7)月球升交点平均纬度：

N＝125.04452-1934.13626×T+0.00207×T²+0.000002×T³

8)太阳偏心率：e_s＝0.01675104-0.0000418×T-0.000000126×T²

而地日中心的平均距离为一定值d_s＝1.49×10⁸km；

由地日中心的平均距离计算公式、太阳偏心率计算公式、太阳平均纬度计算公式和太阳近地点平均纬度计算公式，可计算地日中心的瞬时距离：

D_s＝d_s/(1+e_s×cos(h-ps)+e_s²×cos(2(h-ps)))

由太阳偏心率计算公式、太阳平均纬度计算公式和太阳近地点平均纬度计算公式，可以计算太阳的真经度：L_s＝h+2e_s×sin(h-ps)+1.25es²×sin(2(h-ps))

而黄道和天赤道的交角为：w_s＝23.452

再据潮位站点的地理北纬和东经(NLat，ELon)，可知由观测地点与格林尼治经度差产生的时间差订正值：±L＝(Long+M/60)/15

其中Long为观测点的经度，M为观测时刻的分钟值，东经取正号，西经取负号；记录的潮 汐观测时刻与格林尼治0时时间差订正值：W＝H+M/60，其中H为观测时刻的小时值；积日订正值：N₁＝(W-L)/24

根据计算年首至计算日零时之间的天数D₁，月球升交点平均纬度计算公式，积日订正值计算公式，可计算得到太阳赤纬的计算公式：

SDec＝0.3723+23.2567×sinθ+0.1149×sin2θ-0.1712×sin3θ

-0.758×cosθ+0.3656×cos2θ+0.0201×cos3θ

(θ＝2∏×57.2958×(D₁+N₁-N)/365.2422，N₂＝79.6764+0.2422(Y-1985))

由儒略日的计算公式，春分点的格林时角可由如下公式计算得到：

GHAγ＝100.075+0.985647348(JD-2400000.5-GMT/24-33282)

+15.04107×GMT

其中GMT为世界时间；

春分点格林时角与太阳平均纬度之差即得到太阳格林时角：GHAs＝GHAγ-h由此计算得到太阳地方时角的计算公式：

LHAs＝GHAs+λ(ELon)或LHAs＝GHAs-λ(WLon)

其中加代表潮位站是东经，减代表潮位站位于西经；

由太阳时角计算公式、太阳赤纬计算公式和验潮站的纬度，太阳高度可用如下公式计算得到：Hs＝arcsin(sin(NLat)×sin(SDec)+cos(NLat)×cos(SDec)×cos(LHAs))根据球面三角形原理，可得到太阳方位的计算公式：

由地月平均距离d_M＝3.84×10⁵km、月球轨道偏心率e＝0.0549、月球平均纬度、太阳平均纬度和月球近地点平均纬度计算公式，即可以计算第①个影响因子一地月瞬时距离计算公式：

D_m＝d_M/(1+e×cos(s-p)+(e²×cos(2(s-p))+0.14025e×cos(s-2h+p))

+0.00559504cos(2(s-h)))

也可计算月球的真经度：

L＝s+2e×sin(s-p)+1.25(e²×sin(2(s-p))+0.37249331e×sin(s-2h+p))

+(0.00975096+0.35064375e²)×sin(2(s-p))+0.01188946e×sin(s-2h-p)

+0.02692613e×sin(2s-h+ps) 

由黄道和白道的交角i＝5.15，可计算得到白赤交角的计算公式：

I＝arccos(cos(w_s)×cos(i)-sin(w_s)×sin(i)×cos(N))

春分点格林时角与月球平均纬度之差即得到月球的格林时角：GHAm＝GHAγ-s从而结合验潮站的经度计算得到月球的地方时角的计算公式：

LHAm＝GHAm+ELon

由白赤交角、月球真经度、月球时角，如下即为月球赤纬的计算公式：

从而月球高度可被计算为：

Hm＝arcsin(sin(NLat)×sin(MDec)+cos(NLat)×cos(MDec)×cos(LHAm))

月球方位也可随之得出计算公式：

②利用上述公式，在MATLAB 7.8平台编写程序计算影响潮位值的各天文因素，加上台风引起的增水值，即获取了支持向量机网络的输入变量，共11个，其对应的潮位值为网络的输出值，确定其输入输出量后，将其归一化到[-1，1]之间，其具体的归一化方法为：日地距离和月地距离采用最大最小值归一化方法，

Ds归一化的值＝(Ds-(max(Ds)+min(Ds))/2)/((max(Ds)-min(Ds))/2)；

Dm归一化的值＝(Dm-(max(Dm)+min(Dm))/2)/((max(Dm)-min(Dm))/2)；

太阳赤纬的变化范围为[-23.5，23.5]，月球赤纬的变化范围为[-29，29]，故其归一化方法为SDec归一化的值＝SDec/23.5；

MDec归一化的值＝MDec/29；

其他采用采用类似的方法进行归一化：

LHAs归一化的值＝LHAs/180；

LHAm归一化的值＝LHAm/180；

Hs归一化的值＝Hs/90；

Hm归一化的值＝Hm/90；

As归一化的值＝(As-180)/180；

Am归一化的值＝(Am-180)/180；

增水值归一化＝(增水值-min(增水值))/(max(增水值)-min(增水值))；

然后将获取的数据N组分成两组，N1组用作训练样本，N2组用作测试样本。

③确定支持向量机(SVM)的核函数为径向基函数(RBF)，其核函数的参数-惩罚参数C和函数的宽度参数g，采用网格搜索法求最优，网格搜索法即是C和g分别取M和N个值，对M×N个(C，g)个组合分别训练不同的SVM，然后再估计其学习精度，从M×N个组合中找到学习精度最小的值作为最优参数。

④支持向量机网络的结构图(见附图1) 

三、步骤(3)利用MATLAB 7.8工具箱svmtrain函数对网络进行训练，svmpridict函数对网络进行测试，具体方法为：

①用svmstrain训练的方法是，N1组训练数据和11个变量，组成N1×11归一化后的输入样本矩阵，对应的输出潮位矩阵为N1×1，设置径向基核函数的最优值后，运行svmtrain函数进行网络的训练，训练过程将以达到精度要求为停止条件，训练结束后将由系统自动生成一个SVM模型。

②用svmpredict函数对得到的模型进行测试，输入测试样本矩阵为N2×11。

③使用训练好和测试好的模型对同一潮位站的潮汐进行预测，选取需预报的时间M天，结合潮位站的地理位置，编程计算出月地距离，月球赤纬，月球时角，月球高度，月球方位，日地距离，太阳赤纬，太阳时角，太阳高度，太阳方位的值，即得到模型的输入矩阵M×11，运行svmpredict函数对这M天的潮汐值进行预测。

四、步骤(4)潮汐预测的具体算法如下：

①取得训练集T＝{(x₁，y₁)，(x₂，y₂)，…..，(x_N，y_N)}，N为训练集的样本数x_i为影响潮高的特征数据，即包含上述10个引潮力的相关参数和台风引起的增水值，y_i为对应时刻的潮高值。

②对SVM的模型参数进行初始化，将拉格朗日乘子α_i， 和阈值b赋以随机的初始值

③对训练数据构造相应的目标函数

其中K(x_i，x_j)为径向基核函数，求得最优解 

④选取α的一个正分量 它所对应的第j号样本即支持向量。

并据此计算阈值

⑤根据求得的 和b_j^*代入回归函数 用测试样本进行预测。

⑥计算误差函数，如果误差绝对值小于预先设定的ε，则学习结束，输出潮高的预测值，否则返回步骤③继续学习。

⑦画出算法流程图(见附图2)