本文以GOCE卫星1b数据为例，介绍GOCE重力梯度数据的读取与坐标转换（GRF—LNOF）。

数据介绍

数据选取2013年3月1日至3月31日共计31天、分辨率为1秒的观测数据，数据观测记录2678400个。一个观测记录包含以下31项数据，数据格式如下：

（1）：GPS时间/s；

（2-3）：纬度、经度/°；

（4）：地心半径/m；

（5-10）：XX、YY、ZZ、XY、XZ、YZ六个方向分量上的组合梯度值；

（11-16）：六个方向分量梯度值的观测中误差；

（17-22）：六个方向分量的标识码；

（23-31）：GRF框架转换为LNOF坐标系的旋转矩阵9个元素。
数据读取

fid=fopen('GGC_GRF_DataFile_GOCO03s_20130301T000000_20130331T235959_0002.DAT');
data=zeros(2678400,31) %定义数据矩阵

%读取数据
for n=1:2678400
    line =fgetl(fid);
    if ~ischar(line), break, end
  if isempty(line), continue, end
  cells = textscan(line,'%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%d%d%d%d%d%d%f%f%f%f%f%f%f%f%f');
  for k=1:31  
      data(n,k)=cells{k};
  end
end
fclose(fid);


坐标转换

GOCE卫星重力梯度观测数据是在重力梯度仪参考框架（GRF）下的，为了后续计算，需要将其转换成惯性坐标系（IRF）、地固坐标系（EFRF）或者局部指北坐标系（LNOF）。此处以转换成LNOF为例。

%声明一个存储LNOF系下的梯度矩阵
 Vlnof=zeros(2678400,6);
 
 %计算局部指北坐标系（LNOF）下的GOCE重力梯度Vlnof 
  for  n=1:2678400
     Vgrf=[data(n,5) data(n,8) data(n,9);data(n,8) data(n,6) data(n,10);data(n,9) data(n,10) data(n,7)];
     R=[data(n,23) data(n,24) data(n,25);data(n,26) data(n,27) data(n,28);data(n,29) data(n,30) data(n,31)];
     Vnew=R*Vgrf*R';
     Vlnof(n,1)=Vnew(1,1);Virf(n,2)=Vnew(2,2);Virf(n,3)=Vnew(3,3);Virf(n,4)=Vnew(1,2);Virf(n,5)=Vnew(1,3);Virf(n,6)=Vnew(2,3); 
  end

李宏毅2020机器学习作业——Linear Regression
一、问题描述
​	1. 数据描述
​		本次作业采用的数据是使用台湾某一座城市的观测记录，分别为train.csv,test.csv，其中train.csv
共有4320行, 27列，包含了该地点18种化合物一年12个月，每月钱20天，每天24小时的数值监控，test.csv则是从剩下资料当中取出的连续10小时为一组，前九个小时的观测值当做特征，第十小时的pm2.5当做answer,一共取出240组不重复的数据

2. 问题要求

​		预测test.csv中240组测试集中的PM2.5
二、数据处理


导入必要的库


import pandas as pd
import numpy as np
import math
data = pd.read_csv('train.csv')

读取对应文件


数据处理的必要性


首先观测我们需要训练的feature与answer在表格中的关系



很容易发现1月1日的0-8小时18种化合物的所有数据都可以作为训练的数据，而第9小时的pm2.5的30数值作为这第一个训练数据的answer，从上表我们可以看出需要处理有两个部分，

那个NR值需要处理为数值型，转化为0 2.
数据的格式需要进行处理，每月每天每日的数据不连续，不易于后续处理



数据预处理
#保留完全是数据的行列并导出为numpy
data = data.iloc[:, 3:]
data[data == 'NR'] = 0
raw_data = data.to_numpy()
print(raw_data)

output
[['14' '14' '14' ... '15' '15' '15']
 ['1.8' '1.8' '1.8' ... '1.8' '1.8' '1.8']
 ['0.51' '0.41' '0.39' ... '0.35' '0.36' '0.32']
 ...
 ['36' '55' '72' ... '118' '100' '105']
 ['1.9' '2.4' '1.9' ... '1.5' '2' '2']
 ['0.7' '0.8' '1.8' ... '1.6' '1.8' '2']]

转化成numpy更有利于后期对test与train数据集进行一些转化的骚操作
下面才是真正的骚操作
将原始的数据按照月份重组，因为每个月只取了20天，20天又是连续的，那么每个月就会有(20×24)×18组数据集可供训练，每年就是12×(20*24)×18，而通过观测表格中的数据，你会发现每天24小时后数据在下一天的表格里，而16-24的数据明明可以以下一天的第0小时的pm2.5作为标签训练的！
所以就要把原来train的4320×24的数据按照每月重组成12个18×480的数据，就是把一个月里，第二天的数据放到第一天的相邻的旁边，效果图如下



#对data进行调整，将4320*24重组为12*18*480
month_data = {}
for month in range(12):
    sample = np.empty([18,480])
    for day in range(20):
        sample[:,day * 24 : ( day + 1 ) * 24] = raw_data [ 18 * ( 20 * month + day ) : 18 * ( 20 * month + day + 1 ),: ]
    month_data[month] = sample

这一波操作秀吗？建立一个字典month_data存储每个月的数据，那么现在思考一下train_data与label的关系了，train_data就是九个小时的18种化合物的观测数据，如图1的绿色框框的内容，lable即红色框框的内容，而整个数据集中month_data[i]的列数即为每个月的小时数24×20，行数当然是18(化合物数量)，那么其中的train_data组数即为480-9组，也就以为每个月有471个样本可以训练，其中每个样本的形状就是绿色框框显示的18×9，每年就有12×471个样本可以训练，将一整年的train_data放入一个数组x中，那么x的维度即为(12×471,18×9)，那么将一整年的label放入y数组中，其维度为(12×471,1)
x = np.empty([12*471,18*9],dtype = float)
y = np.empty([12*471,1],dtype = float)
for month in range(12):
    for day in range(20):
        for hour in range(24):
            #一个月的最后一天的数据不能与下个月相连
            if day == 19 and hour>14:
                continue
            x[month * 471 + day * 24 + hour, :] = month_data[month][:,day * 24 + hour : day * 24 + hour + 9].reshape(1,-1)
            y[month * 471 + day * 24 + hour, 0] = month_data[month][9,day * 24 + hour + 9]
print(x)
print(y)

output
[[14.  14.  14.  ...  2.   2.   0.5]
 [14.  14.  13.  ...  2.   0.5  0.3]
 [14.  13.  12.  ...  0.5  0.3  0.8]
 ...
 [17.  18.  19.  ...  1.1  1.4  1.3]
 [18.  19.  18.  ...  1.4  1.3  1.6]
 [19.  18.  17.  ...  1.3  1.6  1.8]]
[[30.]
 [41.]
 [44.]
 ...
 [17.]
 [24.]
 [29.]]

4. 数据处理
​	将数据标准化

$x[i][j] = \frac{x[i][j] - mean_x[j]} {std_x[j]}$

code
#标准化
mean_x = np.mean(x,axis = 0)
std_x = np.std(x,axis = 0)
for i in range(len(x)):
    for j in range(len(x[0])):
        if std_x[j] != 0:
            x[i][j] = (x[i][j] - mean_x[j]) / std_x[j]
#将训练集分成训练-验证集，用来最后检验我们的模型
x_train_set = x[: math.floor(len(x) * 0.8), :]
y_train_set = y[: math.floor(len(y) * 0.8), :]
x_validation = x[math.floor(len(x) * 0.8): , :]
y_validation = y[math.floor(len(y) * 0.8): , :]
print(x_train_set)
print(y_train_set)
print(x_validation)
print(y_validation)
print(len(x_train_set))
print(len(y_train_set))
print(len(x_validation))
print(len(y_validation))

output
[[-1.35825331 -1.35883937 -1.359222   ...  0.26650729  0.2656797
  -1.14082131]
 [-1.35825331 -1.35883937 -1.51819928 ...  0.26650729 -1.13963133
  -1.32832904]
 [-1.35825331 -1.51789368 -1.67717656 ... -1.13923451 -1.32700613
  -0.85955971]
 ...
 [ 0.86929969  0.70886668  0.38952809 ...  1.39110073  0.2656797
  -0.39079039]
 [ 0.71018876  0.39075806  0.07157353 ...  0.26650729 -0.39013211
  -0.39079039]
 [ 0.3919669   0.07264944  0.07157353 ... -0.38950555 -0.39013211
  -0.85955971]]
[[30.]
 [41.]
 [44.]
 ...
 [ 7.]
 [ 5.]
 [14.]]
[[ 0.07374504  0.07264944  0.07157353 ... -0.38950555 -0.85856912
  -0.57829812]
 [ 0.07374504  0.07264944  0.23055081 ... -0.85808615 -0.57750692
   0.54674825]
 [ 0.07374504  0.23170375  0.23055081 ... -0.57693779  0.54674191
  -0.1095288 ]
 ...
 [-0.88092053 -0.72262212 -0.56433559 ... -0.57693779 -0.29644471
  -0.39079039]
 [-0.7218096  -0.56356781 -0.72331287 ... -0.29578943 -0.39013211
  -0.1095288 ]
 [-0.56269867 -0.72262212 -0.88229015 ... -0.38950555 -0.10906991
   0.07797893]]
[[13.]
 [24.]
 [22.]
 ...
 [17.]
 [24.]
 [29.]]
4521
4521
1131
1131

三、数据处理(training)
​	因为本次作业只是简单的回归模型，所以对于已有的训练数据集x，满足公式$x*w+b=y$即可求出预测值y，而w的参数就是我们需要求解的值，我们将x的列数扩展一个维度作为参数b，可以得知参数w的维度为(18×9+1,1),给w初始化后，如何优化参数，那就是老生常谈的内容——梯度下降，但本次作业对学习率的调节方式值得学习。


损失函数

$L=\sqrt{\frac{\sum_{i=1}^{n}(xw-y)^2}{n}}$


w的梯度值

$dw = \frac{x(xw-y)}{nL}$


adagrad
adagrad的参数更新方式与传统的梯度下降方式的参数更新方式不同，如你所知，传统的参数更新方式为：

$w^{t+1}:=w^t - a*dw^t$

再观看了李宏毅老师的gradient decent_1的视频课后，发现对于含有多个参数的模型优化，这时如果所有参数还采用同一学习率乘以梯度作为参数移动的步长，这时参数可以移动的步长取决于参数的梯度，而在每个参数距离这个参数的最优解之间的距离并不取决于它的梯度，而是它的梯度除上它的二阶导，如下所示：

$w^{t+1}:=w^{t}-\frac{a}{\sqrt{\sum_{i=0}^{t}(dw^i)^2}}dw^t$

上式的a均为学习率，若是想深剖细节可去学习李宏毅机器学习的gradient decent_1的视频课，其中的上述式子的均方根的值可近似为L对w的二阶导的数值，将上述公式应用与本题，可得以下公式：

$w^{t+1}:=w^{t}-\frac{a}{\sqrt{\sum_{i=0}^{t}(dw^i)^2+b}}dw^t$

其中b为调节参数
#因为存在偏差bias，所以dim+1
dim = 18 * 9 + 1
# w维度为163*1
w = np.zeros([dim,1])
# x_train_set维度为 4521*163
x_train_set= np.concatenate((np.ones([len(x_train_set),1]),x_train_set),axis = 1).astype(float)
#设置学习率
learning_rate = 10
#设置迭代数
iter_time = 30000
#让dw值初始化
adagrad = np.zeros([dim,1])
eps = 0.0000000001
for t in range(iter_time):
    loss = np.sqrt(np.sum(np.power(np.dot(x_train_set,w)-y_train_set,2))/len(x_train_set))
    if(t%100 == 0):
        print("迭代的次数：%i ， 损失值：%f"%(t,loss))
        gradient = (np.dot(x_train_set.T,np.dot(x_train_set,w)-y_train_set))/(loss*len(x_train_set))
        adagrad += (gradient ** 2)
        w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
#保存参数w
np.save('weight.npy',w)

迭代的次数：0 ， 损失值：27.239592
迭代的次数：100 ， 损失值：598.991742
迭代的次数：200 ， 损失值：96.973083
迭代的次数：300 ， 损失值：240.807182
迭代的次数：400 ， 损失值：71.607934
迭代的次数：500 ， 损失值：212.116933
迭代的次数：600 ， 损失值：117.461546
......
迭代的次数：29400 ， 损失值：19.663413
迭代的次数：29500 ， 损失值：19.694863
迭代的次数：29600 ， 损失值：19.729828
迭代的次数：29700 ， 损失值：19.764402
迭代的次数：29800 ， 损失值：19.794056
迭代的次数：29900 ， 损失值：19.814439

到此为止，模型已经训练完成，接下来对测试集进行测试
四、模型求解
​	首先对test.csv进行数据预处理
#读取数据
testdata = pd.read_csv('test.csv',header = None)
test_data = testdata.iloc[:,2:]
test_data[test_data == 'NR'] = 0
test_data = test_data.to_numpy()
test_x = np.empty([240,18*9],dtype = float)
#数据格式转化，转化为240*162
for i in range(240):
    test_x[i,:] = test_data[18*i:18*(i+1),:].reshape(1,-1)
#数据标准化
for i in range(len(test_x)):
    for j in range(len(test_x[0])):
        if std_x[j] != 0:
            test_x[i][j] = (test_x[i][j] - mean_x[j]) / std_x[j]
#扩展出bais
test_x = np.concatenate((np.ones([240,1]),test_x),axis = 1).astype(float)
print(test_x)

output
[[ 1.         -0.24447681 -0.24545919 ... -0.67065391 -1.04594393
   0.07797893]
 [ 1.         -1.35825331 -1.51789368 ...  0.17279117 -0.10906991
  -0.48454426]
 [ 1.          1.5057434   1.34508393 ... -1.32666675 -1.04594393
  -0.57829812]
 ...
 [ 1.          0.3919669   0.54981237 ...  0.26650729 -0.20275731
   1.20302531]
 [ 1.         -1.8355861  -1.8360023  ... -1.04551839 -1.13963133
  -1.14082131]
 [ 1.         -1.35825331 -1.35883937 ...  2.98427476  3.26367657
   1.76554849]]

最终使用模型进行预测
#在验证集上进行验证
w = np.load('weight.npy')
x_validation= np.concatenate((np.ones([len(x_validation),1]),x_validation),axis = 1).astype(float)
Loss = np.sqrt(np.sum(np.power(np.dot(x_validation,w)-y_validation,2))/len(x_validation))
print ("the Loss on val data is %f" % (Loss))
#预测
ans_y = np.dot(test_x, w)
print('预测PM2.5值')
print(ans_y)

output
the Loss on val data is 18.427185
预测PM2.5值
[[-15.78367116]
 [ -2.32261409]
 [ 59.74234153]
 [ -2.69635112]
 [ 39.23820506]
 [ 13.8801302 ]
......
 [ 54.58686368]
 [ 44.25725877]
 [ 24.03437216]
 [ 50.94458229]
 [ 95.33833139]
 [ 67.19160767]
 [ -9.79522482]
 [ 26.46710166]]

可以看出模型预测的结果并不优秀，甚至出现了负值，毕竟我们只使用了最简单的回归模型，这也就是我们需要为什么学习更复杂的模型的原因


五、反思
​	通过这次作业的不断copy，我学习了机器学习模型从数据处理到预测结果这一系列的整个流程，其中在本次学习中我认为有两个难点，一是我对开始数据预处理的操作的不熟练，所以我称它为骚操作，说它骚，是因为它似乎额外生成了数据，把一个月的所有数据拼接在一起，有效利用了每天之间的数据，也因为它使用numpy进行处理的方式真的很骚，也看出我对numpy数据处理的不熟悉，日后需要进行复习，第二个难点是对含有矩阵或者向量的方程进行求导，并使用代码的方式把他们表示出来不熟悉，这一点希望在以后的练习中能完全掌握。