1.4 时间序列遥感数据趋势分析

Posted 2021-05-02 GIS新视野

1.4 时间序列遥感数据趋势分析

1、分析原理

将三维矩阵逐元素逐页按照时间顺序提取，即相当于提取时间序列遥感影像（多波段）逐像元逐时刻的像元值，提取出的多个值即为该像元时间序列数据，经过转换和计算，得到时间序列分析结果，并将其写入栅格中。示意图如下：

2、时间序列遥感数据趋势分析

趋势分析中采用regress函数对时间序列建立一元线性回归方程，从而确定出回归系数（斜率与截距）以及p检验值。regress调用格式如下：

[B, Bint, E, Eint, Stats]=regress(Y, X, α)

其中，X表示自变量观测值，Y为因变量观测值，α为显著水平，默认为0.05，B为回归系数的最小二乘估计向量，Bint为回归系数的（1-α）%的区间估计，E为残差向量，Eint为残差向量的区间估计，Stats为返回R2统计量、F统计量、F统计量对应的sig概率值和均方差MSE。在B向量中，B(1)为截距，B(2)为斜率

示例：

x=[15.2 10.4 21.2 18.6 26.4 23.4 13.5 16.7 24 19.1];

y=[28.6 19.3 40.5 35.6 48.9 45 29.2 34.1 46.7 37.4];

[m,n]=size(x); %计算自变量矩阵行列数

X=[ones(n,1),x']; %在自变量矩阵中添加常数向量并转置

Y=y'; %因变量向量转置

[B,Bint,E,Eint,Stats]=regress(Y,X);

从三维矩阵中提取出的逐像元时间序列维数为1行1列T页，仍为三维形式，需要将其重塑转换为一列或者一行，在对其进行计算。矩阵重塑见“语法基础”中的“矩阵操作”。

%假设时间序列栅格生成的三维矩阵、投影信息和参考矩阵变量分别为data、proj与R

[m,n,p]=size(data); %m:行数，n:列数，p:页数

x=1:p; %序列的时间标签

X=[ones(p,1),x']; %构造X

for i=1:m %页循环

for j=1:n %行循环

Y=reshape(y(i,j,:),[p,1]); %三维转一列构造Y

[B,Bint,E,Eint,Stats]=regress(Y,X);

slope(i,j)=B(2); %斜率

pval(i,j)=Stats(3); %p检验值

intercept(i,j)=B(1); %截距值

end

end

%将slope、pval、intercept写入TIFF文件方法，参考1.2