【Python】海洋潮汐与水位控制 回归分析法计算验潮站长期平均海面

一、原理以及算法步骤

1.1回归分析法

  回归分析法，也称线性关系最小二乘法，基本原理是假定两站的短期距平具有比例关系。设比例为k,则有：

 根据之前所学的关于距平的公式：

将（1.2）代入（1.1）中，得：

上式中k未知，需要进一步假设两站的长期平均海面之间为线性比例关系，比例系数仍然为k，则有：

上式中C为位置常数项。将上式带入式（1.3）中，整理得到短期平均海面之间的关系如下 

 对比式（1.4）和（1.5）可知，此时两站的长期平均海面与短期平均海面具有相同的线性关系。

1.2算法步骤

1. 构造观测方程组

  计算两站每天的日平均海面，设为MSLA，MSLB

,其中i=1,2,...,n。式（1.5）是估计k与C的观测方程，将每天的日平均海面代入（1.5）得到观测方程：

 n天的日平均海面序列构建观测方程组如下：

 根据间接平差的原理，观测方程组为：

对照（1.7）与（1.8），式（1.8）中相应矩阵的具体形式如下：

2.组建法方程及求解

假设各观测（每天的如平均海面）互相独立，则观测值权阵可设为单位阵。按照间接平差原理，法方程为：

 由上式解得：

二、程序实现及实现步骤

2.1实现步骤

先由三年的A站，B站水位数据存储到一个excel文件中，使用python计算出年平均海面，再计算每天的所有日平均水位，计算十天、十五天和三十天的平均水位数据，代入上述方程进行解算。简易流程图如下：

 2.2程序实现及其结果

1.完整代码

import pandas as pd
import numpy as np
import xlwt
datadayavgA = []
hourdataA = []
sumA = 0
tendaysavgA = []
fifteendaysavgA = []
result1=[]
result2=[]
result3=[]
tendaysumA = 0
fifteendaysumA = 0
dataA = pd.read_excel(io="A站水位数据.xlsx", sheet_name="1", usecols=[1])#读取某一列的数值
df_liA = dataA.values.tolist()
for s_liA in df_liA:
      hourdataA.append(s_liA[0])

for i in range(0, 1096):

  for m in range(0 + 24 * i, 24 + 24 * i):

    if m == 23 + 24 * i:
            datadayavgA.append(sumA /24)
            sumA = 0
    else:
            sumA +=  hourdataA[m]

tendaysdataAdivide=[]
fifdaysdataAdivide=[]
for i in range(0,109):
    data=[]
    for m in range (0+10*i,10+10*i):
        if m==10+10*i:
            data.clear()
        else:
            data.append(datadayavgA[m])
    tendaysdataAdivide.append(data)
for i in range(0,73):
    data = []
    for m in range(0 + 15 * i, 15 + 15 * i):
        if m == 15 + 15 * i:
            data.clear()
        else:
            data.append(datadayavgA[m])
    fifdaysdataAdivide.append(data)
MSLLA=0
def calMSLL(data):
    sum=0
    for i in range (0,len(data)):
        sum+=data[i]
    return sum/len(data)
MSLLA=calMSLL(datadayavgA) #通过所有的数据计算出长期MSLA
datadayavgB = []
hourdataB = []
sumB = 0
tendaysavgB = []
fifteendaysavgB = []
tendaysumB = 0
fifteendaysumB = 0
dataB = pd.read_excel(io="B站水位数据.xlsx", sheet_name="1", usecols=[1])
df_liB = dataB.values.tolist()
for s_liB in df_liB:
      hourdataB.append(s_liB[0])

for i in range(0, 1096):

  for m in range(0 + 24 * i, 24 + 24 * i):

    if m == 23 + 24 * i:
            datadayavgB.append(sumB /24)
            sumB = 0
    else:
            sumB += hourdataB[m]
tendaysdataBdivide=[]
for i in range(0,109):
    data=[]
    for m in range (0+10*i,10+10*i):
        if m==10+10*i:
            data.clear()
        else:
            data.append(datadayavgB[m])
    tendaysdataBdivide.append(data)
fifdaysdataBdivide=[]
for i in range(0,73):
    data = []
    for m in range(0 + 15 * i, 15 + 15 * i):
        if m == 15 + 15 * i:
            data.clear()
        else:
            data.append(datadayavgB[m])
    fifdaysdataBdivide.append(data)
MSLLB=calMSLL(datadayavgB) #基准计算
tenresult = []
tendayssigma=[]
for i in range(0,109):
    one = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    B = np.mat([tendaysdataBdivide[i], one])
    L = np.mat(tendaysdataBdivide[i])
    sector1 = np.linalg.inv(B * B.T) * B * L.T#计算x
    res = sector1.tolist()
    tenresult.append(res[0])
    tenresult.append(res[1])
tengujiMSLB=[]
for i in range(0,109):
 calMSLB=tenresult[0+2*i][0]*MSLLA+tenresult[1+2*i][0]
 tengujiMSLB.append(calMSLB-MSLLB)
print("以十天为一时段回归分析的真误差","\n",tengujiMSLB)
#print("以十天为一时段回归分析的单位权中误差\n",tendayssigma)
fifdayssigma=[]
fifresult=[]
for i in range (0,73):
    one = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    B = np.mat([fifdaysdataAdivide[1], one])
    L = np.mat(fifdaysdataBdivide[1])
    sector1 = np.linalg.inv(B * B.T) * B * L.T
    res = sector1.tolist()
    fifresult.append(res[0])
    fifresult.append(res[1])
fifgujiMSLB=[]
for i in range(0,73):
 calMSLB=fifresult[0+2*i][0]*MSLLA+fifresult[1+2*i][0]
 fifgujiMSLB.append(calMSLB-MSLLB)
print("以十五天为一时段回归分析的真误差","\n",fifgujiMSLB)

thirdayssigma=[]
thirdaysdataBdivide=[]#三十天为一个时段

for i in range(0,36):
    data = []
    for m in range(0 + 30 * i, 30 + 30 * i):
        if m == 30 + 30 * i:
            data.clear()
        else:
            data.append(datadayavgB[m])
    thirdaysdataBdivide.append(data)

thirdaysdataAdivide=[]#三十天为一个时段
for i in range(0,36):
    data = []
    for m in range(0 + 30 * i, 30 + 30 * i):
        if m == 30 + 30 * i:
            data.clear()
        else:
            data.append(datadayavgA[m])
    thirdaysdataAdivide.append(data)

thirresult=[]
for i in range (0,36):
    one = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
    B = np.mat([thirdaysdataAdivide[1], one])
    L = np.mat(thirdaysdataBdivide[1])
    sector1 = np.linalg.inv(B * B.T) * B * L.T
    print(sector1)
    thirresult.append(res[0])
    thirresult.append(res[1])
thirgujiMSLB=[]

for i in range(0,36):
 calMSLB=thirresult[0+2*i][0]*MSLLA+thirresult[1+2*i][0]

 thirgujiMSLB.append(calMSLB-MSLLB)
#print("以三十天为一时段回归分析的真误差","\n",thirgujiMSLB)


xls=xlwt.Workbook()
sht1=xls.add_sheet('sheet')
sht1.write(0,0,label='十天真误差')
sht1.write(0,1,label='十五天真误差')
sht1.write(0,2,label='三十天真误差')
for i in range(len(tengujiMSLB)):
    sht1.write(i+1,0,label=tengujiMSLB[i])
for i in range(len(fifgujiMSLB)):
    sht1.write(i+1,1,label=fifgujiMSLB[i])
for i in range(len(thirgujiMSLB)):
    sht1.write(i+1,2,label=thirgujiMSLB[i])
sht1.write(0,3,label='MSLB的长期海平面')
sht1.write(1,3,label=MSLLB)
xls.save("结果分析1.xls")

2.程序结果截图

 三、结果分析

计算其与B站的长期平均水位的真误差，得到分析如下

图 3.1 真误差统计图

计算与B站的长期平均水位的相对误差，得到分析下图：

 图 3.2 相对误差统计分析图

3.2总结

从图3.1中可知，十天真误差的值约为32厘米，十五天真误差的值约为1厘米，三十天的真误差绝对值更小。由此可知，当时段为10天时，其误差较大，当时段为十五天甚至更多时，误差将达到1cm级别，可以认为其误差相对较小。

  从图3.2中可知，时段为十天的相对误差达到了百分之十以上，即可信度约为百分之九十以下，十五天、三十天的相对误差达到了百分之一甚至更小。

  由此可得出结论，当时段时间跨度越长时，其精度会更高，且十天回归分析的结果明显有较大差异，因此建议当进行回归分析法求解相邻站的长期海平面时应使用十天以上的观测时段进行计算，得到的精度更高。