一种基于遥感数据反演PAR漫衰减系数垂直变化的方法与流程

一种基于遥感数据反演par漫衰减系数垂直变化的方法
技术领域
1.本发明属于海洋信息服务业技术领域，特别涉及一种基于遥感数据反演par漫衰减系数垂直变化的方法。


背景技术：

2.光合有效辐射（photosynthetically active radiation；简称par）指的是太阳辐射中对植物光合作用有效的光谱成分，其波长范围与可见光（400-700 nm）基本重合。作为光合作用的能量来源，par的分布对植物的生长（即初级生产力、碳固定）起着至关重要的作用。海洋初级生产力约贡献了全球碳固定的一半；海洋浮游植物的光合作用不仅发生在表层，而且在整个真光层均可发生。因此，了解par在海表面分布及其垂直穿透能力，即par的三维分布，对于评估海洋初级生产力，进而厘清海洋碳循环、衡量全球气候变化具有重要意义。
3.par在海表面的分布与日地距离、纬度、大气衰减等有关，而这些影响可以通过计算得出。因此，厘清par三维分布的关键在于了解海表面的par进入水体后的垂直变化，而这受到par的漫射衰减系数（diffuse attenuation coefficient，即k
par
）的调控。根据lambert-beer定律，par的垂直分布特征可以通过以下公式表述：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（9）其中，par(0-)和par(z)分别指的是海表面和深度z处的par，代表从海表面到深度z以浅的平均k
par
。作为表观光学特性（aop）参数，k
par
不仅与固有光学特性（iop；例如吸收系数和散射系数）有关，而且与水下光场的变化有关；即使在垂直均质化的水体里，k
par
也不是一个常数，而是随着深度增加而减小。因此，了解par分布的关键在于掌握k
par
在海表面以及随深度的变化规律，也就是k
par
的三维分布特征。
4.尽管k
par
很容易通过航次使用辐照仪或光量子通量仪等测量计算得出，但航次数据缺乏连续性、大面性、等时性。卫星水色遥感则可以弥补航次观测的不足，但作为非直接观测，需要建立遥感数据与k
par
之间的联系。构建漫衰减系数k
par
的遥感反演模型是当前海洋学和卫星遥感的热点，但仍存在着诸多不足。目前对于k
par
的反演模型构建仍比较缺乏且精度不足，尤其是缺乏k
par
垂直剖面的遥感推导，无法实现k
par
的垂直变化的反演。
5.目前遥感的漫衰减系数，即，通常仅代表一个光学深度内的平均k
par
值，由于k
par
在垂直方向上有明显的变化，在透光层内变化可能达到两倍之多，仅使用一个平均值来推导k
par
的垂直分布不精确，不足以准确评估水下光场的变化，也无法实现k
par
三维空间变化的精确反演。因此，构建一种新的par漫衰减系数垂直变化的方法，对于精确评估海洋初级生产力以及深入研究全球碳循环具有非常重要的意义。
6.当前国际上基于卫星遥感数据来估算par漫衰减系数(k
par
)通常有3种方法：1)基于叶绿素a浓度（chl_a）的方法；2)基于固有光学量（iop）的方法；3)基于遥感反射率(r
rs
)的波段比方法。前两种算法需要使用间接遥感产品（例如叶绿素a浓度和固有光学特性中的吸收系数和后向反射系数）来估算k
par
，因此会增加额外的不确定度。并且，卫星遥感数据中的
叶绿素a浓度和固有光学量仅局限在表层，在分层明显的水体中，k
par
垂直变化难以进行估算。辐射传输理论（例如hydrolight
®ꢀ
和monte carlo 模拟）可以通过输入垂直剖面的叶绿素a浓度和固有光学量来拟合k
par
垂向变化，但这些方法不适合运用卫星遥感的数据来反演k
par
。
7.第三种方法，即基于遥感反射率(r
rs
)的波段比方法，被广泛应用于估算全球海洋的漫衰减系数(mueller, 2000)。但该算法通常用于反演490 nm处的下行光辐照度的漫衰减系数（即k
490
），而没有专门针对k
par
的算法。例如，对于中分辨率成像光谱仪 (modis)，遥感的表层k
490 (即)反演算法如下：（10）其中：x = log[r
rs
(488)/r
rs
(555)]为了计算遥感的表层k
par
（即），常采用以下经验公式来推导：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（11）使用上述方法估算k
par
垂直分布存在两个问题。其一，该方法也需要通过中介遥感产品来计算表层光合有效辐射漫衰减系数（），会增加额外的不确定度。其次，缺乏将和k
par
垂向变化相关联的模型。 仅代表一个光学深度以浅的平均k
par
。在一个光学深度以下，k
par
通常较小很多，可能仅为其值的一半。因此，如果未考虑k
par
的垂直变化而假定所有水深的k
par
均等于，会导致光的穿透深度被严重低估。因此，建立k
par
垂向变化的模型，对于精确量化par的水下光场尤为重要。


技术实现要素：

[0008]
为解决上述技术问题，本发明提供了一种基于遥感数据反演par漫衰减系数垂直变化的方法，以达到实现par漫衰减系数垂直变化的反演的目的，对于精确量化par的水下光场具有非常重要的意义。
[0009]
为达到上述目的，本发明的技术方案如下：一种基于遥感数据反演par漫衰减系数垂直变化的方法，包括如下步骤：步骤一，通过光谱测量设备采集现场观测数据，根据现场观测数据计算遥感反射率、各物理深度的par和各物理深度的par漫衰减系数，利用表层par漫衰减系数和遥感反射率，通过数据拟合建立表层par漫衰减系数与遥感反射率的关系模型；步骤二，根据计算的各物理深度的par漫衰减系数，结合定义的光学深度，寻找par漫衰减系数随光学深度变化的规律，通过数据拟合建立某一物理深度以浅的平均par漫衰减系数与表层par漫衰减系数的关系模型；步骤三，根据建立的两个关系模型，结合卫星遥感反射率数据实现par漫衰减系数垂直变化的反演。
[0010]
上述方案中，所述现场观测数据包括：光谱测量设备测量波段范围为400～700 nm的上行光辐亮度、下行光辐照度和海表面入射光辐照度。
[0011]
上述方案中，步骤一中，根据现场观测数据计算的方法如下：（1）各物理深度的par漫衰减系数计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)其中，为物理深度为z、波长为的par漫衰减系数，为物理深度为z、波长为的下行光辐照度，z为物理深度；（2）离水辐亮度计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)其中，为波长为的离水辐亮度，为波长为的水面之下上行光辐亮度，为菲涅尔反射系数，为波长为的海水的菲涅尔反射比；（3）各物理深度的par计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)其中，为物理深度为z的par，为物理深度为z、波长为的下行光辐照度，h为普朗克常数，c为光速；（4）遥感反射率计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)其中，为波长为的遥感反射率，为波长为的离水辐亮度，表示刚好在水面之上、波长为的下行光辐照度。
[0012]
上述方案中，所述表层par漫衰减系数与遥感反射率的关系模型如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)x=log[r
rs
(488)/r
rs
(555)]其中，a0、a1、a2、a3为常数，为表层par漫衰减系数，r
rs
(488)和r
rs
(555)分别为波长488nm和波长555nm的遥感反射率。
[0013]
上述方案中，所述步骤二的具体方法如下：在透光层范围内定义一个光学深度f，以海表面光强的百分比来表示，其对应的物理深度为zf，f = par(zf)/par(0-)
×
100%其中，分别表示海表面和物理深度zf处的par，通过公式（3）计算得到；将计算的各物理深度的par漫衰减系数依光学深度f排序，当par漫衰减系数从表面下降到物理深度zf时，比较对应物理深度zf以浅的平均par漫衰减系数与表层par漫衰减系数之间的关系，建立的关系模型如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)其中，b0、b1为系数：为系数：
其中，c0、c1、c2、d0、d1、d2为常数；为指定的光学深度，为已知量；表层par漫衰减系数通过公式（5）计算得到；根据定义，
ꢀꢀꢀꢀ
(7)因此，物理深度zf计算如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)由此，得到物理深度zf以浅的平均par漫衰减系数与表层par漫衰减系数的关系模型。
[0014]
一种基于遥感数据反演par漫衰减系数垂直变化的装置，包括存储器和处理器；所述存储器，用于存储计算机程序；所述处理器，用于当执行所述计算机程序时，实现如上所述的一种基于遥感数据反演par漫衰减系数垂直变化的方法。
[0015]
一种计算机可读存储介质，所述存储介质上存储有计算机程序，当所述计算机程序被处理器执行时，实现如上所述的一种基于遥感数据反演par漫衰减系数垂直变化的方法。
[0016]
通过上述技术方案，本发明提供的一种基于遥感数据反演par漫衰减系数垂直变化的方法具有如下有益效果：1、本发明通过现场实测的光学和海洋学数据，能够构建区域优化的k
par
的三维遥感反演新方法，该方法不止适用于西太平洋，也适用于其它海域。
[0017]
2、本发明的方法可以反演各个深度的k
par
，克服了现有技术中只能研究表层水的k
par
的缺陷。
[0018]
3、本发明提出的新方法将为k
par
的三维分布及其对海洋初级生产力的评估提供新的思路和科学认识。
附图说明
[0019]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。
[0020]
图1为本发明实施例所公开的一种基于遥感数据反演par漫衰减系数垂直变化的方法流程示意图。
[0021]
图2为本发明提出的方法的研究区域和海洋现场观测数据站位分布图。
[0022]
图3为本发明中提出的遥感反射率波段比与遥感的par漫衰减系数之间的关系图。
[0023]
图4a为表层光合有效辐射漫衰减系数（）和490 nm处下行辐照度漫衰减系数（）之间的关系图。
[0024]
图4b为 log(
ꢀ‑
0.0166)和log[rrs(488)/rrs(555)]之间的关系（数据来源于断面t1-t3实际测量的数据）。
[0025]
图5a为运用本发明计算的t1-t3断面不同深度以浅par平均漫衰减系数（）和
实际观测数据的对比图。
[0026]
图5b为运用本发明计算的物理深度zf和实际观测数据的对比图。
[0027]
图6a为运用本发明计算的的有效性验证图。
[0028]
图6b为运用本发明计算的zf的有效性验证图。
具体实施方式
[0029]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
[0030]
本发明提供了一种基于遥感数据反演par漫衰减系数垂直变化的方法，如图1所示，包括如下步骤：步骤一，通过光谱测量设备采集现场观测数据，根据现场观测数据计算遥感反射率、各物理深度的par和各物理深度的par漫衰减系数，利用表层par漫衰减系数和遥感反射率，通过数据拟合建立表层par漫衰减系数与遥感反射率的关系模型；现场观测数据（包含t1-t3断面）包括：提取使用美国satlantic公司的 hyperpro ii水下高光谱仪测量波段范围为400-700 nm的上行光辐亮度lu、下行光辐照度ed和海表面入射光辐照度es。
[0031]
根据现场观测数据计算的方法如下：（1）各物理深度的par漫衰减系数计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)其中，为物理深度为z、波长为的par漫衰减系数，为物理深度为z、波长为的下行光辐照度，z为物理深度；（2）离水辐亮度计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)其中，为波长为的离水辐亮度，为波长为的水面之下上行光辐亮度，为菲涅尔反射系数，为波长为的海水的菲涅尔反射比；（3）各物理深度的par计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)其中，为物理深度为z的par，为物理深度为z、波长为的下行光辐照度，h为普朗克常数，c为光速；（4）遥感反射率计算公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)其中，为波长为的遥感反射率，为波长为的离水辐亮度，表示刚好在水面之上、波长为的下行光辐照度。
[0032]
使用t1-t3断面现场观测的数据进行拟合，如图2所示，图2中背景图为航次期间（2017年4月14日至5月13日）的海表温（sst）的平均分布；t1、t2、t3和t4为采样断面，
○
代表
有采集par垂直剖面的站位，t1-t3断面实际观测的站位数n=19；ke、npsg、wtp分别代表黑潮延伸体、副热带环流系、西热带太平洋三个主要水团，tz-1和tz-2为过渡区域；箭头表示该区域的主要表面洋流。t1-t4断面的现场观测站位信息如表1所示。
[0033] 表1 t1-t4断面现场观测站位信息建立表层par漫衰减系数与遥感反射率的关系模型如图3所示，图3 为使用断面t1-t3实际观测的数据建立的log[r
rs
(488)/r
rs
(555)]和log()之间的关系。下实线是依据本发明得出的最好的拟合线；虚线是均方根误差为
±
0.034的边界线；上实线是使用公式(10)和（11）得到的两者的关系。本发明得到的拟合曲线如下：
ꢀꢀꢀꢀꢀꢀ
(5)拟合得到的上述关系模型的相关系数r
2 = 0.952。
[0034]
x=log[r
rs
(488)/r
rs
(555)]其中，为表层par漫衰减系数，r
rs
(488)和r
rs
(555)分别为波长488nm和波长555nm的遥感反射率。
[0035]
上述公式中得到log( ) 的均方根误差是
±
0.036，相应的的不确定度为
±
8%。相比于使用(10)和(11)组合公式，得到的log()的均方根误差为
±
0.075，相应的
不确定度为
±
18%。因此，本发明中提出的反演方法中，通过遥感图像中的遥感反射率r
rs
波段比和表层par漫衰减系数直接建立关系，不确定度减少了一半。
[0036]
从图2很明显看出，使用上述公式(10)和(11)组合导致表层par漫衰减系数的高估。这种偏差主要来源于公式(11)，如图4a，通过 （范围：0.03-0.10 m-1
）计算，明显看出高估 24%(图4a)。通过公式(10)计算的 不确定度是比较小的，计算log( )得到的均方根误差为
±
0.053。从公式(10)的实际观测数据回归曲线中也可以看出没有明显的弯曲(图4b)。图4a和图4b中的实线分别是公式(11)和(10)的回归线。图4b虚线表示公式(10)中回归线的均方根误差范围(
±
0.055)。
[0037]
步骤二，根据计算的各物理深度的par漫衰减系数，结合定义的光学深度，寻找par漫衰减系数随光学深度变化的规律，通过数据拟合建立某一物理深度以浅的平均par漫衰减系数与表层par漫衰减系数的关系模型。
[0038]
具体方法如下：在透光层范围内定义一个光学深度f，以海表面光强的百分比来表示，其对应的物理深度为zf，f = par(zf)/par(0-)
×
100%其中，分别表示海表面和物理深度zf处的par，可通过公式（3）计算得到；将计算的各物理深度的par漫衰减系数依光学深度f排序，当par漫衰减系数从表面下降到物理深度zf时，比较对应物理深度zf以浅的平均par漫衰减系数与表层par漫衰减系数之间的关系，建立的关系模型如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)其中，b0、b1为系数：其中， 为指定的光学深度，为已知量；表层par漫衰减系数通过公式（5）计算得到；根据定义，
ꢀꢀꢀꢀ
(7)因此，物理深度zf计算如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)由此，得到物理深度zf以浅的平均par漫衰减系数与表层par漫衰减系数的关系模型。
[0039]
步骤三，根据建立的两个关系模型，结合卫星遥感反射率数据实现par漫衰减系数垂直变化的反演。
[0040]
具体为，联合公式(5)和(6)，给定光学深度f，其物理深度zf以浅的平均par漫衰减系数能够通过遥感反射率r
rs
计算得到。对应的物理深度zf可以通过公式(8)进一步计算得到。由此，得到物理深度为zf时的平均par漫衰减系数与表层par漫衰减系数的关系，实现par漫衰减系数垂直变化的反演。
[0041]
使用本发明的方法计算得到的log()的均方根误差(rmse)约为
±
0.051（范围：
±
0.039
‑‑±
0.061），对应f变化从50% 到1%。相应的不确定度为
±
12%。计算得到的和zf如图5a、图5b所示。图5a和图5b中，实际观测数据来源于t1-t3断面的所有站位，图5a和图5b不同形状的散点分别代表光强为表层值的30% (
○
)，10% (
∆
)，3% (
◊
) and 1% (
□
)。实线为1:1线，虚线为
±
12%偏离线。上述结果表明，只要通过卫星遥感图像给定足够准确的遥感反射率r
rs
。在透光层任何光强度区域，利用卫星遥感图像数据可以精确计算空间垂直剖面的漫衰减系数及对应的漫衰减系数穿透深度，即物理深度zf，计算的不确定度仅为
±
12%。
[0042]
其次，为了进一步验证本发明提出来的新算法的区域有效性，将本发明基于遥感反射率数据r
rs
(λ)计算得出的和zf与真实观测的数据进行对比（选取区域为延赤道区域的t4断面和美国夏威夷大学aloha 站位的历史观测数据）。验证结果表明：当光强变化从f=50%到1%，本发明的方法计算得到的log()的均方根误差(rmse)为
±
0.043。相应的
ꢀꢀ
(图6a) 和 z
f (图6b)的不确定度为
±
10%。图6a和图6b 类似于图5a和图5b，数据来自美国夏威夷大学aloha站位和 t4断面站位观测数据。结果表明，本发明的方法基于遥感数据计算的和zf与实际观测的数据具有良好的一致性。图6a、图6b中的实线为1：1线，虚线为
±
12%偏离线。
[0043]
本领域内的技术人员应明白，本发明的实施例可提供为方法、装置或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、cd-rom、光学存储器等）上实施的计算机程序产品的形式。
[0044]
本发明是参照本发明实施例的方法、设备（系统）、和计算机程序产品的流程图和／或方框图来描述的。应理解可由计算机程序指令实现流程图和／或方框图中的每一流程和／或方框，以及流程图和／或方框图中的流程和／或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和／或方框图一个方框或多个方框中指定的功能的系统。
[0045]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程和／或方框图一个方框或多个方框中指定的功能。
[0046]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和／或方框图一个方框或多个方框中指定的功能的步骤。
[0047]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。