利用MATLAB进行逐像元的VFC线性回归(一) 一些准备工作

MATLAB线性回归

1.regress函数

线性回归其实就是通过拟合的方法求出离散点的一元线性方程，大多数采用的是最小二乘方法。最后能求出 $k\text{和}b$。
即 $$y=kx+b$$ $$[b,bint,r,rint,stats]=regress(y,X)$$
  这里直接用最为常见的5个参数的regress 其中b是一个一行两列的向量，第一个返回的是常数项，第二个返回的是 $k$ 。$stats$ 由4个数构成，第一个是 $R^2$，用来表示这个回归模型是否良好，第二个数是 $F$ 分布自由度对应的 $F$ 值，主要用来做 $F$ 检验用，通过这个值大于自由度下的 $F_{0.05}(m,n)$ 以及 $F_{0.01}(m,n)$ 进行比较做F检验，如果利用第三个参数 $p$ 值来做检验，如果 $0.01<p<0.05$ 表示回归的比较显著。如果$p<0.01$则非常显著。
  例如下面的c是10年的VFC植被覆盖度，x为年份，建议用1-n年份，这样回归后的系数会比较大一定
  实际上，就是一个假设检验的问题，我做一个 $H_0$ 的假设和一个 $H_1$ 的备择假设，而 $F$ 分布是在 $y$ 轴右边的，只要找到合适的区间，判断是否在区间内即可，这部分内容需要参考概率论的知识。

c = [0.729587018	0.7546345	0.718326032	0.725992024	0.579145432	0.741261721	0.674983203	0.649813294	0.442957938	0.716242433]';  
x = linspace(1,10,10)';	
X = [ones(size(c)),x]; 							
[b,bint,r,rint,stats] = regress(c,X);              
b1 = [b(2),b(1)];										
k = num2str(b1(1));
pianyi = num2str(b1(2));
y1 = polyval(b1,x1);
plot(x1,c','r*','markersize',8);
R2 = num2str(stats(1));
box off;                                           %只显示两条坐标轴
hold on;
plot(x1,y1,'b-');                                   
text(5,0.6,['y =',k,'x','+',pianyi],'Interpreter','latex','fontsize',12);
str = ['$R^2$ = ',R2];                                                                   
text(2010,0.87,str,'Interpreter','latex','fontsize',12);			%LaTeX的解释器
title('平度站点VFC10年线性回归的结果','fontsize',20,'fontname','楷体');	
xlabel('年份','fontsize',15);
ylabel('VFC','fontsize',15);
set(gca,'XTickLabel',...
    {'2009','2010','2011','2012','2013','2014','2015','2016','2017','2018'});

比如上面这段代码利用线性回归求平度站点10年来的回归系数，需要注意的是，在用regress的时候你需要在前面加入ones(size())用来产生常数，同时b返回的顺序也需要注意，这段代码主要是用来做线性回归并且再作图的。

做完图的效果如上。也可以看到其实一元线性回归的样子其实不是特别好,毕竟是特别简单的函数式来拟合不是特别理想

2.Geotiffread && 等等

[a,R] = geotiffread('回归后的图.tif');
info = geotiffinfo('回归后的图.tif');
geotiffwrite('test.tif',a,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

这三句话是之后要用的核心，对于遥感数据或者地学数据来说也是核心，其中[a,R]中的a是图像的矩阵,R是一些地理框架，而info返回的是地理信息。。

info是一个结构体，R是特定的数据类型

这个是R数据的截图，里面有各种各样的参数，可以忽略之。

这个是info里面的信息，里面有我们经常用到的包括投影信息等等。其中这个GeoTIFFTags非常的重要，对于想要写出包含地理信息的图像的话必须要加上它里面的GeoKeyDirectoryTag。否则会报错，我之前在做一个最大值滤波的尝试的时候就遇到这个问题，后来查到了。
但是写出的数据由于是geotiff格式的，可以用ENVI打开。

在单个像元的基础和有了这些地理函数的基础上，就可以做逐个像元的线性回归。主要是借鉴了某大神的方法，自己改进了一下。