北冰洋固有光学特性遥感反演研究

宫丽娇，黄 珏*，吴 铭

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590；2.浙江大学 海洋学院，浙江 舟山 316021)

0 引言

水体固有光学特性(Inherent Optical Properties，IOPs)仅受水体自身组分影响，不随光照条件变化而变化，是生物光学模型的重要输入参数，对水色遥感的研究具有重要意义。水体总吸收系数(a)是水体固有光学特性的组成成分之一，由浮游植物吸收系数(aph)、有色溶解和碎屑颗粒物吸收系数(adg)以及纯水吸收系数(aw)组成。了解总吸收系数的各组成成分，有助于更好地认识生态系统与生物地球化学相互作用的问题，值得深入研究[1]。aph是浮游植物生物量的指标，可纳入初级生产模型[2-3]；由于光谱相似，非浮游植物颗粒和有色溶解有机物的吸收系数通常相加并记为有色溶解和碎屑颗粒物吸收系数。有色溶解和碎屑颗粒物光化学活性显著，在地球化学和碳循环以及海洋中的其他过程扮演着重要的角色[4-6]。

国内外学者对不同地区水体光学特性反演模型进行了诸多的研究，其中，半分析模型是应用最广泛的模型[7-10]。Hoge和Lyon[11-13]开发了线性矩阵求逆模型(The Linear Matrix Inversion Algorithm，LMI)，后经Boss和Roesler[14]修改；Garver和Siege[15]开发了半分析海洋水色模型(Garver-Siegel-Maritorena，GSM)，Maritorena等[16]开发了模拟退火程序，改善GSM模型性能；Lee等[6]基于辐射传输理论，开发了Quasi-Analytical Algorithm(QAA)来反演水体固有光学特性，并通过解析函数评估模型反演固有光学特性的可行性[17]。3种反演模型得到了广泛的验证与应用[1,5,15-16,18-23]。

随着全球气候变暖，北极生态系统和生物地球化学循环正在不断变化[24-25]，研究北冰洋IOPs的需求日益增加。国际海洋水色协调小组[15]在2006年利用实测和模拟光学数据集对各种反演模型的性能进行了比较，但是目前的反演模型通常是基于低纬度的海洋数据集开发出来的，北极地区气候严寒，北冰洋常年被冰雪覆盖，仅在6—10月部分海冰融化可通航，恶劣的地理条件使得北冰洋的水色遥感研究困难重重。现有研究证明在大陆架区域内的北冰洋，生物光学特性通常很复杂，与低纬度地区观察到的关系大不相同[26-32]。

与楚科奇海和波弗特海的实测数据相比，Wang等[27]的研究表明GSM模型反演的adg(443)误差较大；Matsuoka等[30]根据北冰洋的IOPs优化GSM模型，基于实测数据的波谱形状，将数据集分为2类，优化模型反演的adg(443)的误差分别为0.35和0.50，但是基于有限数据分析得到的结论具有偶然性。Zheng等[31]发现，QAA_V5模型在低纬度地区和北极地区，均低估了浮游植物的吸收，特别是在北极地区，模型反演的aph出现大量负值。最新形式的QAA_V6模型，根据遥感反射率Rrs(670)将水体分为2类并选用不同的参考波段进行计算，可提高模型在复杂水体中的性能，但目前还未应用于北冰洋。

本文基于北冰洋3次巡航取样的数据集，比较LMI、GSM和QAA_V6模型在高纬度水域反演IOPs的性能，并对误差原因进行分析，结合模型特性，进行区域参数优化以适用于北冰洋，为北冰洋光学遥感提供更全面的参考。

1 数据与方法

1.1 研究区域

北冰洋是以北极点为中心的水域，其海岸线十分曲折，形成了许多浅而宽的边缘海，其中，波弗特海面积为47.60万平方千米，楚科奇海为59.50万平方千米，均在70°N以上，气温低，海面几乎全年被冰雪覆盖。楚科奇海通过白令海峡与营养丰富的太平洋连通，具有较高的生产力；波弗特海受麦肯齐河影响，其养分利用率低，但是在阿蒙森湾入口处形成的冰穴延长了浮游植物的生长季节，因而在波弗特海东部观察到相对较高的初级生产力[11]。

2009年8月，MALINA探险队在波弗特海巡航取样，2010年6月和2011年7月，作为美国国家航空航天局项目“ICESCAPE”的一部分，在楚科奇海进行巡航取样。3次采样共83个站点，除去几个观测异常的站点，本文使用了77个站点的数据，站位分布如图1所示。使用的数据集同时测定了水体固有光学特性和表观光学特性，是迄今为止在北冰洋进行的较为全面的光学测量结果，这对全面评估模型的性能尤为重要。

图1 北极地区测量点位分布Fig.1 Distribution of measurement points in the arctic

1.2 数据获取

1.2.1 遥感反射率测量

Rrs(λ)被定义为上行辐照度Lu(λ，z= 0+)≡Lw(λ)和下行辐照度Ed(λ，z= 0+)≡Es(λ)的比值，其中，Lw(λ)为离水辐亮度，Lw(λ)为水面的入射辐照度。在ICESCAPE航行中，使用剖面反射率辐射仪(Biosperical)测量了包括412，443，490，510，555和665 nm在内的18个光谱波段的垂直辐射剖面。在MALINA航行中，使用便携式水体光学剖面测量系统C-OPS(Biosperical Instruments Inc.)测量了包括412，443，490，510，555和670 nm在内的18个光谱波段。开阔海域在波长为555 nm或者更长处有峰值，而沿海水域则在412～490 nm处有峰值，在这2种数据中，各随机选取4个光谱数据绘图，如图2所示。

图2 遥感反射率光谱Fig.2 Spectra of remote-sensing reflectance

1.2.2 吸收系数和散射系数测量

吸收系数ap和ad是通过测量放置在Perkin Elmer Lambda 35紫外/可见分光光度计积分球内的样品过滤器的光谱吸光度获得的，详细步骤参考文献[33-34]。将样品过滤器放在样品光束检测器的一侧，在0°，90°和180°测量3次，得到ap和ad三次扫描的平均值，ag通过UltraPath仪器测量，其中aph=ap-ad，adg=ag+ad。

水体的后向散射系数bb是用垂直剖面模式放置的原位仪器测量的。在2次ICESCAPE巡航中，bb由2套Hydract-6测量仪器分别进行测量。

1.3 实验方法

Hoge和Lyon[11-13]开发的LMI模型利用412，443和555 nm三个波长的遥感反射率，构建一个由aph，adg和bbp组成的线性矩阵方程，导出IOPs。LMI模型通过IOPs的共同变化来限制误差，光谱斜率S的变化范围在0.01～0.02，幂律指数Y的变化范围在0～2，适用于全局优化。

Garver和Siegel[15]基于遥感反射率与a和bbp之间的二次关系开发了GSM模型，引入了叶绿素a比吸收系数a*ph(λ)，针对全球海洋叶绿素浓度反演，确定了GSM模型的几个重要参数：S=0.020 6，Y=1.037 0，在443～555 nm处，a*ph(λ)分别为0.006 65，0.055 82，0.020 55，0.019 10和0.010 15；Maritorena等[30]采用全局最小化技术，改进Garver和Siegel的模拟退火模型。GSM模型需要满足3个条件才能运行：0

Lee等[6]开发的QAA模型，基于辐射传输方程的理论分析和数值模拟，导出遥感反射率与IOPs之间的关系。QAA模型主要分为2部分：第1部分模型反演得到总吸收系数a和后向散射系数bbp；第2部分是根据解析方程将a分解为aph和adg。QAA模型已更新了好几个版本，本次研究采用最新版本的QAA_V6，步骤如下所示。

将水面遥感反射率Rrs(λ)转换为水面下的遥感反射率rrs(λ)[6]：

rrs(λ)=Rrs(λ)/[0.52+1.7Rrs(λ)]。

(1)

Gordon等[9]通过模拟得到了遥感反射率与固有光学特性的关系：

rrs(λ)=g0u(λ)+g1u(λ)，

(2)

u(λ)=bb/(a+bb)，

(3)

式中，a表示水体总吸收系数；bb表示后向散射系数；g0和g1为模型参数，分别取值0.089和0.124 5。结合式(1)和式(2)，u也可以表示为：

(4)

QAA_V6第一部分，选定合适的参考波段：

若Rrs(670)<0.001 5sr-1，选取555 nm作为参考波段，a(λ0)为：

(5)

a(λ0)=a(555)=aw(555)+10-1.146-1.1366χ-0.469 χ2,

(6)

式中，aw(555)是555 nm处纯海水的吸收系数。

若Rrs(670)≥0.001 5sr-1，选取670 nm作为参考波段，a(λ0)为：

(7)

式中，aw(670)是670 nm处纯海水的吸收系数。

bb(λ)可表示为[6]：

bb=bbw(λ)+bbp(λ)，

(8)

(9)

式中，bbw(λ)和bbp(λ)为纯海水和悬浮颗粒物的后向散射系数。

当a(λ0)，u(λ0)和bbw(λ0)已知，bbp(λ0)为：

(10)

光谱斜率η可表示为：

(11)

据此，可以计算其他波长的bbp(λ)。

将bb(λ)带入式(3)，可计算其他波段的a(λ)：

(12)

QAA_V6第二部分，据解析方程分解a为：

a(λ)=aph(λ)+adg(λ)+aw(λ)，

(13)

adg可表示为：

adg(λ)=ag(443)exp[-S(λ-443)]，

(14)

式中，光谱斜率S表示为：

(15)

ag(443)可表示为：

(16)

式中，ζ和ξ为：

(17)

ξ=exp[S(442.5-415.5)]。

(18)

据此，可以计算其他波长的adg和aph，即：

adg(λ)=ag(443)exp[-S(λ-443)]，

(19)

aph(λ)=a(λ)-adg(λ)-aw(λ)。

(20)

1.4 反演结果评估

根据确定系数(R2)、均方根误差(Root Mean Square Error，RMSE)和平均相对误差(Mean Relative Error，MRE)来评估模型性能：

(21)

(22)

(23)

2 实验结果分析

本文共77个数据，但是研究发现，LMI模型反演的3个吸收系数均有负值，约占了总数据的1/7；GSM模型受3个条件的约束，有2个无效值；QAA_V6模型反演的aph在5个波段出现多个负值，与Zheng等[31]用QAA_V5模型反演得到的结论相似，本文基于有效数据值对3个模型进行评估。

采用对数坐标，如图3所示。

图3 412，443，490，510和555 nm处LMI,GSM和QAA模型反演和实测a的散点图

图3比较了LMI、GSM和QAA_V6三个模型在412，443，490，510和555 nm五个波段处反演a的性能。3个模型反演的a基本散落在1∶1的拟合线附近，其中，QAA模型反演结果最为理想，如表1所示，在412～555 nm处，R2均大于0.80，RMSE在0.13～0.63之间波动，MRE在0.33～0.68之间波动，但是在555 nm处，图3(e)高值区域，明显高估；GSM模型严重低估a(443)，而且在a(412)处，误差较大；LMI模型反演的a在含量较高区域出现负值，且相关性随波长的增加而降低，在555 nm处相关性急剧下滑，其RMSE和MRE分别为0.03，0.13。

QAA模型反演的adg与实测数据相关性最好，如图4所示。

表1 反演a(λ)，adg(λ)和aph(λ)的统计

图4 412，443，490，510和555 nm处LMI，GSM和QAA模型反演和实测adg的散点图

在adg含量高的区域，反演结果均高估，且对比GSM和LMI，RMSE和MRE精度较低，特别是在412 nm处，RMSE高达1.03；GSM模型导出的adg表现也较好，但均表现为低估，在555 nm处，其MRE为0.76；LMI反演负值出现在含量较高的区域，除掉10个无效值后，LMI模型反演的adg，其R2均大于0.70，误差波动也较少。

相比于adg，3个模型反演的aph误差都较大。QAA模型反演的aph，其R2均低于0.60，且MRE误差大，555 nm处高达19.52；如图5所示，GSM模型反演的aph基本围绕1∶1拟合线附近，分布较为分散，R2均小于0.40，但是精度较稳定；LMI模型反演的aph，其R2在0.30左右，RMSE在0.05～0.16间波动，MRE在1.5～3.2间波动，反演结果均高于实测数据且误差较大。从表1可以看出，3个模型均在555 nm处的反演结果较差。

图5 412，443，490，510和555 nm处LMI，GSM和QAA模型反演和实测aph的散点图

综合来看，QAA模型反演结果最为理想，GSM模型次之，LMI模型在北冰洋的反演结果较不理想，而且3个模型在反演a，adg和aph时，均在555 nm处得到较差的反演结果。

3 讨论与改进

3.1 讨论

GSM模型是应用于全球大尺度的反演模型，但是应用于高纬度地区时并不准确。分析数据发现，从开阔海域到沿海海域，adg(443)的值变化很大，范围分别是0.02～0.22 m-1，0.49～1.31 m-1，约占总吸收系数的77+26%，占据了主导地位，与低纬度地区的情况大不同，这与Matsuoka等[32]的研究相一致。模型反演的a和adg在各个波段整体表现为在高值区域高估，在低值区域低估，尤其低估a(443)，但是反演的aph基本分散1∶1拟合线。而且不同于前人在低纬度水域的研究[16]，GSM模型没有显著低估aph(412)。可以看出，模型反演的a和adg的变化趋势相似，而adg的反演精度依赖于参数S和Y的精度。由于颗粒的散射和高CDOM吸收，光在开阔海域的穿透深度很浅，但在沿海区域很深。因此，需要为海洋和沿海水域分别设置S和Y，Matsuoka[30]等研究了波弗特海南部海域的光谱特征，发现开阔海域在波长为555 nm或者更长处有峰值，而沿海水域则在412～490 nm处有峰值，如图2所示，依次将77个数据划分为2簇，对海洋和沿海水域分别设置S和Y参数，进行灵敏度分析，优化GSM模型。依次将本次77个数据代入，发现反演的aph的R2明显增加，在412 nm处，R2提高为0.56；在443 nm处，R2提高为0.63；在490 nm处，R2提高为0.63；在510 nm处，R2提高为0.58；在555 nm处，R2提高为0.53，尽管R2有所提高，但是MRE误差较大，在1.06～1.96区间波动；在412～555 nm范围内，模型反演的adg，其R2稳定在0.60左右；但是反演的a(443)明显改善，其R2由0.62提高为0.96，RMSE由0.25降为0.08，MRE由0.58降为0.19。改进的GSM模型在应用于北冰洋时，尽管性能有所提高，但是整体的精度仍不高，可能是受到数据集的限制。

QAA_V6模型先根据经验公式反演获得a(λ0)和bb(λ0)，再根据半分析模型推算得到adg(λ)和aph(λ)。分析可知，在a和adg时，相关性较好，但是精度有待改进；在反演aph时，结果并不理想。这是因为经验公式通常只适用于光学特性与模型开发中使用的光学特性相似的水域，因此验证模型在北冰洋的适用性，对模型进行大量的评估，并根据结果优化模型性能是必要的。

3.2 QAA_V6模型改进

本文从77个站点中随机选取23个站点，对QAA_V6进行改进，再用其余站点进行验证。

QAA_V6模型是针对开阔海域和沿海海域提出的，根据遥感反射率的不同，选取不同的参考波段。图3中的a(555)明显高估，特别是高吸收区域。Lee等[6]发现在高吸收的浑浊水域，使用大于555 nm的参考波段精确度更高，但是，高纬度地区不同于低纬度地区，气候严寒，极昼极夜现象使得人类与植物的活动痕迹较少，Zheng等[31]基于QAA_V5反演的北冰洋的a(555)与实测数据吻合良好，在670 nm处的相关性却较差。根据张亭禄等水质分类，77个站点中一类水体约占71%，而且实测a(555)>0.24 m-1的站点约占8%，水质清澈，再加上Rrs(670)远小于Rrs(555)，过小的数值使得实测数据的检测较为困难，而实测数据的准确度是优化模型的关键。QAA_V6优化模型根据Austin和Petzold的方法，选取对水质敏感的a(443)i来估算a(555)，参数根据吸收系数进行调整[6]：

(24)

412，443，490，510和555 nm处QAA_V6和改进模型反演和实测a的散点图如图6所示。

由图6可以看出，改进后的a(555)结果明显改善，特别是在高值区域，R2由0.81提高为0.91，RMSE由0.20降为0.03，MRE由0.68降为0.13。延续QAA_V6模型的分类方法，将数据分为2簇，选取不同的参考波段。在ICESCAPE中，可以根据实测的后向散射系数直接拟合η值，共14个数据值，反演的bb(443)相关性由0.66到0.84；MALINA数据缺少实测后向散射系数，将实测的总吸收系数代入到QAA_V6模型，让η在0～3之间迭代，直到反演与实测的a在5个波段中的平均相对误差小于0.20，迭代停止，得到参考η值。当Rrs(670)<0.001 5 sr-1，η与rrs(443)/rrs(555)的R2由0.47提高为0.57；当Rrs(670)≥0.001 5 sr-1，η与rrs(443)/rrs(555)的R2提也提高为0.57。根据a(555)反演的其他4个波段的总吸收系数也紧紧围绕在1∶1拟合线附近。

QAA_V6模型将a(λ)划分为adg(λ)和aph(λ)依赖参数ζ和ξ。通过拟合光谱比率(rrs(443)/rrs(555))与实测数据aph(412)/aph(443)求得ζ，adg(412)/adg(443)求取S，进而得到ξ：

(25)

ξ=exp[S(442.5-415.5)]。

(26)

412，443，490，510和555 nm处QAA_V6和改进模型反演和实测adg的散点图如图7所示。优化模型反演的adg(443)改善明显，尤其是在高值区域，基于adg(443)反演的adg(λ)与实测数据的相关性也取得较为满意的结果，R2均大于0.85。

图6 412，443，490，510和555 nm处QAA_V6和改进模型反演和实测a的散点图

图7 412，443，490，510和555 nm处QAA_V6和改进模型反演和实测adg的散点图

QAA_V6模型反演的aph出现了大量负值，约占总数据的20%～43%，在统计分析中会被删除。本文对误差源进行敏感性分析，按照控制变量的原则，对参数进行验证，发现参数η是QAA模型反演出现负值的主要原因，与Zheng等[31]的结果相一致，模型反演的bb(λ)一般会高估短于551 nm波长的后向散射和吸收系数，并可能在北冰洋的反演中产生负值。优化后的aph负值个数约占总数据的18%～35%。412，443，490，510和555 nm处QAA_V6和改进模型反演和实测aph的散点图如图8所示。由图8可以看出，aph(555)整体下降，更接近于实测数据，且在412～510 nm处，反演结果均有所改善，R2提高为0.60左右，但是在555 nm处，R2为0.29，仍不理想。优化模型反演的a和adg与实测数据较为吻合，但是aph和实测数据之间产生了较大的偏差，很大原因是受气候影响，在北冰洋中，有色溶解和碎屑颗粒物对总吸收系数的贡献更大。

图8 412，443，490，510和555 nm处QAA_V6和改进模型反演和实测aph的散点图

QAA_V6模型优化前后3个吸收系数的评价指标如图9所示。由图9可以看出，adg(λ)的R2有所下降，但整体保持在0.85以上，aph(λ)的R2明显提高；除了adg(555)的MRE稍增大，RMSE和MRE均取得了令人满意的结果。其中，adg(412)，adg(443)，aph(555)和a(555)的误差大幅下降。

3个模型反演的3个吸收系数中，浮游植物吸收系数与实测数据差异最大，很大可能是由于模型的性能。早先的研究表明，局部生物光学模型应优先于全局模型[36]，但是优化参数，提高模型的整体性能并不简单，因为模型是高度非线性的，有大量的潜在的联系。此外，400～520 nm是浮游植物中多种辅助色素的吸收谱段，导致浮游植物吸收系数与叶绿素a的相关性降低[35]。输入数据Rrs(λ)的准确性对模型反演的结果至关重要，但是，实测数据通常包含由各种实验和环境因素引起的大量不确定性，这也是模型反演产生误差的一个不可避免的来源。

4 结束语

基于北冰洋的数据，测试LMI，GSM和QAA_V6三个反演模型在北冰洋的适用性。研究发现：① QAA_V6模型反演结果最好，GSM次之，LMI模型的反演结果出现多个负值，可能不适合北冰洋的反演；② 3个模型均在555 nm处得到较差的反演结果；③ 根据北冰洋特有的光学特性优化的GSM模型，其反演结果仍不如QAA_V6模型；④ 用555 nm代替670 nm作为参考波段，QAA_V6优化模型反演的a与实测数据的相关性更高。

精确的参数是模型优化成功的关键。本文建立了适用于北冰洋的QAA_V6优化模型，为水域反演提供了更全面的参考，但水域的光学性质变化很大，受到实测点数和采样环境的限制，未来还需要更多的北冰洋的实测数据来提高模型的性能。