一种基于辐射传输模型反演的叶片氮含量估算方法

1.一种基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，包括：
采集两组位于不同大洲的森林及草地生态系统的叶片样本，分别获取两组所述叶片样本对应的光谱反射率实测值和叶片氮含量实测值，分别定义为第一数据组和第二数据组；
将镜面反射光谱辐射纳入叶片辐射传输模型PROSPECT建立反演模型；
所述反演模型中的短波红外光谱波段按照预设长度间隔分为若干个光谱子集，以不同波段的所述光谱子集作为初始条件，分别将所述第一数据组和所述第二数据组对应的所述光谱反射率实测值代入独立的两个所述反演模型进行初步反演，得到各自的第一叶片氮含量估测值，计算两个所述第一叶片氮含量估测值与对应的所述叶片氮含量实测值之间的第一均方根误差，以精度最高的所述第一均方根误差为筛选条件，通过后向迭代递归从所述光谱子集中选取各自的半离散最优光谱集合；
将两组所述半离散最优光谱集合交叉代入已进行过所述初步反演的反演模型进行交叉验证，得到各自的第二叶片氮含量估测值，计算两个所述第二叶片氮含量估测值与对应的所述叶片氮含量实测值之间的第二均方根误差，若所述第二均方根误差处于预设值内，输出所述反演模型。
2.如权利要求1所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，所述反演模型的建立过程包括：
将所述镜面反射光谱辐射纳入所述叶片辐射传输模型PROSPECT进行正演，获得的光谱反射率模拟值用于建立代价函数；
将所述光谱反射率实测值代入所述代价函数进行初步反演，得到叶片结构参数估测值和镜面反射光谱辐射估测值；
将所述叶片结构参数估测值和所述镜面反射光谱辐射估测值作为先验知识重新代入所述代价函数，通过归一化处理所述先验知识，确定所述先验知识在反演中所占权重，所述权重用于补全所述代价函数。
3.如权利要求2所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，所述光谱反射率模拟值的计算方法为：
其中，λ为光谱波长， ， ，Nstruc, Cab, Cxc, Cw, Narea, Carea, Rs分别代表入射光线角度，辐射高度角度，叶片结构参数，叶绿素含量，类胡萝卜含量，水分含量，叶片氮含量，含碳化合物含量和镜面反射光谱辐射。
4.如权利要求3所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，所述代价函数为：
其中，θ为模型参数向量， 和 代表代价函数的光谱波长起始值和终点值；BRFmes(λ)为光谱反射率实测值；BRFsim(λ)为光谱反射率模拟值；通过所述代价函数运算得到所述叶片结构参数估测值 和所述镜面反射光谱辐射估测值 。
5.如权利要求4所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，经过补全的所述代价函数为：
其中， 为叶片结构参数最大值， 为叶片结构参数最小值， 为
镜面反射光谱辐射最大值， 为镜面反射光谱辐射最小值， 为叶片结构参数的权重，为镜面反射光谱辐射的权重。
6.如权利要求1所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，所述短波红外光谱波段设置为1400‑2399nm，按照50nm的长度间隔分为20个所述光谱子集。
7.如权利要求6所述的基于辐射传输模型反演的叶片氮含量估算方法，其特征在于，所述后向迭代递归从所述光谱子集中选取各自的半离散最优光谱集合的过程包括：以50nm的长度间隔分为20个所述光谱子集作为初始条件，所述第一数据组和所述第二数据组分别进行所述初步反演得到各自20个所述光谱子集的第一叶片氮含量估测值，将20个所述光谱子集的第一叶片氮含量估测值分别与对应数据组中的叶片氮含量实测值进行均方根误差计算，并搜索得到第一叶片氮含量估测值精确度最高的光谱子集，将搜索到的光谱子集加入当前的最优光谱集合，并通过结合当前的最优光谱集合与剩余的19个光谱子集，形成新的
19个光谱子集，每个光谱子集包括两个50nm间隔的光谱子集，利用19个不同光谱子集重新进行所述初步反演，并搜索得到当前阶段下第一叶片氮含量估测值精确度最高的光谱子集，依次迭代类推，直至所有的50nm长度间隔光谱子集依次加入到最优光谱集合，计算每个迭代过程中最优光谱集合得到的第一叶片氮含量估测值的均方根误差，并根据均方根误差排序得到均方根误差最小值所对应的光谱集合。

一种基于辐射传输模型反演的叶片氮含量估算方法\n技术领域\n[0001] 本发明涉及植被遥感技术领域，尤其涉及一种基于辐射传输模型反演的叶片氮含量估算方法。\n背景技术\n[0002] 氮元素是植物生长的重要限制因子，可以通过潜在的气候反馈作用影响陆地生态系统的碳变化，同时也是生态系统过程模型的重要输入参数。叶片氮能够调节植物的多项生理过程，主要包括光合作用，叶片呼吸和蒸腾作用，并与冠层和植株水平性状紧密联系比如光能利用效率，木质生长和净初级生产力等。因其在生物多样性和生态系统功能的重要作用，叶片氮含量也被列为关键生物多样性指标之一，是遥感和生态领域均重点关注的研究目标。\n[0003] 叶片氮含量只占叶片干重的0.2‑6.4%，比例很低，基于遥感估测叶片氮含量主要利用了蛋白质和叶绿素的吸收特征，例如可见光、近红外和短波红外区域的光谱信号。由于叶片氮的光谱吸收波段大多位于短波红外区域，常被干物质和水分更强更明显的吸收特征所掩盖，因此叶片氮的遥感反演存在着很高的不确定性。\n[0004] 过去利用高光谱遥感反演叶片氮含量的研究主要采用统计回归方法，这种方法通过建立叶片氮含量和光谱反射率数据的回归方程用来预测，方法简单易用，但是常常与叶片氮实际的吸收特征不一致，导致模型的机理解释度较低；另外，由于回归方程是由数据拟合得到，难以用到其它的独立数据集，导致模型的通用性较低。再者，高光谱遥感反演叶片氮含量方法的辐射传输模型通过基于物理光学定律模拟太阳辐射与叶片的相互作用，包括反射、散射和吸收等过程，来模拟叶片的光谱反射率，而利用模型的反向过程反演估测叶片组分参数。这种方法被广泛的应用于例如叶绿素、水分等具有明显光谱吸收特征的组分，但是用于叶片氮的光谱吸收特征较弱，以及不同叶片参数组合会产生相似光谱反射率，导致出现“病态反演”问题，且反演精度较低。\n[0005] 综上，目前仍然缺乏能够高精度估算叶片氮含量的辐射传输模型反演方法。\n发明内容\n[0006] 针对上述问题，本发明提出一种基于辐射传输模型反演的叶片氮含量估算方法，主要解决现有高光谱卫星遥感反演叶片氮含量方法的反演精度低、模型通用性较低的问题。\n[0007] 为解决上述技术问题，本发明的技术方案如下：\n[0008] 一种基于辐射传输模型反演的叶片氮含量估算方法，包括：\n[0009] 采集两组位于不同大洲的森林及草地生态系统的叶片样本，分别获取两组所述叶片样本对应的光谱反射率实测值和叶片氮含量实测值，分别定义为第一数据组和第二数据组；\n[0010] 将镜面反射光谱辐射纳入叶片辐射传输模型PROSPECT建立反演模型；\n[0011] 所述反演模型中的短波红外光谱波段按照预设长度间隔分为若干个光谱子集，以不同波段的所述光谱子集作为初始条件，分别将所述第一数据组和所述第二数据组对应的所述光谱反射率实测值代入独立的两个所述反演模型进行初步反演，得到各自的第一叶片氮含量估测值，计算两个所述第一叶片氮含量估测值与对应的所述叶片氮含量实测值之间的第一均方根误差，以精度最高的所述第一均方根误差为筛选条件，通过后向迭代递归从所述光谱子集中选取各自的半离散最优光谱集合；\n[0012] 将两组所述半离散最优光谱集合交叉代入已进行过所述初步反演的反演模型进行交叉验证，得到各自的第二叶片氮含量估测值，计算两个所述第二叶片氮含量估测值与对应的所述叶片氮含量实测值之间的第二均方根误差，若所述第二均方根误差处于预设值内，输出所述反演模型。\n[0013] 本发明的有益效果为：通过采集广泛地理区域、物种和生长时期的植物获取叶片光谱反射率数据，采用先验知识和选取模型反演半离散最优光谱集合，并进行反演策略的交叉验证，从而实现高精度高通用性反演估测叶片氮含量。\n附图说明\n[0014] 图1为本发明实施例公开的基于辐射传输模型反演的叶片氮含量估算方法的流程图；\n[0015] 图2a为第一数据组初步反演后第一叶片氮含量估测值与对应的叶片氮含量实测值之间的第一均方根误差的关系示意图；\n[0016] 图2b为第二数据组初步反演后第一叶片氮含量估测值与对应的叶片氮含量实测值之间的第一均方根误差的关系示意图；\n[0017] 图3a为第一数据组采用辐射传输模型反演和半离散最优光谱集合的关系示意图；\n[0018] 图3b为第二数据组采用辐射传输模型反演和半离散最优光谱集合的关系示意图；\n[0019] 图4a为第二数据组的半离散最优光谱集合用于第一数据组的结果示意图；\n[0020] 图4b为第一数据组的半离散最优光谱集合用于第二数据组的结果示意图。\n具体实施方式\n[0021] 为使本发明的目的、技术方案及优点更加清楚、明确，下面结合附图和具体实施方式对本发明的内容做进一步详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部内容。\n[0022] 本实施例提出了一种基于辐射传输模型反演的叶片氮含量估算方法，通过采集广泛地理区域、物种和生长时期的植物获取叶片光谱反射率数据，采用先验知识和选取模型反演半离散最优光谱集合，并进行反演策略的交叉验证，从而实现高精度高通用性反演估测叶片氮含量。该方法如图1所示，包括：\n[0023] S1，采集两组位于不同大洲的森林及草地生态系统的叶片样本，分别获取两组叶片样本对应的光谱反射率实测值和叶片氮含量实测值，分别定义为第一数据组和第二数据组。\n[0024] 在本实施例中，叶片样本基于我国华南地区6个野外站点和北美东部18个野外站点，主要包括亚热带、热带、温带森林和草地生态系统，针对乔灌草的优势物种，获取顶部阳叶枝条，采集约10克鲜叶。植物叶片样本包括不同地理区域、物种、光照条件和生长阶段，旨在覆盖光谱反射率实测值和叶片氮含量实测值的变化区间。最终，得到两个数据集包括我国华南地区93物种363个叶片样本（定义为第一数据组）和北美东部161物种384叶个片样本（定义为第二数据组）。\n[0025] 具体的，利用地物光谱仪和叶片夹测量植物叶片样本的光谱反射率实测值，该测量方式具有固定的入射光线角度和辐射高度角度，能极大的减小测量中存在不确定性，得到高精度和高一致性的叶片光谱反射率。其中，叶面积较大的阔叶树种可直接测量，窄叶草本则通过排列多个叶片完成测量。将已完成光谱测量的叶片样本烘干研磨后，利用元素分析仪得到叶片氮含量实测值。反射率实测值的获取过程包括：\n[0026] S2，将镜面反射光谱辐射纳入叶片辐射传输模型PROSPECT建立反演模型。\n[0027] 反演模型的建立过程包括：\n[0028] S201，将镜面反射光谱辐射纳入叶片辐射传输模型PROSPECT进行正演，获得的光谱反射率模拟值用于建立代价函数；\n[0029] 光谱反射率模拟值的计算方法为：\n[0030]\n[0031] 其中，λ为光谱波长，  ，  ，Nstruc, Cab, Cxc, Cw, Narea, Carea, Rs分别代表入射\n2 2\n光线角度，辐射高度角度，叶片结构参数，叶绿素含量（µg/cm），类胡萝卜含量（µg/cm），水\n2 2 2\n分含量（g/m），叶片氮含量（g/m），含碳化合物含量（g/m）和镜面反射光谱辐射。\n[0032] 在本实施例中，公式（1）中的叶片结构参数，叶绿素含量，类胡萝卜含量，水分含量，叶片氮含量，含碳化合物含量和镜面反射光谱辐射根据实测数据集取值范围，将叶片参数Nstruc, Cab, Cxc, Cw, Narea, Carea, Rs的初始值设为[1.5,30,5,100,2.26,10,0.2]，最小值设为[1,0.5,0.5,10, 0.02,0.1,‑0.2],最大值设为[3,100,20,440,6.77,200,0.6]。\n[0033] S202，将光谱反射率实测值代入代价函数进行初步反演，得到叶片结构参数估测值和镜面反射光谱辐射估测值；\n[0034] 代价函数为：\n[0035]\n[0036] 其中，θ为模型参数向量， 和 代表代价函数的光谱波长起始值和终点值；BRFmes(λ)为光谱反射率模拟值实测值（即，在波长λ处的实测双向反射率因子）；BRFsim(λ)为光谱反射率模拟值（即，在波长λ处利用辐射传输模型参数得到的模拟双向反射率因子）；通过代价函数运算得到叶片结构参数估测值 和镜面反射光谱辐射估测值 。\n[0037] 在本实施例中，公式（2）的光谱波长λ根据全波段光谱反射率设置为400‑2500nm。\n[0038] S203，将叶片结构参数估测值和镜面反射光谱辐射估测值作为先验知识重新代入代价函数，通过归一化处理先验知识，确定先验知识在反演中所占权重，权重用于补全代价函数。\n[0039] 经过补全的代价函数为：\n[0040]\n[0041] 其中， 为叶片结构参数最大值， 为叶片结构参数最小值， 为\n镜面反射光谱辐射最大值， 为镜面反射光谱辐射最小值，为叶片结构参数的权重，为镜面反射光谱辐射的权重。\n[0042] S3，反演模型中的短波红外光谱波段按照预设长度间隔分为若干个光谱子集，以不同波段的光谱子集作为初始条件，分别将第一数据组和第二数据组对应的光谱反射率实测值代入独立的两个反演模型进行初步反演，得到各自的第一叶片氮含量估测值，计算两个第一叶片氮含量估测值与对应的叶片氮含量实测值之间的第一均方根误差，以精度最高的第一均方根误差为筛选条件，通过后向迭代递归从光谱子集中选取各自的半离散最优光谱集合。\n[0043] 在本实施例中，短波红外光谱波段设置为1400‑2399nm，按照50nm的长度间隔分为\n20个光谱子集S=[S1,S2,…,S20]，其中，S1为1400‑1449nm，S2为1450‑1499nm，…，S20为2350‑\n2399nm。模型迭代反演以20个不同的光谱子集S=[S1,S2,…,S20]作为初始条件。\n[0044] 以第一数据组为例，参照步骤S202和S203记载的公式（2）、（3）利用单个光谱子集分别进行初步反演得到针对20个光谱子集的第一叶片氮含量估测值，将20个光谱子集得到的第一叶片氮含量估测值分别与第一数据组中的叶片氮含量实测值进行均方根误差（RMSE）计算，并搜索得到第一叶片氮含量估测值精确度最高的光谱子集Si；其次，将搜索到的光谱子集Si加入当前的最优光谱集合Smax，并通过结合当前的最优光谱集合Smax与剩余的\n19个光谱子集S i，形成新的19个光谱子集Si,j，其中每个光谱子集Si,j包括两个50nm间隔~ ~ ~\n的光谱子集Si；再者，利用这19个不同光谱子集Si,j重新进行步骤S203的反演，并搜索得到~\n当前阶段下第一叶片氮含量估测值精确度最高的光谱子集；依次迭代类推，直至所有的\n50nm长度间隔光谱子集依次加入到最优光谱集合Smax；最后，计算每个迭代过程中最优光谱集合得到的第一叶片氮含量估测值的均方根误差（RMSE），并根据均方根误差排序得到均方根误差最小值所对应的光谱集合Sopt，也即为叶片氮含量最终的半离散最优光谱集合，每个光谱子集进入最优光谱集合的顺序代表其对模型反演的重要性。\n[0045]\n[0046] 表1表明，第一数据组的最优集合为2150‑2199, 2100‑2149, 1850‑1899, 1950‑\n1999, 2200‑2249, 2000‑2049nm，而第二数据组得到的半离散最优光谱集合为2150‑2199, \n1500‑1549, 2100‑2149, 1800‑1849nm，这两者共同含有2150‑2199nm和2100‑2149nm两个光谱子集，这些光谱子集主要位于反射率光谱的短波红外范围，并与叶片氮N‑H等化合键的光谱吸收特征诸如1510、1940、1980、2130、2180和2240nm等保持一致。\n[0047] 图2a,2b分别示意了第一数据组和第二数据组的第一均方根误差与对应的叶片氮含量实测值之间的关系。X轴表示后向迭代递归选取模型反演得到半离散最优光谱子集的顺序，Y轴表示迭代次序。“+”号代表当前迭代次序下得到RMSE最小值的光谱子集；星号表示所有迭代过程中得到RMSE最小值的光谱子集。第一均方根误差整体随迭代次序逐步先降低直到得到最小值，之后添加光谱子集第一均方根误差将升高。第一数据组需要六个光谱子集构成其半离散最优光谱集合，而第二数据组则需要四个光谱子集。\n[0048] 图3a,3b表明，采用辐射传输模型反演和半离散最优光谱集合，第一数据组的叶片\n2 2\n氮含量估测均方根误差为0.65 g/m ，而第二数据组的误差为0.57 g/m ，略低于第一数据组。同时可以得出，这两个数据组的叶片氮含量估测误差均较小，其相对均方根误差均在\n20%以内。\n[0049] S4，将两组半离散最优光谱集合交叉代入已进行过初步反演的反演模型进行交叉验证，得到各自的第二叶片氮含量估测值，计算两个第二叶片氮含量估测值与对应的叶片氮含量实测值之间的第二均方根误差，若第二均方根误差处于预设值内，输出反演模型。\n[0050] 在本实施例中，步骤S3所获取的对应第一数据组和第二数据组的半离散最优光谱集合交叉输入到已进行过初步反演的反演模型进行交叉验证。以第一数据组为例，亦即利用第一数据组得到的半离散最优光谱集合用于第二数据组在步骤S3中进行过初步反演的反演模型，反演得到第二叶片氮含量估测值，并将第一数据组经交叉反演所得的第二叶片氮含量估测值与第二数据组测得的叶片氮含量实测值计算第二均方根误差，通过第二均方根误差衡量反演模型的反演精度和反演策略的通用性。\n[0051] 图4a,4b表明，通过辐射传输模型反演交叉验证得到，利用第二数据组的半离散最\n2\n优光谱集合用于第一数据组，叶片氮含量估测均方根误差为0.61 g/m ，而利用第一数据组\n2\n的半离散最优光谱集合用于第二数据组，估测误差为0.59 g/m 。由此可见，反演策略交叉验证得到的叶片氮含量估测均方根误差类似于利用各自数据集的反演半离散最优光谱集合得到的叶片氮含量估测均方根误差，也即辐射传输模型反演策略具有极高的通用性和可推广性。\n[0052] 需要说明的是，上述的第一数据组和第二数据组、第一叶片氮含量估测值和第二叶片氮含量估测值，第一均方根误差和第二均方根误差中所用的“第一”、“第二”等仅用于表明数据的属性。如，第一数据组代表来自中国的叶片样本所测得的光谱反射率实测值和叶片氮含量实测值，第二数据组代表来自美国的叶片样本所测得的光谱反射率实测值和叶片氮含量实测值；第一叶片氮含量估测值代表初步反演中第一数据组和第二数据组分别计算得到的叶片氮含量估测值；第二叶片氮含量估测值代表两组半离散最优光谱集合在交叉反演中针对第二数据组和第一数据组分别计算得到的叶片氮含量估测值；第一均方根误差代表第一数据组通过初步反演得到的第一叶片氮含量估测值与第一数据组中的叶片氮含量实测值进行均方根误差，也同时代表第二数据组通过初步反演得到第一叶片氮含量估测值与第二数据组中的叶片氮含量实测值进行均方根误差；第二均方根误差代表第一数据组通过交叉反演得到的第二叶片氮含量估测值与第二数据组中的叶片氮含量实测值进行均方根误差，也同时代表第二数据组通过交叉反演得到的第一叶片氮含量估测值与第一数据组中的叶片氮含量实测值进行均方根误差。\n[0053] 上述实施例只是为了说明本发明的技术构思及特点，其目的是在于让本领域内的普通技术人员能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡是根据本发明内容的实质所做出的等效的变化或修饰，都应涵盖在本发明的保护范围内。