1、引　言

地球表面约3/4被海洋覆盖，海洋在大气辐射传输、水循环和全球气候变化中起着至关重要的作用(陈兴峰 等，2011)。在海洋遥感中卫星接收的辐射信息主要包括大气散射、波浪镜面反射、海面白沫辐射和离水辐射等，其获取和处理方式直接影响海洋辐射信息探测的精度。特别是太阳光入射海表特定区域因波浪菲涅尔反射形成海洋耀光，对载荷成像质量影响巨大。该强反射区域通常称为海洋耀光区，具有显著的偏振辐射特性。研究海洋耀光辐射及偏振辐射特性与这些参数的关系，建立海洋耀光的检测方法，可以为卫星遥感影像的后续处理和高质量数据生产提供重要技术支持(毛志华 等，1996)。

研究光学载荷接收到的耀光区辐射信息主要考虑来自水体贡献(Cox和Munk，1954a)。为进一步研究耀光区的辐射特性，作者在合肥董铺水库进行了水面耀光光谱实验研究(图1)，可以看出耀光区的形状、大小和辐射特性与水面风速风向、太阳及观测几何等相关。通过对比可见光—近红外波段耀光区和非耀光区的反射率差异(图2)，发现近红外(NIR)波段对耀光区辐射信息更加敏感，可有效用于耀光检测。通常大气层顶的偏振辐射信息受大气分子、云和气溶胶的影响，尤其在近红外波段更为突出。研究大气分子、云和气溶胶对偏振敏感性不同的特性(Li 等，2006；顾行发 等，2015；伽丽丽 等，2016)，使偏振遥感在海洋和大气研究中更具优势。法国科学家从1996年—2004年先后发射3颗POLDER(Polarization and Directionality of Earth’s Reflectances)偏振载荷用于对地观测研究，在海气系统研究中提出了一些新方法。Bréon和Henriot(2006)在不考虑大气分子和气溶胶影响的情况下，利用NIR波段的POLDER数据研究了海洋耀光区辐射特性与海表波浪之间的关系，建立了海浪斜率分布的辐射特性模型。Goloub等人(1995, 1997)利用NIR波段偏振辐射信息进行了海洋上空云微物理特性研究和气溶胶参数反演(Bréon 等，1995；Deuzé 等，2001)。美国国家航空航天局为提高全球气溶胶特性探测精度，研发了APS(Aerosol Polarimetry Sensor)载荷。中科院安光所科研团队研制了星载多角度偏振探测仪 (DPC)，高分五号卫星(GF-5)用于全球尺度的对地偏振遥感研究(顾行发 等，2015)。DPC在轨运行时，获取的观测数据需要经过一系列处理才能进行云和气溶胶特性反演，海洋耀光检测是数据处理的重要环节。目前国内外学者针对海洋耀光的检测提出了一系列方法，主要有Cox-Munk模型法(Cox和Munk，1954b)和经验值法。现有大多数载荷参考MODIS(Moderate Resolution Imaging Spectroradiometer) 40°耀光角经验值去除耀光(Remer 等，2005)，法国的POLDER载荷联合粗糙海表偏振模型与经验临界值对耀光进行检测，其业务化处理采用30°经验值判别(Toubbe 等，1999)。以上方法虽然能一定程度检测出海洋耀光，但是采用上述经验临界值用于DPC载荷的海洋耀光检测并不十分准确，会导致一些像元漏检或误判为耀光像元，反演过程中像元的有效利用率降低。因此如何动态获取耀光角判别临界值是实现海洋耀光动态检测的关键，本文将针对此研究难点开展研究。

考虑DPC的主要技术指标参照法国POLDER载荷，本文以PARASOL(Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar)卫星搭载的POLDER3载荷一级数据产品为研究对象，结合海气辐射传输相关理论，考虑大气分子和气溶胶的影响，选取近红外(865 nm)偏振通道数据，提出基于NIR偏振辐射卫星数据的海洋耀光动态检测(OGDD)模型。该模型首先对海洋耀光的相关参量(耀光角和偏振辐射特性)进行计算，再利用耀光区偏振辐射特性规则表现的OGDD算法，获取准确的耀光角判别临界值，最后进行算法实验验证和结果分析。

2、海气辐射传输相关理论

在平面平行大气近似下，假设大气的光学和微物理特性水平均一不变(潘德炉 等，1997)。晴朗天气下海洋上空太阳辐射入射大气层，一部分经气体分子及气溶胶散射和吸收形成上行辐射，一部分下行到达海表，经海表吸收和镜面反射重新返回大气，另一部分透过海表，在海水次表面受叶绿素、悬浮颗粒和水溶物质等后向散射经海面折射再次返回大气，返回大气的辐射再次透过大气到达大气层顶(TOA)。可得卫星载荷入瞳处接收到海洋耀光区的总辐射量计算公式

${L_{\rm{tot}}} = t({L_{\rm{ray}}} + {L_{\rm{aer}}} + {L_{\rm{wat}}} + {L_{\rm{f}}}{\rm{ + }}{L_{\rm{glit}}})$

(1)

式中， ${L_{\rm{tot}}}$ 为卫星载荷接收的总辐射； ${L_{\rm{ray}}}$ 为大气分子散射辐射； ${L_{\rm{aer}}}$ 为气溶胶散射辐射； ${L_{\rm{wat}}}$ 为离水辐射； ${L_{\rm{f}}}$ 为白沫辐射； ${L_{\rm{glit}}}$ 为耀光反射辐射； $t$ 为海表和载荷之间大气透过率。

通常海洋耀光区的海气辐射传输按能量贡献大小分为耀光辐射、大气程辐射、离水辐射和白沫辐射等，其中后3者称为非耀光辐射。

    (2.1) 非耀光辐射

大气程辐射由大气分子散射和气溶胶散射构成。考虑单次散射，其偏振光大气散射近似表示为(Toubbe 等，1999)

$\rho _{\rm{pol}}^{\rm{X}} = \frac{{{\text{π}} L_{\rm{pol}}^{\rm{X}}}}{{E\cos {\theta _s}}} = \frac{{{\text{π}} \cdot {\tau ^{\rm{X}}}q_\lambda ^{\rm{X}}(\Theta)}}{{4E\cos {\theta _s}\cos {\theta _v}}}$

(2)

式中， $L_{\rm{pol}}^{\rm{X}}$ 为瑞利或气溶胶偏振辐射(X = ray为瑞利散射；X = aer为气溶胶散射)； $E$ 为太阳光辐照度； ${\tau ^{\rm{X}}}$ 为分子或气溶胶光学厚度； $q_{\textit{λ}} ^{\rm{X}}(\Theta)$ 为偏振相函数； ${\textit{λ}} $ 为波长； ${\theta _v}$ 为观测天顶角； ${\theta _s}$ 为太阳天顶角。

离水辐射包括可见光离水辐射和荧光离水辐射。即一部分太阳光透过海表向下传输时受水色因子(叶绿素、悬浮颗粒等)后向散射经水面折射出射，另一部分受叶绿素影响产生荧光出射海表(潘德炉 等，1997)。

白沫辐射由海洋表面泡沫覆盖比率决定。海洋白沫贡献同样采用各向同性朗伯反射 ${\rho _f}$ ，泡沫内在反射率可由海表风速 $v$ 经验公式(不依赖波长)计算得到(Koepke，1984；Gordon和Wang，1994)。

    (2.2) 耀光辐射

由上述海气辐射传输理论可知，相对于海洋耀光区卫星载荷接收到的总辐射能量，离水辐射和白沫辐射贡献较小，因此主要考虑贡献较大的粗糙海表太阳反射辐射。海面受海风的影响变得粗糙不平，可视其为微面元结构。微面元反射辐射利用理想菲涅尔反射相关理论计算，而粗糙海表反射辐射则利用Cox-Munk海表偏振反射模型来获取(Cox和Munk, 1955, 1956)。

         2.2.1. 菲涅尔反射

考虑理想平静海面，太阳光入射海表将发生折射和反射，反射光的反射率和偏振度可由Fresnel公式计算。根据Fresnel定理，设 ${A}$ 为入射辐射振幅，可分解为平行和垂直入射面的 ${{A}_\parallel }$ 、 ${{A}_ \bot }$ 两个分量，则反射辐射的两个分量为

$\left\{ \begin{aligned} & {{R}_\parallel } = \frac{{{m^2} \cdot \cos {\theta _i} - \sqrt {{m^2} - {{\sin }^2}{\theta _i}} }}{{{m^2} \cdot \cos {\theta _i} + \sqrt {{m^2} - {{\sin }^2}{\theta _i}} }}{{A}_\parallel } \\ & {{R}_ \bot } = \frac{{\cos {\theta _i} - \sqrt {{m^2} - {{\sin }^2}{\theta _i}} }}{{\cos {\theta _i} + \sqrt {{m^2} - {{\sin }^2}{\theta _i}} }}{{A}_ \bot } \end{aligned} \right.$

(3)

式中， ${\theta _i}$ 为太阳光入射角或反射角， $m = {{{n_2}}/{{n_1}}}$ ， ${n_2}$ 和 ${n_1}$ 分别为海水和空气折射率。可得Fresnel反射率 $\, {\rho _{fr}}$ 和偏振反射率 $\, \rho _{fr}^P$

$\left\{ \begin{gathered} {\rho _{fr}} = {{|r_\parallel ^2 + r_ \bot ^2|}/2} \\ \rho _{fr}^P = {{|r_\parallel ^2 - r_ \bot ^2|}/2} \\ \end{gathered} \right.$

(4)

考虑入射海表的太阳光线部分被海水吸收 $(n = {n_r} - i{n_i})$ ，其偏振反射率系数 $\, \rho _{fr}^P(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v})$ 可由式(5)计算得到(Born和Wolf，1975；Liou，2002)

$\begin{aligned} & \rho _{fr}^P(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v}) = \\ & \frac{1}{2}\left[ { - \frac{{{{(\cos {\theta _i} - u)}^2} + {v^2}}}{{{{(\cos {\theta _i} + u)}^2} + {v^2}}} + } \right. \\ & \left. {\frac{{{{[(n_r^2 - n_i^2)\cos {\theta _i} - u]}^2} + {{(2{n_r}{n_i}\cos {\theta _i} - v)}^2}}}{{{{[(n_r^2 - n_i^2)\cos {\theta _i} + u]}^2} + {{(2{n_r}{n_i}\cos {\theta _i} + v)}^2}}}} \right] \end{aligned} $

(5)

式中， ${u^2} = \displaystyle\frac{1}{2}\left\{ \left| {n_r^2 - n_i^2 - {{\sin }^2}{\theta _i}} \right| + \sqrt {{{(n_r^2 - n_i^2 - {{\sin }^2}{\theta _i})}^2} + 4 n_r^2 n_i^2} \right\} $ , ${v^2} \!=\! \displaystyle\frac{1}{2}\left\{ - \left| {n_r^2 - n_i^2 - {{\sin }^2}{\theta _i}} \right| +\!\! \sqrt {{{(n_r^2 - n_i^2 - {{\sin }^2}{\theta _i})}^2} + 4 n_r^2 n_i^2}\right \} $ , $\cos {\theta _i} = \sqrt {\displaystyle\frac{1}{2}[1 + \cos {\theta _s}\cos {\theta _v} + \sin {\theta _s}\sin {\theta _v}\sin ({\varphi _s} - {\varphi _v})]} $ , $\sin {\theta _i} = \sqrt {\displaystyle\frac{1}{2}[1 - \cos {\theta _s}\cos {\theta _v} - \sin {\theta _s}\sin {\theta _v}\sin ({\varphi _s} - {\varphi _v})]} $ , ${n_r}$ 表示海水折射指数， ${n_i}$ 表示海水消光系数， ${\theta _i}$ 为光线入射海浪微面元的入射角。

         2.2.2. 粗糙海表反射

实际上，因受海风的影响海表变得粗糙不平，理想菲涅尔反射不再适用。Cox和Munk通过航空影像研究太阳光入射粗糙海表的反射辐射，结合海表风速建立了微面元斜率各向同性和异性的高斯分布粗糙海表模型。该模型将粗糙海表看作由不同倾斜斜率的微面元组成，太阳光在单个微面元呈镜面反射。考虑这样一个坐标系(P;X, Y, Z)，P为观测点，Z为海拔，PY为太阳方向，PX为垂直太阳主平面方向。面元倾斜度由 ${Z_x}$ 和 ${Z_y}$ 组成(Bréon，1993)

$\left\{ \begin{gathered} {Z_x} = \partial Z/\partial X = \sin (\alpha)\tan (\beta) \\ {Z_y} = \partial Z/\partial Y = \cos (\alpha)\tan (\beta) \\ \end{gathered} \right.$

(6)

$\beta {\rm{ = }}\tan{^{ - 1}}{(Z_x^2 + Z_y^2)^{{1/2}}}$

(7)

式中， $\alpha $ 、 $\beta $ 分别为面元法线方向方位角和倾斜角。由太阳几何(天顶角 ${\theta _s}$ , 方位角 ${\varphi _s}$ )和观测几何(天顶角 ${\theta _v}$ ，方位角 ${\varphi _v}$ )可知，倾斜面元满足条件

$\left\{ \begin{aligned} & {Z_x} = - \sin {\theta _v}\sin ({\varphi _s} - {\varphi _v})/(\cos {\theta _s} + \cos {\theta _v}) \\ & {Z_y} = (\sin {\theta _s} + \sin {\theta _v}\cos ({\varphi _s} - {\varphi _v}))/(\cos {\theta _s} + \cos {\theta _v}) \end{aligned} \right.$

(8)

$\left\{ \begin{aligned} & {{Z'}_x} = \cos ({\varphi _s} - {\varphi _w}) \cdot {Z_x} + \sin ({\varphi _s} - {\varphi _w}) \cdot {Z_y} \\ & {{Z'}_y} = - \sin ({\varphi _s} - {\varphi _w}) \cdot {Z_x} + \cos ({\varphi _s} - {\varphi _w}) \cdot {Z_y} \end{aligned} \right.$

(9)

式中， ${\varphi _s} - {\varphi _w}$ 为坐标系从太阳方位角 ${\varphi _s}$ 旋转到风向方位角 ${\varphi _w}$ 的旋转角度。

倾斜面元的概率分布用Gram-Charlier级数展开(Cox和Munk, 1954a, b，1955，1956)

$\begin{aligned} P({{Z'}_x},{{Z'}_y}) = & {(2{\text{π}} {{\sigma '}_x}{{\sigma '}_y})^{ - 1}}\exp (\frac{{ - {\xi ^2} - {\eta ^2}}}{2})\{ 1 - \\ & \frac{1}{2}{c_{21}}({\xi ^2} - 1)\eta - {\textstyle{1 \over 6}}{c_{03}}({\eta ^3} - 3\eta ) + \\ & \frac{1}{{24}}{c_{40}}({\xi ^4} - 6{\xi ^2} + 3) + \\ & \frac{1}{{24}}{c_{04}}({\eta ^4} - 6{\eta ^2} + 3) + \\ & \frac{1}{4}{c_{22}}({\xi ^2} - 1)({\eta ^2} - 1)\} \end{aligned}$

(10)

式中， $\xi = {Z'_x}/{\sigma '_x}$ ， $\eta = {Z'_y}/{\sigma '_y}$ ， ${\sigma '_x}$ ， ${\sigma '_y}$ 分别为 ${Z'_x}$ 和 ${Z'_y}$ 的均方根，斜偏系数 ${c_{21}}$ 、 ${c_{03}}$ 和峰度系数 ${c_{40}}$ 、 ${c_{04}}$ 、 ${c_{22}}$ 都用风速的线性函数来表示。

可得粗糙海表的耀光反射率 ${R_{\rm{glit}}}$ 计算公式

${R_{\rm{glit}}} = \frac{{{\text{π}} {L_{\rm{glit}}}}}{{{E_0}\cos {\theta _s}}} = \frac{{{\text{π}}{\rho _{fr}}(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v})}}{{4\cos {\theta _s}\cos {\theta _v}{{\cos }^4}\beta }}P({Z'_x}, {Z'_y})$

(11)

式中， $\,{\rho _{fr}}(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v})$ 表示Fresnel反射系数( $n$ 为海水复折射指数)，同样可以将Fresnel偏振反射率 $\rho _{fr}^P(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v})$ 代入式(11)计算出粗糙海表的耀光偏振反射率 $R_{\rm{glit}}^P$ 。

由上可知，TOA海洋耀光辐射不仅需要考虑太阳辐射入射到海洋耀光区经海水吸收的粗糙海表反射，还要兼顾上行辐射透过大气层时的分子散射和气溶胶散射，即卫星载荷入瞳处的总辐射只考虑海面耀光辐射和大气程辐射。从式(1)和式(11)可推出海洋耀光区的TOA总偏振反射率 $R_{tot}^P$ 计算公式

$R_{\rm{tot}}^P = t(R_{\rm{glit}}^P{\rm{ + }}R_{\rm{ray}}^P + R_{\rm{aer}}^P)$

(12)

式中， $R_{\rm{ray}}^P$ 、 $R_{\rm{aer}}^P$ 分别为大气分子和气溶胶贡献(与波长、气压或海拔相关)。

3、海洋耀光动态检测模型

海洋耀光区TOA总偏振辐射不仅受海表风速风向及大气的影响，而且海表微面元倾斜形成的耀光角与太阳及观测几何相关。因此，检测海洋耀光首先必须获取耀光角值，再结合海气辐射传输理论推导海洋耀光相关参量计算方法，对获取的各参量值进行动态分析。

    (3.1) 海洋耀光相关参量计算          3.1.1. 耀光角计算

太阳光入射粗糙海表微面元时发生折射和反射，定义入射光和反射光之间的夹角(考虑光线方向)为散射角 ${\varTheta _{\rm{sca}}}$ ，由球面几何学可得

${\varTheta _{\rm{sca}}}{\rm{ = }}{\text{π}} {\rm{ - co}}{{\rm{s}}^{ - 1}}{\rm{(cos}}{\theta _s}{\rm{cos}}{\theta _v}{\rm{ + sin}}{\theta _s}{\rm{sin}}{\theta _v}{\rm{cos}}\phi)$

(13)

式中，相对方位角 $\phi = {\varphi _s} - {\varphi _v}$ 。同理太阳光入射海表倾斜微面元形成耀光角 ${\varTheta _{\rm{glit}}}$ ，可推出计算公式(Cox和Munk，1954a)

${\varTheta _{\rm{glit}}}{\rm{ = co}}{{\rm{s}}^{ - 1}}{\rm{(cos}}{\theta _s}{\rm{cos}}{\theta _v}{\rm{ - sin}}{\theta _s}{\rm{sin}}{\theta _v}{\rm{cos}}\phi)$

(14)

利用耀光角计算公式可对卫星遥感数据进行耀光检测。MODIS海洋耀光的检测主要以40°耀光角经验值判别，满足 ${\varTheta _{\rm{glit}}}$ 小于40°的海洋像元被标记为耀光(Remer 等，2005)。法国Toubbé为了对POLDER载荷进行在轨定标，提出了一种耀光区域检测及校正方法，并对耀光校正角 ${\theta _{\rm{glit}}}$ 和 ${\phi _{\rm{glit}}}$ (天顶和方位)做了如下条件限定(Toubbe 等，1999；Hagolle 等，2004；Harmel和Chami，2013)

$\left\{ \begin{gathered} {\theta _s} - 20^\circ < {\theta _{\rm{glit}}} < {\theta _s} + 20^\circ \\ - 40^\circ < {\phi _{\rm{glit}}} < 40^\circ \\ \end{gathered} \right.$

(15)

参考上述海洋耀光检测方法，本文OGDD模型的耀光角参量采用式(14)计算获取，耀光角临界值以[20°, 40°]为参考进行动态调整。

         3.1.2. 海洋耀光偏振辐射特性计算

(1) 耀光像元NIR波段大气层顶反射率和偏振反射率计算。考虑晴空条件下，TOA海洋耀光偏振辐射特性主要受海表耀光辐射、大气分子辐射、气溶胶辐射(不考虑白沫和离水辐射)影响。根据式(2)、式(11)和式(12)可得到海洋耀光NIR波段(865 nm) TOA反射率 $R_{\rm{glit}}^{TOA}$ 和偏振反射率 $PR_{\rm{glit}}^{TOA}$ 计算公式

$\left\{\begin{aligned} & R_{{\rm{glit}}}^{TOA} = t\left\{ {\frac{{{\text{π}} {\rho _{fr}}(n,{\theta _s},{\theta _v},{\varphi _s},{\varphi _v})}}{{4\cos {\theta _s}\cos {\theta _v}{{\cos }^4}\beta }}P({{Z'}_x},{{Z'}_y}) + } \right. \\ & \qquad\quad \left. {\frac{{{\text{π}} {\tau _r}{P_r}(\varTheta )}}{{4E\cos {\theta _s}\cos {\theta _v}}} +\frac{{{\text{π}} {\tau _a}{P_a}(\varTheta )}}{{4E\cos {\theta _s}\cos {\theta _v}}}} \right\}\\ & PR_{{\rm{glit}}}^{TOA} = t\left\{ {\frac{{{\text{π}} \rho _{fr}^P(n,{\theta _s},{\theta _v},{\varphi _s},{\varphi _v})}}{{4\cos {\theta _s}\cos {\theta _v}{{\cos }^4}\beta }}P({{Z'}_x},{{Z'}_y}) + } \right. \\ & \qquad\qquad \left. {\frac{{{\text{π}} {\tau _r}{q_r}(\varTheta )}}{{4E\cos {\theta _s}\cos {\theta _v}}} +\frac{{{\text{π}} {\tau _a}{q_a}(\varTheta )}}{{4E\cos {\theta _s}\cos {\theta _v}}}} \right\} \end{aligned} \right. $

(16)

式中， ${\tau _r}$ 和 ${\tau _a}$ 分别表示分子和气溶胶光学厚度， ${P_r}(\varTheta)$ 和 ${P_a}(\varTheta)$ 分别表示分子和气溶胶相函数， ${q_r}(\varTheta)$ 和 ${q_a}(\varTheta)$ 分别表示分子和气溶胶偏振相函数， $\varTheta $ 为光线入射方向和散射方向之间的矢量夹角。

当满足海洋像元耀光角值 ${\varTheta _{\rm{glit}}}$ 最小，且几何约束条件满足太阳天顶角等价于观测天顶角 ${\theta _s}{\rm{ = }}{\theta _v}$ ，相对方位角满足 $\varphi {\rm{ = }}{\text{π}} $ ，则将该像元判别为耀光中心 $\varTheta _{\rm glit}^{cet}$ 。由式(16)可推导出耀光中心像元的NIR波段(865 nm)TOA反射率 $cetR_{\rm{glit}}^{TOA}$ 计算公式

$\begin{aligned} cetR_{\rm{glit}}^{TOA} = & t\left\{ {\frac{{{\text{π}}{\rho _{fr}}(n, {\theta _s}, {\theta _v}, {\varphi _s}, {\varphi _v})}}{{4{{\cos }^2}{\theta _s}{{\cos }^4}\beta }}P({{Z'}_x}, {{Z'}_y})} \right. + \\ & \left. {\frac{{{\text{π}} ({\tau _r}{P_r}(\varTheta) + {\tau _a}{P_a}(\varTheta))}}{{4E{{\cos }^2}{\theta _s}}}} \right\} \end{aligned} $

(17)

同样将 $\,\rho _{fr}^P$ 代入式(17)可以计算耀光中心像元的NIR波段(865 nm)偏振反射率 $cetPR_{\rm{glit}}^{TOA}$ 。联合已知的太阳及观测几何 ${\rm{(}}{\theta _{\rm{s}}}{\rm{, }}{\theta _v}, {\varphi _s}{\rm{, }}{\varphi _v})$ 、海表风速风向 $v$ 、气溶胶光学厚度 $\tau $ 等信息，根据式(16)、式(17)计算得到NIR波段(865 nm)耀光中心像元 $pix(I, J)$ 的模拟TOA反射率 $cetR_{\rm{glit}}^{TOA}(I, J)$ 和偏振反射率 $cetPR_{\rm{glit}}^{TOA}(I, J)$ 。耀光区像元的海表风速 $v$ 从ECMWF(European Centre for Medium range Weather Forecasts)获取，ECMWF数据网格和分辨率不同于POLDER3一级产品，本文采用最邻近插值获得该像元对应经纬度位置的风速风向。气溶胶光学厚度 $\tau $ 通过设定固定步长循环获取。

(2) 海洋耀光区NIR波段偏振辐射特性均值计算。假定待检测图像海洋耀光区域为 $Apix(0, 0, 0;I, J, Z)$ ， $I, J$ 分别为图像的行数和列数， $Z$ 为观测图像层数。即POLDER3载荷对地观测成像时，依托卫星运动使其对地面同一目标进行多角度观测(观测角度数为13—15次)。获取9个波段的观测数据(6个非偏，3个偏振)，按波段分组存储，一个波段数据组由以全球网格(行×列)形式存储的多角度观测数据组成(Deschamps 等，1994)，将角度数为 $Z$ 的观测数据简称为 $Z$ 观测层。该耀光区域的NIR波段(865 nm)TOA平均反射率 $\overline R_{\rm{glit}}^{TOA}$ 和平均偏振反射率 $\overline {PR} _{\rm{glit}}^{TOA}$ ，可由如下公式计算

$\overline R_{\rm{glit}}^{TOA} = (\sum\limits_{I = 0}^M {\sum\limits_{J = 0}^N {R_{\rm{glit}}^{TOA}(I, J, Z)} })/(I \times J)$

(18)

$\overline {PR} _{\rm{glit}}^{TOA} = (\sum\limits_{I = 0}^M {\sum\limits_{J = 0}^N {PR_{\rm{glit}}^{TOA}(I, J, Z)} })/(I \times J)$

(19)

式中， $I, J, Z$ 分别表示待检像元所处的行列号和层数。

NIR波段(865 nm)耀光中心TOA反射率及偏振反射率的实测值与模拟值之间绝对偏差 $(\varDelta _{I, J, Z}^R, \varDelta _{I, J, Z}^{PR})$ 分别可以通过如下公式计算：

$\left\{ \begin{aligned} & \varDelta _{I, J, Z}^R = TR_{\rm{glit}}^{TOA}(I, J, Z) - cetR_{\rm{glit}}^{TOA}(I, J, Z) \\ & \varDelta _{I, J, Z}^{PR} = TPR_{\rm{glit}}^{TOA}(I, J, Z) - cetPR_{\rm{glit}}^{TOA}(I, J, Z) \end{aligned} \right.$

(20)

式中， $TR_{\rm{glit}}^{TOA}, TPR_{\rm{glit}}^{TOA}$ 分别为卫星载荷的实测TOA反射率和偏振反射率。

    (3.2) 本文方法原理

通过上述海洋耀光相关参量的计算方法，不仅可以准确获取晴空条件下耀光像元的NIR波段(865 nm)TOA反射率和偏振反射率，还可以对耀光区进行NIR波段偏振辐射特性均值计算。由于上述参量与耀光角之间存在一定关联，可建立它们之间关系表达式。根据球面几何条件建立耀光角临界值 $\varTheta _{\rm{glit}}^{lin}$ 与耀光区角度均值 $\overline \varTheta _{\rm{glit}}^Z$ 的线性关系表达式

$\left\{ \begin{aligned} & \varTheta _{\rm{glit}}^{lin} = \beta \cdot (\sum {\bar \varTheta _{\rm{glit}}^{{Z}}})/Z \pm \varDelta \\ & \overline \varTheta _{\rm{\rm{glit}}}^Z = (\sum {\sum {{\varTheta _{\rm{glit}}}(I, J, Z)} })/(I \times J) \end{aligned} \right.$

(21)

式中， $I, J, Z$ 分别表示单层图像像元行列位置和层数， $ \pm \varDelta $ 表示耀光角临界值增减量。 $\, \beta $ 为耀光角系数，其值与载荷所处位置相关，海拔越高值越大，地球曲率越大值越大。

由前面已知耀光角临界值满足 $20^\circ < \varTheta _{\rm{glit}}^{lin} < 40^\circ $ ，为便于数据计算和统计分析，将耀光角划为等距 $k$ (取整)间隔区，建立TOA耀光反射率 $R_{\rm{glit}}^{TOA}$ 和耀光区限定角度 ${\overline \varTheta _{\rm{glit}}}$ 反射率均值 $\overline R_{\rm{glit}}^{TOA}$ 的表达式

$R_{\rm{glit}}^{TOA}(I, J, Z;{\varTheta _{\rm{glit}}};k) = \alpha \cdot \overline R_{\rm{glit}}^{TOA}(I, J, Z;{\overline \varTheta _{\rm{glit}}};k) \pm \varDelta '$

(22)

式中， $\alpha $ 为反射率系数，其值同样与载荷位置相关， $ \pm \varDelta '$ 表示耀光反射率增减量。

在耀光区域范围内，以角度均值 ${\overline \varTheta _{\rm{glit}}}$ 为循环准则，计算在角度数值区间 $(c, c + k)$ 的NIR波段(865 nm) TOA反射率均值的斜率K

$K = \frac{{\overline R_{\rm{glit}}^{TOA}(I, J, Z;{{\overline \varTheta }_{\rm{glit}}}(c + k)) - \overline R_{\rm{glit}}^{TOA}(I, J, Z;{{\overline \varTheta }_{\rm{glit}}}(c))}}{k}$

(23)

研究发现耀光区域内辐射特性呈规则变化：从耀光中心(晴空区)开始随着 $c$ 的逐渐增加，反射率均值和偏振反射率均值对应逐渐递减，当循环到一定阶段，斜率K呈现出负正交替，反射率均值和偏振反射率均值渐渐趋于稳定。

    (3.3) OGDD算法流程

卫星在轨运行时，其目标辐射特性会随着时间和空间变化而发生一定改变，不同载荷间的差异特别明显，既有仪器的原因也可能是由算法的差异而引起。为统一实现单一或多个载荷不同轨道高度、不同经纬度及不同仪器的海洋耀光动态检测，本文结合偏振辐射特性规则表现，设计OGDD算法流程(如图3)。该流程分两个阶段，第1阶段获取待检测参考层，第2阶段动态获取耀光角临界值检测耀光。具体步骤如下：

步骤1　POLDER3一级产品预处理。首先读取POLDER3一级产品数据，再逐像元提取经纬度信息，并对照同经纬度的海陆标识码，从而获取海洋像元，最后形成研究区域的海洋一级衍生数据；

步骤2　读取研究海洋区域的一级产品衍生数据，获取每个像元的太阳及观测几何 $({\theta _{\rm{s}}}, {\varphi _s}, {\theta _v}, {\varphi _v})$ ；

步骤3　利用耀光角计算式(14)获得每个像元的耀光角 ${\varTheta _{\rm{glit}}}$ , 生成海洋像元耀光角数据文件并保存；

步骤4　利用辐射传输软件和辐射特性计算公式，结合外部数据(风速、风向、气溶胶光学厚度等)，模拟得到晴空条件下TOA海洋耀光区NIR波段的 $R_{\rm{glit}}^{TOA}$ 、 $PR_{\rm{glit}}^{TOA}$ 、 $cetR_{\rm{glit}}^{TOA}(I, J)$ 、 $cetPR_{\rm{glit}}^{TOA}(I, J)$ 、 $\overline R_{\rm{glit}}^{TOA}$ 、 $\overline {PR} _{\rm{glit}}^{TOA}$ ；

步骤5　用最邻近法进行数据配准，获取同样几何条件下，待检测海区耀光中心像元实测 $TR_{\rm{glit}}^{TOA}$ 、 $TPR_{\rm{glit}}^{TOA}$ 。若气溶胶光学厚度 $\tau $ 以步长 $s$ 逐渐递增，循环次数为m $(m = 0, 1, 2, \cdots, n)$ ，满足如下条件：

(1) $\varDelta _{I, J, Z}^R(\tau) > 0$ 且 $\varDelta _{I, J, Z}^R(\tau) > (\displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^R(\tau)})/\sum Z $ ，或 $\varDelta _{I, J, Z}^{PR}(\tau) > 0$ 且 $\varDelta _{I, J, Z}^{PR}(\tau) > \displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^{PR}(\tau)}/\sum Z $ ，判断中心像元 $\varDelta _{I, J, Z}^R$ 或 $\varDelta _{I, J, Z}^{PR}$ 离散程度较大，继续循环；

(2) $\varDelta _{I, J, Z}^R(\tau) > 0$ 且 $\displaystyle\frac{1}{2}(\displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^R(\tau)})/\sum Z < \varDelta _{I, J, Z}^R(\tau) < $ $\displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^R(\tau)}/\sum Z $ ，或 $\varDelta _{I, J, Z}^{PR}(\tau) > 0$ 且 $\displaystyle\frac{1}{2}(\displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^{PR}(\tau)})/ $ $\displaystyle\sum Z < \varDelta _{I, J, Z}^{PR}(\tau) < \displaystyle\sum\limits_Z {\varDelta _{I, J, Z}^{PR}(\tau)}/\sum Z $ ，判断 $\varDelta _{I, J, Z}^R$ 或 $\varDelta _{I, J, Z}^{PR}$ 离散程度不大，将该观测层数据作为待检测参考层，进入下一阶段处理；

步骤6　利用Buriez云检测算法(Buriez 等，1997)剔除待检测层海区的有云像元；

步骤7　选定待检测参考层，计算角度区间 $(c, c + k)$ 的晴空TOA反射率均值 $\overline R_{\rm{glit}}^{TOA}(I, J, Z;$ ${\overline \varTheta _{\rm{glit}}}(c + k))$ ；

步骤8　TOA反射率均值的斜率K变化检测。计算角度区间 $(c, c + k)$ 与 $(c + k, c + 2k)$ 反射率均值函数的斜率 $K(c)$ 和 $K(c + k)$ 。若满足斜率 $K(c) < 0$ 且 $K(c + k) > 0$ ，标记首次斜率规则变化点，记为耀光角临界均值拐点 $\overline \varTheta _{\rm{glit}}^Z(c + k)$ ；

步骤9　经过系列优化筛选后，得到Z层的耀光角临界均值拐点集合 $\{ (Z, \overline \varTheta _{\rm{glit}}^Z(c + k))\} $ ，根据式(21)可以得到耀光角临界值 $\varTheta _{\rm{glit}}^{lin}$ ；

步骤10　对所有像元进行耀光角规则判别，满足 ${\varTheta _{\rm{glit}}} < \varTheta _{\rm{glit}}^{lin}$ 条件的像元都标记为耀光像元 $pix(I, J, Z;{\varTheta _{\rm{glit}}})$ 。

4、实验数据及结果分析     (4.1) 数据源选择及耀光角计算

法国于2004年12月成功发射PARASOL/POLDER3卫星载荷后，其观测波段由3个偏振和6个非偏振波段组成(表1), 沿太阳同步轨道由南向北进行成像观测(图4)，对同一目标进行连续多次重复观测(多角度观测)。由法国国家空间研究中心(CNES)面向全球用户提供3级业务化产品，现其数据产品可以从法国里尔第一大学的ICARE数据服务中心(www.icare.univ-lille1.fr[2017-03-16])下载。本文针对高分五号卫星DPC载荷数据预处理中如何动态检测海洋耀光进行研究，数据源选择POLDER3一级产品作为实验对象，空间分辨率为6×6 km, 投影方式为Sinusoidal投影(Riédi 等，2010)。筛选日期为2005-03-04的单条带历史数据，文件名为P3L1TBG1006137SD和P3L1TBG1006137SL，该实验数据位于印度洋区域，耀光特征明显且天气相对晴朗少云，能有效减少云层对耀光的干扰。该条带实验数据总共进行了129次多角度成像观测，依据波段生成存储数据组，每个波段组由3维结构组成(行×列×观测层)。用Anapol软件(POLDER和PARASOL数据可视化及分析的IDL图形界面工具，可从ICARE网站下载)打开该条带数据并预览(如图5(a)所示)，发现该条带数据海洋区域靠近赤道以南部分有大量的云层覆盖，而赤道以北靠近非洲东海岸附近的印度洋海区相对晴朗少云，因此截取该海区一个相对晴朗的实验区域作为研究对象。该区域共有29个观测层(55—83层)，坐标经度为49.6°E—66.0°E，纬度为6.6°S—7.0°N，如图5(b)所示。

从图5(b)截取的实验海区可知海洋耀光随着卫星轨道方向移动，不同观测角度下耀光位置不同，利用OGDD模型中的计算公式对所有观测层进行有云耀光角计算和实测偏振辐射信息提取，可得到如图4所示研究区域NIR波段(865 nm)的反射率值(以69层为例)。

图6中，蓝色线表示等角线(耀光角)，最内环线表示耀光角值为10°，线与线间隔10°，反射率值以百分数形式表示。从图中可以看到耀光角值在20°范围内确定全部为耀光，在20°—30°耀光对比并不是很明显，且耀光与薄云间的辐射信息会有重叠。分析实验海区NIR波段的伪彩色图，发现耀光区耀光角值由内向外逐渐增加，可看成由等间距值的环线(等角线)构成，且等角线递增的角度值与反射率递减趋势一致，通过耀光区像元的偏振辐射特性与耀光角间的规则变化这一规律可以构建它们之间的线性关系表达(式(21))。剔除实验海区受云层干扰较大的区域，减少有云海区的范围，将优化后的数据利用OGDD模型进一步处理。

    (4.2) 耀光偏振辐射特性计算及分析

从图6实验海区中优化选取赤道以北云污染较小的相对晴空区(经度为52.7°E—57.8°E, 纬度为1.2°N—6.2°N)，得到如图7所示(第73观测层)少量云污染的实验区域。对该区域所有观测层数据，用耀光角及散射角计算公式，获取每层数据所有像元的耀光角值和散射角值，并以数据表形式存储。图7(a)中的蓝色线为耀光角等角线，红色线为散射角等值线。将实验区域第73层所有像元的耀光角值生成如图7(b)所示的角度值栅格图，纵轴和横轴对应像元所处行列号。图7(b)中蓝色为角度值较小的耀光区，蓝色和黄色分界线在20°左右，大于40°耀光角的像元比较少。若采用MODIS经验临界值来判(40°)别耀光像元，则该实验区几乎全部被标记为耀光像元，显然降低了有效像元的比率。

由前可知，在获取耀光角值后，计算耀光区偏振辐射特性模拟值，并与实测值进行数据配准，比较分析其离散程度确定待检测参考层。逐像元获取对应经纬度坐标的海表风速风向 $v$ ，气溶胶光学厚度 $\tau $ 等参量(Kotchenova和Vermote，2007)，通过大气矢量辐射传输软件6SV1.1(http://6s.ltdri.org[2017-03-16]下载)模拟计算每个像元(865 nm)的TOA反射率 $R_{\rm{glit}}^{TOA}$ 和偏振反射率 $PR_{\rm{glit}}^{TOA}$ 。海表风速风向值通过ECMWF获取，由于ECMWF数据和POLDER3一级产品空间分辨率不同，为准确获取对应位置和时间的风速风向v，必须进行数据配准，采用最近邻域法提取对应经纬度像元的风速风向。气溶胶光学厚度 $\tau $ 值范围为[0.0001, 2.0]，步长0.1。设置观测和太阳几何，生成该区域的模拟TOA反射率及偏振反射率查找表LUT。通过LUT中对应辐射参量可以获得反射率及偏振反射率模拟值。

根据TOA耀光中心偏振辐射特性计算方法，获得如图8(a)所示实验区域第71—75层耀光中心点(观测天顶角依次为：17.2°，17.5°，18.1°，18.5°，19.3°；太阳天顶角依次为：17.1°，17.5°，18.0°，18.5°，19.0°)在风速为5 m/s，气溶胶光学厚度 $\tau $ 值变化(从LUT获取0.0001、0.1001、0.2001、0.3001)的4组反射率及偏振反射率模拟值 $(cetR_{\rm{glit}}^{TOA}, cetPR_{\rm{glit}}^{TOA})$ 与实测值 $(TR_{\rm{glit}}^{TOA}, TPR_{\rm{glit}}^{TOA})$ 。

图8(a)中五角星线和五边形线分别表示POLDER3实测的 $TR_{\rm{glit}}^{TOA}$ 、 $TPR_{\rm{glit}}^{TOA}$ 。当海表风速 $v$ 为5 m/s，气溶胶光学厚度 $\tau $ 值分别取0.0001、0.1001、0.2001、0.3001时，随观测天顶角的增加不同观测层的模拟 $cetR_{\rm{glit}}^{TOA}$ 和模拟 $cetPR_{\rm{glit}}^{TOA}$ 有缓慢递增趋势。观测天顶角不变，而气溶胶光学厚度逐渐增加时，模拟 $cetR_{\rm{glit}}^{TOA}$ 和模拟 $cetPR_{\rm{glit}}^{TOA}$ 值有下降的趋势，与实测值离散程度越来越大。因在完全晴空条件下，大气分子散射和吸收相对总辐射贡献较小，随着 $\tau $ 的增加( $\tau < 2$ )，气溶胶对太阳辐射散射和吸收，使TOA反射率值和偏振反射率值变小。图8(b)左边为 $\tau $ 取不同值时，不同观测天顶角下的耀光中心实测 $TPR_{\rm{glit}}^{TOA}$ 与模拟 $cetPR_{\rm{glit}}^{TOA}$ 的绝对偏差 $\Delta _{{cet}}^{PR}(\tau)$ ，右边为实测 $TR_{\rm{glit}}^{TOA}$ 与模拟 $cetR_{\rm{glit}}^{TOA}$ 的绝对偏差 $\Delta _{{cet}}^R(\tau)$ 。由TOA反射率实测值与模拟值相对偏差值( $\Delta _{{cet}}^R(\tau)$ , $\Delta _{{cet}}^{PR}(\tau)$ )，发现第72、73层值绝对偏差较大(最大值达6.7)，表明离散程度最大，因为这两层耀光中心反射辐射与薄云反射辐射叠加，实测值会有较大增加，第74层受海风误差(实际海风大于5 m/s)的影响耀光中心实测值会比模拟值低(绝对偏差最小值达–2.3)，当 $\tau $ 值取0.0001时第71、75层由于未受云和海风干扰模拟值和实测值比较接近，其离散程度较小，故可将其选择作为耀光区偏振辐射特性计算的参考层。

    (4.3) OGDD算法结果分析

考虑耀光区通常有大量的云层覆盖，云层对TOA反射率和偏振反射率值有较大影响。由前面分析耀光中心偏振辐射特性方法获取耀光区的参考层，再利用OGDD算法来判别耀光角临界值。由于云层和大气干扰了辐射信息，选择受大气干扰较少的NIR波段(865 nm)来提取辐射特性均值。先对分析区域进行去云处理(Buriez 等，1997)。为了使统计更具有普遍性，利用以上原则选取如图3(b)所示的位于印度洋西部靠近沙特的实验海区对耀光区进行OGDD算法结果分析。该选取的实验海区数据共有29个观测层，编号从55层始至83层。利用耀光区偏振辐射特性计算方法对每层的数据进行计算和结果分析，可以获得不同间隔的数据统计情况。由前述可知将第71、75层选作参考层，本实验选定第71层数据作为参考层进行计算(取 $\beta = 1$ )，表2为第71层晴空条件下耀光临界值约束下耀光区各个参量(包括区域内总耀光像元数Tpix(N)、间隔区含云耀光像元数pix_c(N1)、间隔区无云像元数pix(N2)、间隔区无云像元数比总像元数pix(N2)/Tpix(N)、无云耀光像元数比40°区域内总像元数pix_c(N1)/pix_40(N3)、区域内平均TOA反射率 $\overline R_{\rm{glit}}^{TOA}(i, j, z;{\overline \varTheta _{\rm{glit}}};k)$ )随角度变化的计算结果统计情况。

对实验海区29层数据进行初步筛选，统计第62—78层所有参量数据，获取每个层在不同耀光角临界值约束条件下的无云像元的TOA反射率均值(图9(a))。在耀光角20°—21°范围内第62—78层(73层除外)由于受云的干扰较大，该临界值区域内的无云像元TOA反射率均值为0。第62、63、77、78层由于处于卫星观测视场进入和离开的观测层，耀光角选择区域未进入观测视场导致无云像元TOA反射率均值为0的情况比较多。由于受观测视场范围和云干扰的影响，使得耀光角临界值的获取较为困难。

据上可知，考虑载荷在对地观测时，其成像质量受视场范围限定和云干扰，需对观测层及TOA反射率数据进一步优化，依据TOA反射率值有效(非0)原则删除第62、71、72、74—78层TOA反射率均值，获取图9(b)所示筛选优化后的第63—73层(不含71, 72层)30°—40°范围的TOA反射率均值。从第71、75参考层中计算无云像元(限定耀光角大于20°) TOA反射率均值，取其最大值2.9%，将其作为所有观测层反射率均值的参考值(图9(b)水平方向红色虚线)。对每层小于参考值的TOA反射率均值保留，再对每层保留的反射率均值进行变化趋势分析。计算耀光角均值 $\overline \varTheta _{\rm{glit}}^{\textit{z}}$ 在 $(c, c + k)$ 范围对应的TOA反射率均值规律变化的斜率K，以K首次正负交替变化提取每层的耀光角临界值拐点，生成集合 $\{ ({\textit{z}}, \overline \varTheta _{\rm{glit}}^{\textit{z}}(c + k))\} $ (如图9(b)中所示凸出显示的点标识)。结果集为{(63, 33)，(64, 35)，(65, 35)，(66, 33)，(67, 35)，(68, 33)，(69, 34)，(70, 35)，(73, 35)}，对结果集内参考层所有耀光角临界值拐点求均值，得到实验海区的耀光角临界值34°(图9(b)垂直方向蓝线虚线)。最后将所有观测层的耀光角临界值 $\varTheta _{\rm{glit}}^{lin}$ 都动态调整为34°，根据调整后的 $\varTheta _{\rm{glit}}^{lin}$ 对实验海区进行耀光像元标记。

为客观评价OGDD方法效果，选取图5(b)中的实验海区，分别用OGDD方法与POLDER3方法及MODIS方法对实验海区进行58—81层耀光检测。以耀光角临界值 $\varTheta _{\rm{glit}}^{lin}$ 分别取20°(起始临界值)、30° (POLDER3业务化算法经验值)、34°(本文OGDD算法获取临界值)、40°(Modis业务化算法经验值)以及不限定耀光角(Ulimit) 5种耀光角临界值情形进行统计，获取剔除耀光像元后晴空像元TOA反射率均值{ $\overline R_{\rm{clear}}^{TOA}({\varTheta _{\rm{glit}}} > \varTheta _{\rm{glit}}^{20^\circ })$ , $\overline R_{\rm{clear}}^{TOA}({\varTheta _{\rm{glit}}} > \varTheta _{\rm{glit}}^{{\rm{POLDER}}})$ , $\overline R_{\rm{clear}}^{TOA}({\varTheta _{\rm{glit}}} > \varTheta _{\rm{glit}}^{{\rm{OGDD}}})$ , $\overline R_{\rm{clear}}^{TOA}({\varTheta _{\rm{glit}}} > \varTheta _{\rm{glit}}^{{\rm{MODIS}}})$ , $\overline R_{\rm{clear}}^{TOA}(\varTheta _{\rm{glit}}^{U{\rm{limit}}})$ }。结合本文提出的OGDD算法在耀光角[20°, 40°]范围内，以TOA反射率均值的斜率K规则变化动态调整耀光角临界值为34°，剔除实验海区临界值内耀光像元和所有云污染像元后，获取剩余晴空像元的 $\overline R_{\rm{clear}}^{TOA}({\rm{OGDD}})$ 。统计结果表明OGDD方法相比现有的MODIS和POLDER耀光检测方法有明显算法优势(图10)。

(1) 不限定耀光角及耀光角满足 ${\varTheta _{\rm{glit}}} > 20^\circ $ 时(绿线和紫红线)，印度洋区域晴空无云像元受耀光像元的影响程度接近，58—81层的 $\overline R_{\rm{clear}}^{TOA}(\varTheta _{\rm{glit}}^{U{\rm{limit}}})$ 与 $\overline R_{\rm{clear}}^{TOA}({\varTheta _{\rm{glit}}} > 20^\circ)$ 值非常接近且大部分重合，表明耀光角20°外未被标记为耀光的像元和不剔除耀光像元对晴空反射率的影响几乎相同。因此 $\varTheta _{\rm{glit}}^{lin}{\rm{ = }}20^\circ $ 得到的晴空像元 $\overline R_{\rm{clear}}^{TOA}$ 受耀光的干扰仍然较大，会干扰晴空区域气溶胶反演，无法提供合格数据。

(2) 利用POLDER3经验值满足 ${\varTheta _{\rm{glit}}} > 30^\circ $ (蓝线)，剔除实验海区的耀光和有云像元后，30°区域外晴空像元 $\overline R_{\rm{clear}}^{TOA}$ 与 ${\varTheta _{\rm{glit}}} > 40^\circ $ 时近似完全晴空条件下 $\overline R_{\rm{clear}}^{TOA}$ 有一定差异。以30°耀光角作为耀光检测的临界值并不十分准确，得到的晴空像元仍受到部分耀光的影响，使得标记为晴空像元多于实际晴空像元。

(3) 利用MODIS经验值满足 ${\varTheta _{\rm{glit}}} > 40^\circ $ (黑线)，剔除实验海区的耀光和有云像元后，晴空无云像元的NIR波段 $\overline R_{\rm{clear}}^{TOA}$ 随观测层的变化趋于稳定。大于40°得到晴空像元的 $\overline R_{\rm{clear}}^{TOA}$ 为理想晴空条件下均值，许多晴空像元标记为耀光被剔除，有效像元利用率较低。

(4) 通过本文OGDD算法获取 $\varTheta _{\rm{glit}}^{lin}{\rm{ = 34}}^\circ $ 来检测耀光(红线)，显然当 ${\varTheta _{\rm{glit}}} > 34^\circ $ 时区域外晴空像元 $\overline R_{\rm{clear}}^{TOA}$ 非常接近 ${\varTheta _{\rm{glit}}} > 40^\circ $ 时几乎完全晴空下的 $\overline R_{\rm{clear}}^{TOA}$ (黑线)。与MODIS经验值耀光检测进行对比，耀光临界值相对减少15%，标记为耀光的像元减少了30%，增加了像元的利用率。相比POLDER3方法，临界值相对增加13%，标记为耀光的像元增加了17%，降低了耀光像元对晴空区气溶胶反演的干扰。

因此使用OGDD方法可以动态获取印度洋海区的最优化耀光角临界值，通过该值能比较准确的检测出海洋耀光，从而大大的提高了有效像元的比率，为后续云和气溶胶特性反演提供更加准确的数据产品。

5、结　论

海洋耀光对海洋上空云和气溶胶研究干扰较大，准确识别耀光对卫星数据的高效处理及高质量生产至关重要。现有耀光检测方法大都利用经验值进行判别且不考虑波段选择。利用这些经验值检测不同载荷的耀光，往往会降低有效像元的利用率，导致云和气溶胶反演误差增大，故不适于DPC的海洋耀光检测。因此，动态获取耀光角判别临界值是实现海洋耀光的动态检测的关键。

本文针对海洋耀光动态检测这一难题，提出了NIR偏振辐射数据的OGDD模型。现有的耀光检测基本不考虑波段的选择，而作者通过水面耀光实验发现耀光和非耀光区NIR反射率差异显著，表明利用NIR波段检测耀光更准确。该模型的第1阶段结合菲涅尔反射和粗糙海表偏振辐射理论，构建了耀光区TOA偏振辐射特性的计算公式。利用其获取相应参量的数值，对比耀光中心参量实测值与模拟值的离散程度，进而确定待检测参考层。在第2阶段剔除有云像元后，建立TOA反射率均值斜率规律变化准则，并利用其进行斜率动态分析，生成每层耀光角临界值集合并求均值，得到待检测耀光区的耀光角临界值，通过动态获取该临界值最终实现海洋耀光的动态检测。

利用OGDD模型可动态获取实验海区的最优化耀光角临界值，并对比分析5种不同耀光角临界值条件下获取的晴空TOA反射率均值。研究发现该模型能有效的区分耀光像元和非耀光像元，大大提高像元的利用率，降低耀光对晴空区气溶胶反演的干扰。此外，该模型还克服了现有耀光检测方式单一、适用性差的问题，可用于不同载荷耀光临界值的动态调整，实现跨平台检测，为高分5号卫星DPC载荷的海洋耀光处理提供理论和技术支持。

志 谢 本文的研究数据和工具由法国ICARE数据服务中心提供，在此表示衷心的感谢！