基于数据驱动的大气PM2.5遥感定量反演(一)

 说明

采用定量遥感领域中一种常用的较为简单的回归算法——线性回归 （Linear Regression）。线性回归通过构建卫星观测和目标参量之间的线性关系， 求解损失函数最小时的参数，实现由卫星观测到目标参量之间的反演。利用线性 回归模型构建 AOD 和 PM2.5 之间的关系，模型结构可表示为： 𝑃𝑀2.5 = 𝛽0 + 𝛽1𝐴𝑂𝐷

---

“PM2.5 浓度.xlxs”为 2015 年站点观测数据，包含站点的经纬度和 PM2.5 年均 值。 “AOD_Wuhan.tif”为 2015 年 MODIS 传感器的 AOD 年均值数据（GeoTIFF 格式；缺省值为-9999，不参与训练和反演），包含地理信息。

---

流程

星地匹配：首先读取数据：使用 geotiffread（）函数读取 AOD 数据，使用 xlsread（） 函数读取 PM2.5 数据。 接着进行坐标转换：使用 map2pix（）函数将地理左标转换到影像空间坐标， 再使用 round（）函数进行取整 然后进行数据对应：将 PM2.5 对应位置上的 AOD 值提取出来 最后剔除缺省值：去掉其中-9999 的值。

关系建模：首先进行线性回归：使用 regress（）函数进行线性回归。 然后进行精度评定：先有拟合的模型计算拟合值，然后计算相关系数和均方 误差 MSE，还画了拟合的结果图来更好的反应线性回归的效果。

时空制图：首先进行反演：利用拟合的模型婴应用到 tif 对应的区域，计算其所对应的 PM2.5 值。 再将计算的结果制图，得到和 AOD 对应的 PM2.5 的图。

 代码

%% 星地匹配
%数据读取
[AOD,R]=geotiffread('AOD_Wuhan.tif'); %AOD 为数据矩阵，R 为地理信息
info=geotiffinfo('AOD_Wuhan.tif');%info 可以辅助坐标转换
[num,txt,raw]=xlsread('PM2.5 浓度.xlsx');%读取 PM2.5,num 为数字，txt 为文字，raw 为
%坐标转换
lon=num(:,3);lat=num(:,2);
[i,j]=map2pix(info.RefMatrix,lon,lat);%i,j 为小数，需取整
i=round(i);j=round(j);
%数据对应
Y=num(:,4);%PM2.5
X=[];
for n=1:size(i,1)
 X(n,1)=AOD(i(n,1),j(n,1));%AOD
end
%去掉缺省值
num1=0;%计数
for n=1:size(X,1)
 if X(n,1)>0
 num1=num1+1;
 X1(num1,1)=X(n,1);
 Y1(num1,1)=Y(n,1);
 end
end
%% 关系建模
%线性回归
X1_=[ones(length(X1),1),X1];%常数项对应 1
[b,~,~,~,~]=regress(Y1,X1_);%b 即为线性回归系数
%精度评定
Y1_fit=X1_*b;%计算拟合值，反演也是类似的操作
r=corr(Y1,Y1_fit);%计算真实值和拟合值的相关系数 r
%MSE
SUM=0;
for n=1:size(X1,1)
 SUM=SUM+(Y1(n,1)-Y1_fit(n,1))*(Y1(n,1)-Y1_fit(n,1));
end
MSE=SUM/n;
%结果图
figure(1)
plot(X1,Y1,'r*');%真值用红色*
hold on
plot(X1,Y1_fit,'b.');%预测值用蓝色.
xlabel('AOD')
ylabel('PM2.5')
%% 时空制图
%反演
[row,col]=size(AOD);
pm=ones(row,col)*nan;
for m=1:row
 for n=1:col
 if AOD(m,n)==-9999 %缺失值跳过，对应 pm2.5 为 nan
 continue
 end
 x_=[1,AOD(m,n)];%常数项对应 1
 pm(m,n)=x_*b;%计算反演值
 end
end
%制图
figure(2)
h=imagesc(pm); %pm 为反演得到的 PM2.5 浓度矩阵,示例里缺失值设为 nan
set(h,'alphadata',~isnan(pm)) ;
colorbar cool;%色带可以自行调整 例如：colorbar cool