基于PROSAIL模型的植被参数敏感性分析
1 核心目的
(1)掌握PROSAIL模型,分析不同参数对应的光谱特征;
(2)学习模型敏感性分析方法,以夏玉米为例,分析PROSAIL模型各参数的敏感性;
2 前期准备
安装Matlab;安装SimLab;获取PROSAIL模型的参数敏感性分析代码工具。
3 模型介绍
本次实习用到的模型为PROSAIL模型,综合课程讲解和文献资料,现对模型介绍如下:
随着光学遥感技术的发展,辐射传输模型致力于理解植物冠层对光的截获,并根据生物物理特征解释植被反射率。PROSAIL模型便是一种经典的辐射传输模型,由PROSPECT叶片光学特性模型和 SAILH 冠层结构模型耦合得到[1],目前已被广泛用于研究太阳辐射域内植物冠层的光谱和方向反射率,以及生物物理化学参数反演等研究。

图 PROSAIL模型参数示意图
PROSPECT模型是基于最初的Allen平板模型发展起来的叶片尺度光学模型,通过将叶片抽象为N层平板介质,模拟入射光在叶片层面的反射、吸收和透射过程。其中,第一层即叶片表皮,被认为是非各向同性的,而在也在的内部,被认为是各向同性的。该模型主要用来模拟叶片400-2500nm的光学特性。SAILH模型是在SAIL模型的基础上,通过引入热点效应的概念发展起来的冠层尺度的辐射传输模型,将植被当作混合介质,假设叶片方位角分布均匀,考虑任意的叶片倾角,模拟冠层的双向反射率[2]。PROSAIL模型将PROSPECT模型输出的叶片反射率透射率作为SAILH模型的输入进行模型耦合,考虑了植被对太阳辐射的吸收、二向反射、叶片反射透射率、叶片结构参数、地表状况等因素的影响,能较为真实地模拟植被冠层情况,并且具有很高的稳定性[3]。
因此,PROSAIL的核心是将PROSPECT模拟的叶片反射/透射特性作为输入传给SAIL,从而由叶片“小单元”的光学表现,拓展到冠层“大群体”的反射率模拟,实现微观到宏观的结合。
4 参数介绍
PROSAIL模型需要的参数主要由PROSPECT 模型和SAILH模型的参数组成。其中,PROSPECT模型的输入参数包括叶片结构参数(N),叶绿素含量(Cab)、类胡萝卜素含量(Car)、褐色素含量(Cbrown)、叶片含水量(Cw)和叶片干物质含量(Cm)共六个参数。而SAILH模型的输入参数包括叶面积指数(LAI)、平均叶倾角(ALA)、热点参数(hspot)、土壤反射率(ρsoil)、天空光散射比(SKYL)、太阳天顶角(tts)、观测天顶角(tto)和相对方位角(psi)共8个参数。
表1 PROSAIL模型输入参数及取值范围
|
参数 |
参数描述 |
单位 |
默认值 |
最小值 |
最大值 |
|
|
PROSPCT模型参数 |
Cab |
叶绿素含量 |
μg·cm-2 |
58 |
5 |
80 |
|
Car |
类胡萝卜素含量 |
μg·cm-2 |
14.5 |
1 |
20 |
|
|
Cm |
干物质含量 |
g·cm-2 |
0.003 |
0.001 |
0.005 |
|
|
Cw |
叶片水分含量 |
cm |
0.013 |
0.003 |
0.03 |
|
|
Cbrown |
褐色素含量 |
/ |
0 |
/ |
/ |
|
|
N |
叶片结构参数 |
/ |
1.5 |
1 |
24 |
|
|
SAILH 模型参数 |
LAI |
叶面积指数 |
m2·m-2 |
3 |
1 |
8 |
|
ALA |
平均叶倾角 |
° |
50 |
30 |
60 |
|
|
hspot |
热点参数 |
/ |
0.25 |
0.25 |
1 |
|
|
ρsoil |
土壤反射率 |
/ |
/ |
/ |
/ |
|
|
SKYL |
太空光散射比 |
/ |
/ |
/ |
/ |
|
|
tts |
太阳天顶角 |
° |
30 |
0 |
90 |
|
|
tto |
观测天顶角 |
° |
0 |
0 |
90 |
|
|
psi |
相对方位角 |
° |
0 |
0 |
180 |
|
本次上机实习,以夏玉米冠层为研究对象,设定 PROSAIL 模型输入参数见表 1,选取了叶片结构参数,叶绿素含量、类胡萝卜素含量、叶片含水量和叶片干物质含量和叶面积指数六个参数进行敏感性分析。通过查阅文献资料,根据玉米的特征,将热点参数默认值设为 0.25,其取值范围为 0.25~1。除此之外,由于褐色素含量与叶绿素含量相冲突,且其含量不易测量,因此将褐色素含量设为 0。
5 参数敏感性分析
5.1 局部敏感性分析
(1)方法介绍
局部敏感性分析采用的方法是 OTA(one factor at a time),基于参数相互独立的假设,通过预先设定一个参数的多个值,固定其他值的不变,来分析参数变化对模型输出结果的影响[2]。本次上机实习对 叶片结构参数,叶绿素含量、类胡萝卜素含量、叶片含水量和叶片干物质含量和叶面积指数六个参数冠层参数进行敏感性分析,每次只改变一个参数的值,其他参数取表中的默认值,观察输入参数的改变对冠层反射光谱在 400-2500 nm 波段的影响。
(2)实习步骤
调用外部模型“main_PROSAIL_5B”,打开Matlab代码。
clc;
% 设置不同的Cw值
LAI_values = [1, 3, 8]; % 你可以根据需要添加更多Cw值
% 颜色数组,用于区分不同的Cw值
colors = ['r', 'g', 'b']; % 红色,绿色,蓝色
hold on;
for i = 1:length(LAI_values)
LAI = LAI_values(i);
% 设置其他参数
TypeLidf = 2; % Ellipsoidal LIDF
if (TypeLidf == 1)
LIDFa = -0.35;
LIDFb = -0.15;
elseif (TypeLidf == 2)
LIDFa = 57;
LIDFb = 0;
end
% 叶片化学成分和结构参数
Cab = 58;
Cm = 0.003;
Cbrown = 0.0;
Car = 14.5;
Cw = 0.013;
N = 1.5;
data = dataSpec_P5B; % 读取光谱数据
% 土壤反射率
Rsoil1 = data(:, 10);
Rsoil2 = data(:, 11);
psoil = 1;
rsoil0 = psoil * Rsoil1 + (1 - psoil) * Rsoil2;
% 树冠结构参数
hspot = 0.25;
tts = 30.0;
tto = 10.0;
psi = 0.0;
% 调用PRO4SAIL模型
[rdot, rsot, rddt, rsdt] = PRO4SAIL(N, Cab, Car, Cbrown, Cw, Cm, LIDFa, LIDFb, TypeLidf, LAI, hspot, tts, tto, psi, rsoil0);
% 计算不同光源的辐射
Es = data(:, 8);
Ed = data(:, 9);
rd = pi / 180;
skyl = 0.847 - 1.61 * sin((90 - tts) * rd) + 1.04 * sin((90 - tts) * rd) * sin((90 - tts) * rd);
PARdiro = (1 - skyl) * Es;
PARdifo = skyl * Ed;
% 计算方向反射率
resv = (rdot .* PARdifo + rsot .* PARdiro) ./ (PARdiro + PARdifo);
% 绘制不同LAI值下的光谱曲线
plot(data(:, 1), resv, 'linewidth', 2, 'Color', colors(i));
end
% 设置图表属性
axis([400 2500 0 1]);
xlabel('Wavelength (nm)');
ylabel('Reflectance');
title('不同LAI值光谱曲线');
grid on;
legend(arrayfun(@(x) sprintf('LAI = %.3f', x), LAI_values, 'UniformOutput', false));
hold off;
源代码只能够一次性输入一个值,因为我们想探索同一个参数值的变化对曲线响应的影响,这里对代码进行调整,改成循环结构。以Cw为例,固定其他参数不变。

图 Matlab图片编辑界面

图 光谱曲线导出结果
其他的参数采用类似的方式进行曲线绘制,得到所有指数的敏感性分析结果。
(3)局部敏感性分析结果
对于叶片中的色素来说,叶绿素含量(Cab)和类胡萝卜素(Car)对冠层反射光谱的影响主要集中在可见光波段。

图 Cab局部敏感性分析
其中,对Cab参数最敏感的为510~690nm,对Car最敏感的为480~560nm,但从光谱曲线来看,整体Cab参数对反射率的影响程度比Car大。且随着Car和Cab等这类叶片色素的含量增加,玉米冠层反射率减小。

图 Cm局部敏感性分析

干物质含量(Cm)影响主要集中在760~1300nm、1560~1860nm、2100~2300nm三个反射峰波段,但是整体影响较小,与反射率呈反相关。叶片水分含量(Cw)参数的影响则涵盖800~2400nm的宽波段范围,在该范围内随着叶片水分含量的增加,叶片冠层的反射率下降。因此叶片的含水量在一个较宽的波段内对植被叶片发射率都有较大影响。

图 N局部敏感性分析
从图中可以看出,叶片结构参数N没有表现出特定波段的敏感性,在400-2500nm的波段范围内整体呈现正相关。

5.2 全局敏感性分析
(1)方法介绍
全局敏感性分析主要用于分析每个参数以及参数之间相互作用对模型结果的影响。本次实习,结合SimLab软件和Matlab软件,使用EFAST方法对相关参数进行敏感性分析。EFAST方法是Sailtelli等结合Sobel法和傅立叶幅度敏感性检验法的优点提出的一种全局敏感性分析方法。该方法稳健、计算高效并且需要的样本数较低。
本次实习,主要使用SimLab软件,并且调用Matlab的全局敏感性分析方法,对玉米冠层的几个参数进行全局敏感性分析。
(2)实习步骤
打开SimLab软件,按照Configure→Create New→点击新建参数→设置Name→设置分布类型为Uniform即均匀分布→设置参数的范围→Apply→OK,完成参数设置。

图 参数设置
其他参数的设置和上面类似。
所有参数设置完成后,点击Accept Factors,完成参数创建

图 参数创建
接下来,需要生成参数样本,设置方法为Fast,点击旁边的Specify Switches→选择第二种方法,设置种子seed为10,Number of Executions为2048,点击OK,完成。

图 EFAST方法生成参数样本
设置输出的样本文件,点击Sample File Name,将文件命名为corn.sam,一定要注意,命名的后缀为.sam,否则无法进行后续的操作

图 样本文件命名设置
完成上面的配置后,点击Generate,便能生成样本文件

图 文件夹中的样本文件

图 样本生成
修改头文件,打开刚生成的corn.sam文件,讲参数样本文件的头部和尾部信息删除,并且保存为新的corn1.sam,这里一共生成了2022个点,共计6个参数,每一列一个参数。

图 修改文件的尾部
调用外部模型,打开Matlab,选择“打开”,调用“global_sensibility.m”分析。
clear
tic
clc
test=importdata("E:\PROSAIL Model\SimLab\corn1.sam");
Cab=test(:,1);
Car=test(:,2);
Cm=test(:,3);
Cw=test(:,4);
N=test(:,5);
LAI=test(:,6);
TypeLidf=2;
% if 2-parameters LIDF: TypeLidf=1
if (TypeLidf==1)
% LIDFa LIDF parameter a, which controls the average leaf slope
% LIDFb LIDF parameter b, which controls the distribution's bimodality
% LIDF type a b
% Planophile 1 0
% Erectophile -1 0
% Plagiophile 0 -1
% Extremophile 0 1
% Spherical -0.35 -0.15
% Uniform 0 0
% requirement: |LIDFa| + |LIDFb| < 1
LIDFa = -0.35;
LIDFb = -0.15;
% if ellipsoidal LIDF: TypeLidf=2
elseif (TypeLidf==2)
% LIDFa = average leaf angle (degrees) 0 = planophile / 90 = erectophile
% LIDFb = 0
LIDFa = 30;
LIDFb = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LEAF CHEM & STR PROPERTIES %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cab = 58; % chlorophyll content (礸.cm-2)
%Car = 10; % carotenoid content (礸.cm-2)
Cbrown = 0.0; % brown pigment content (arbitrary units)
%Cw = 0.0131; % EWT (cm)
%Cm = 0.003662; % LMA (g.cm-2)
%N = 1.518; % structure coefficient
data=dataSpec_P5B;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Soil Reflectance Properties %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rsoil1 = dry soil
% rsoil2 = wet soil
Rsoil1=data(:,10);Rsoil2=data(:,11);
psoil = 1; % soil factor (psoil=0: wet soil / psoil=1: dry soil)
rsoil0=psoil*Rsoil1+(1-psoil)*Rsoil2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4SAIL canopy structure parm %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LAI = 3.; % leaf area index (m^2/m^2)
hspot = 0.25; % hot spot
tts = 30.; % solar zenith angle (?
tto = 0.; % observer zenith angle (?
psi = 0.; % azimuth (?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALL PRO4SAIL %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rdot: hemispherical-directional reflectance factor in viewing direction
% rsot: bi-directional reflectance factor
% rsdt: directional-hemispherical reflectance factor for solar incident flux
% rddt: bi-hemispherical reflectance factor
ref=[];
for i=1:2022
[rdot,rsot,rddt,rsdt]=PRO4SAIL(N(i),Cab(i),Car(i),Cbrown,Cw(i),Cm(i),LIDFa,LIDFb,TypeLidf,LAI(i),hspot,tts,tto,psi,rsoil0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% direct / diffuse light %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the direct and diffuse light are taken into account as proposed by:
% Francois et al. (2002) Conversion of 400?100 nm vegetation albedo
% measurements into total shortwave broadband albedo using a canopy
% radiative transfer model, Agronomie
% Es = direct
% Ed = diffuse
Es=data(:,8);Ed=data(:,9);
rd=pi/180;
skyl = 0.847- 1.61*sin((90-tts)*rd)+ 1.04*sin((90-tts)*rd)*sin((90-tts)*rd); % % diffuse radiation
PARdiro = (1-skyl)*Es;
PARdifo = (skyl)*Ed;
% resv : directional reflectance
resv = (rdot.*PARdifo+rsot.*PARdiro)./(PARdiro+PARdifo);
ref=[ref resv];
end
dlmwrite('Refl_test.txt',[data(:,1),ref],'delimiter','\t','precision',5)
red=ref(645-400+1,:)';
green=ref(555-400+1,:)';
blue=ref(450-400+1,:)';
nir=ref(859-400+1,:)';
swir=ref(1640-400+1,:)';
save('Red_1.dat','red','-ascii');
save('Green_1.dat','green','-ascii');
save('Blue_1.dat','blue','-ascii');
save('NIR_1.dat','nir','-ascii');
save('SWIR_1.dat','swir','-ascii');
toc
添加测试参数,加上Car=test( :2)。

图 代码修改1
调用PROSPECT模型,注释掉的参数,包括叶绿素、水分、干物质和结构以及类胡萝卜素,都是从采样点文件中读取,其他的参数直接设为固定值。

图 代码修改2
调用SAILH模型,讲LAI参数注释掉,将hspot的默认值修改为0.25。

图 代码修改3
将循环的行数变成2022,这与前面的文本文件是一致的。

图 代码修改4
运行代码,最后生成红、绿、蓝、近红外、短波红外这几个反射率的点文件。
由于生成的波段反射率文件用来进行蒙特卡洛模拟,需要对每一份波段反射率文件的开头增加一些信息,以便SimLab软件能够正常识别。

图 编辑波段头文件信息
继续回到SimLab软件,选择Load Sample File,这里加载sam的原文件,注意不是修改了头部和尾部的文件。点击configure(monte carlo),此时出现黄色小箭头,点击“Select Model”,在窗口中选择蓝波段反射率文件,选中后点击OK。

图 蒙特卡洛模拟和敏感性分析1
点击Start(monee carlo),点击Analyse(UA/SA),进行敏感性分析

图 蒙特卡洛模拟和敏感性分析2
按照选择ref变量→点击Add→在弹出的窗口中点击New Variable→点击SA进行敏感性分析步骤进行操作。

图 蒙特卡洛模拟和敏感性分析3

图 蒙特卡洛模拟和敏感性分析4
运行后,我们能够得到两个结果,在Chart中用饼状图展示了各个参数的贡献率,在Iabulated Value中展示了具体的贡献率值。其中,First Order是一阶敏感性指数,反映每个参数对模型结果的直接贡献率;Total Order是总敏感性指数,反映参数直接贡献率以及参数间的交互耦合作用对模型输出结果的贡献率之和。
将分析结果输入到Execl表格中进行记录,以便后续进行绘图分析。

图 全局敏感性分析结果1

图 全局敏感性分析结果2
其他波段的反射率文件也是参照这个过程进行分析,得到所有波段的数据后,整理并绘制柱形图,对比不同波段下,各个参数的全局敏感性如何。
接下来,打开Origin软件,将每个波段的数据输入进表格,使用分组柱形图绘制,调整字体和图形大小等,最终得到每个波段的全局敏感性分析结果。

图 利用Origin进行全局敏感性分析制图
(3)全局敏感性分析结果
根据上面的步骤,对玉米冠层生理结构参数进行全局敏感性分析,对400~2000nm的波段范围结果进行输入输出均匀采样,计算各输入参数的一阶敏感性和全局敏感性,得到了蓝光波段,绿光波段,红光波段,近红外波段以及短波红外波段的结果。

图 蓝光波段参数全局敏感性分析结果
从蓝光波段来看,叶片结构参数N的敏感性最高,一阶和全局敏感性参数能够超过0.5,其次为Cab,其全局敏感性在0.3左右。
从绿光波段来看,叶绿素含量Cab的敏感性最高,一阶和全局敏感性参数能够超过0.5,其次为N,其全局敏感性也能够达到0.5。

图 绿光波段参数全局敏感性分析结果

图 红光波段参数全局敏感性分析结果
红光波段的敏感性排序与绿光波段类似,Cab的敏感性最高,基本能够解释超过一般的反射率的方差变化,这主要由于叶绿素吸收红光反射绿光的特点造成的[4]。
近红外波段中,冠层叶面积指数的敏感性断层式领先其他参数,一阶和全局敏感性指数可超达0.8,主要由于近红外光波段的光主要受到植物叶片的结构影响。

图 近红外波段参数全局敏感性分析结果

图 短波红外波段参数全局敏感性分析结果
短波红外中,叶片含水量Cw的敏感性最高,一阶和全局敏感性指数可超达0.5,这主要由于短波红外波段对水分的吸收非常敏感,水分含量的变化对短波红外波段的反射率影响非常大。其次为叶片结构指数N,敏感性也能够达到0.45。
以上便是结合SimLab软件开展参数敏感性分析的全过程,创作不易,如果对你的学习有帮助,欢迎关注我们!您的肯定是我们更新最大的动力。
参考文献:
- JACQUEMOUD S, VERHOEF W, BARET F, et al. PROSPECT+SAIL models: A review of use for vegetation characterization[J]. Remote Sensing of Environment, 2009, 113: S56-S66.
- 王李娟,牛铮.PROSAIL模型的参数敏感性研究[J].遥感技术与应用,2014,29(02):219-223.
- 苏伟,邬佳昱,王新盛,等.基于Sentinel-2影像与PROSAIL模型参数标定的玉米冠层LAI反演[J].光谱学与光谱分析,2021,41(06):1891-1897.
- 李青.基于PROSAIL模型的玉米冠层参数敏感性分析[J].科技创新与应用,2023,13(23):18-22.
- 肖艳芳,周德民,宫辉力,等.冠层反射光谱对植被理化参数的全局敏感性分析[J].遥感学报,2015,19(03):368-374.
- 何维,杨华.模型参数全局敏感性分析的EFAST方法[J].遥感技术与应用,2013,28(05):836-843.
更多推荐



所有评论(0)