赤潮水体水质遥感参数时空变化分析

陈曦，汪小钦，陈芸芝

(1.福州大学空间数据挖掘和信息共享教育部重点实验室，福州 350002；2.福建省空间信息工程研究中心，福州 350002)

1 引 言

赤潮(Harmful Algal Blooms，简称为HABs)是海水中某些浮游生物在一定的环境条件下，突发性地增殖或集聚引起的生态现象，被列为世界三大近海污染问题之一。卫星遥感具有覆盖范围广、快速、同步的优势，已经被大量用于赤潮的探测监测及预警预报研究和应用。

EOS卫星上搭载的中分辨率成像光谱仪MODIS具有较短的重访周期、较高的空间分辨率以及合理的波段设置，被广泛应用在赤潮监测中。王其茂等[1]、丛丕福等[2]通过分析赤潮水体的光谱特性，建立了多波段比值模型用以提取赤潮信息；李继龙等[3]利用MODIS影像真彩合成、多波段差值比值、海表温度及叶绿素a浓度阈值对2004年、2005年长江口及邻近海域的赤潮进行了探测判定；孙丽雅等[4]分别采用叶绿素浓度固定阈值法、月平均值差值法和海温及多光谱综合法对2008年5月东海水体的异常信息进行提取分析，并比较分析了各种方法的提取结果；Yongmin Kim等[5]对2003年9月30日和10月2日韩国周边海岸的MODIS进行了信息提取，基于赤潮指数MRI提出了两步分离的方法，即海温和667nm波段的波长逐步分离，最后用Moran’s I空间统计的方法探测了赤潮，并用误差矩阵、F检验、Jaccard和Yule系数进行了信息提取的精度验证；Carvatho等[6]利用2002年～2006年的MODIS数据研究了佛罗里达州的赤潮，并开发出了一种混合算法来判定赤潮，通过混合算法与单一的传统经验算法和生物光学方法比较，发现利用混合算法判定的赤潮精度高于单一的传统经验算法和生物光学方法；Son等[7]利用2011年的MODIS数据开发了一种新的赤潮判定方法，这种新方法综合使用了荧光高度法、双波段比值法、叶绿素a反演法和生物光学方法，通过4步判定实现对赤潮的判定，并以韩国南部海域的5次赤潮事件加以验证；Masoud Moradi和Keivan Kabiri[8]研究了2008年11月波斯湾东部发生的赤潮，发现用MODIS荧光高度法提取的赤潮信息结果优于叶绿素浓度异常法。

本文以2011年发生在厦门附近海域的赤潮为例，综合使用叶绿素浓度异常法、海温异常法和多波段比值法，侧重分析赤潮发生前后MODIS数据相关水质参数在时间、空间上的变化情况，辅以比较赤潮信息提取方法，以期为今后赤潮的预报分析提供支持。

2 研究区域与数据

赤潮是浮游生物在短期内迅速繁殖并在海面积聚的结果[9]，其形成条件是有赤潮生物存在、营养物质供应丰富、水温适宜、水体动力环境平稳，所以与赤潮相伴随的往往是叶绿素浓度迅速上升、海表温度上升和表层海水光谱性质的变化。

本文的研究区域是邻近福建的东海海域，是赤潮易发区域。由中国海洋信息网[10]数据获知：2011年7月26日～2011年8月7日福建省厦门同安湾顶琼头海域以及鳄鱼屿以南、集美大桥至五缘湾大桥一带海域发生赤潮，最大面积为105km2，引发赤潮的是中肋骨条藻(Skeletonema costatum)。

根据这些统计数据，从NASA(美国国家航空与航天局)免费下载了MODIS 1km数据，从预览图中选取无云层遮挡、覆盖研究区域的赤潮发生前后的6景数据，包括赤潮发生前(7月23日)、赤潮发生期间(7月26日两景，7月31日)、赤潮发生后(8月8日，8月11日)。

3 数据处理

首先对下载的MODIS L1B数据进行预处理，包括辐射校正、几何校正，利用附带的辐射定标数据进行定标，得到MODIS L2级数据，整个过程使用NASA发布的SeaDAS软件完成。对L2级数据截取24°N～25°N，118°E～120°E区域作为研究区域。

利用生成的L2数据，提取研究区域水质参数进行分析，特征参数包括：海水叶绿素a浓度(Chl-a)、海水表面温度(SST)、海水反射率(Rrs)、归一化离水辐射率(nLw)。采用的技术路线如图1所示。

图1 技术路线

3.1 叶绿素a浓度反演

通过遥感影像反演叶绿素a浓度，就必须建立叶绿素a浓度与海水反射率或者离水辐射率之间的函数关系。根据张春桂[11]等人的研究，通过对比MODIS OC2标准经验算法、OC3标准经验算法和Clark半分析算法，发现OC3标准经验算法比较适合于福建近海叶绿素a浓度反演，所以本文采用MODIS OC3标准经验算法[11]：

Cchla=10(0.283-2.753R+1.457R2+0.659R3-1.403R4)

(1)

R=lg[max(Rrs(443,488))/Rrs(551)]

(2)

其中，Cchla表示叶绿素a浓度，R表示OC3算法的波段比率，Rrs(λ)为波段λ的反射率，其中，443，488，551分别对应MODIS第9、10、12波段的中心波长。提取的叶绿素a浓度时间序列如图2所示。

3.2 海水表面温度反演

海水表面温度(SST)是重要的海洋环境参数，可以辅助判别赤潮的发生与否。通用的海水表面温度反演算法包括[12]：

SST=c1+c2*T31+c3*|T31-T32|+c4*(secθ-1)(T31-T32)

(3)

(4)

公式(3)中，c1-4为回归系数，θ为卫星天顶角，可以通过辐射定标数据获得，Ti为波段i的亮温，通过公式(4)得到；公式(4)中C1=3.74151×10-22(W·m2)，C2=0.0143879(m·K)，λi是i波段的中心波长(m)，I(λi)是热红外波段的光谱辐照度。

公式(3)中，c1～4通过回归得到，其数值为4.9261、0.9888、0.1301、1.627。提取的海水表面温度时间序列如图3所示。

图3 时间序列海水表面温度

3.3 Rrs、nLw比值模型

大量研究表明，正常海洋水体和赤潮水体的光谱存在明显差异[13-16]：赤潮水体在波长440nm～460nm和650nm～670nm分别有两个吸收峰，在波长560nm～580nm和690nm～710nm有两个反射峰，而正常海洋水体没有。

MODIS第3波段(459nm～479nm)信息中含有水体的吸收信息，第4波段(545nm～565nm)信息中含有赤潮水体的反射信息，可以利用这两个波段反射率的比值建立模型。再辅助使用MODIS第9(438nm～448nm)、第11(526nm～536nm)波段的归一化离水辐射率比值来识别赤潮。

Rrs(551)/Rrs(469)>Cr

(5)

nLw531/nLw443>Cnw

(6)

其中，Rrs(469)、Rrs(551)是MODIS第3、4波段的反射率；nLw531、nLw443是MODIS第9、第11波段归一化离水辐射率，Cr、Cnw是赤潮水体识别常数，其值的大小与发生赤潮的海区和赤潮藻种有关。由于波段比值区间范围差异大，无法使用统一的图例，所以本文仅展示一景为例，见图4、图5。

图4 Rrs551/Rrs469

图5 nLw531/nLw443

4 结果分析

4.1 水质参数时间变化规律分析

从图2可以看出7月23日厦门大金门岛以南海域叶绿素a浓度(Chl-a)已经超出正常水平，有少量海域已经达到了10mg/m3以上，无论Terra还是Aqua影像都显示7月26日厦门岛以东海域Chl-a很高，而到7月31日Chl-a迅速降低，低于10mg/m3以下，之后的影像显示Chl-a正常。

从图3可以看到从7月23日到8月11日，沿海岸海水的表面温度都高于远离海岸的地区，相差将近10℃，而7月26日的Terra影像显示远离海岸地区的海水表面温度整体高于其他时相，7月26日的Aqua影像上厦门大金门岛以西、以南海域海水表面温度骤降，8月8日、8月11日的影像海水表面温度较低，保持在25℃左右。

从图4、图5中看出近海水域的比值大于远洋区域，图中发现的零星亮点是比值较大的区域，可以判定为赤潮可能发生的区域。

叶绿素a浓度是最主要的赤潮特征参数，海水表面温度是海水中藻类繁殖的重要客观条件，对叶绿素a浓度、海水表面温度进行统计。首先对特征提取后的图像作掩膜处理，去掉陆地、云及耀斑等无效数据，只留下有效的海洋水体信息，设定0～1000mg/m3为有效叶绿素a浓度值，-2℃～45℃为有效温度，然后依次计算统计叶绿素a浓度、海水表面温度，得到每一时相的叶绿素a浓度和海水表面温度的最小值、最大值、平均值、方差(表1、表2)。

从表1中可以看出：6个时相的叶绿素a浓度平均值、方差变化不大，7月26日Aqua和Terra影像的最大值远远大于其他时相。从表2中可以看出：6个海水表面温度平均值基本上稳定在26℃～28℃，7月26日Aqua影像最大值略大于其他时相，7月26日Terra影像的方差大于其他时相。

表1 MODIS叶绿素a浓度统计数据

表2 MODIS海水表面温度统计数据

综合以上，可以得到叶绿素a浓度在赤潮发生前就开始在厦门大金门岛以南缓慢增加，赤潮发生时厦门岛以东大面积快速增高，赤潮退散后降低；海水表面温度除了赤潮最为严重时先增后减外其余变化平稳；Rrs、nLw比值不能明确显示赤潮发生，仅在赤潮最为严重的时候有零星高值。

4.2 赤潮水体与非赤潮水体特征分析

从上面的图表不能体现赤潮水体和非赤潮水体在赤潮发生前、中、后的变化特征，因此根据中国海洋信息网资料数据[10]在厦门大金门岛屿附近以及泉州附近地区分别截取两个大小基本相同的AOI以代表赤潮水体和非赤潮水体。然后统计各自范围内有效像元的叶绿素a浓度、海水表面温度、Rrs、nLw比值。叶绿素a浓度、海水表面温度有效值同之前设置的一致，在Rrs、nLw有效值设置上，为了剔除两个波段中陆地、云层等非水体区域，即比值小于等于1的区域，设1.000001～100000为有效值。赤潮水体和非赤潮水体的特征参数时间变化如图6所示，其中，7月26日数值是Terra和Aqua的平均值。

图6 赤潮与非赤潮水体特征参数时间变化趋势

从图6可以看出，赤潮水体的各项数值均大于非赤潮水体；除海水表面温度以外，其余三种特征的曲线很相似，在7月26日有一个明显的峰值，7月31日以后数值平稳回落。

非赤潮水体水质参数在赤潮发生前、中、后变化不大，拐点不明显；赤潮水体水质参数赤潮发生前、中、后变化明显，且在赤潮发生时有明显的数值拐点。

图6(b)海水表面温度的曲线与其他特征曲线有所差异，通过进一步分析表明：海水表面温度不是赤潮发生的决定性因素，根据温度变化无法判断是否发生赤潮，上升流也可能带来温度的变化，而海表温度法对表面温度变化不大的赤潮也不适用[17]。此外，根据霍文毅等[18]的研究得知，中肋骨条藻在海水温度0℃～37℃内均可生长，其最适宜增殖温度为24℃～28℃，作为对比，东海原甲藻(Prorocentrum donghaiense)的适温范围是16℃～23℃[19]，其爆发时的水温为17℃～20℃[20-21]，因此利用海水表面温度可以辅助判别赤潮藻类类型。

4.3 赤潮分布区域提取与分析

通过以上对赤潮和非赤潮水体特征的变化分析，可以综合利用Chl-a、Rrs比值、nLw比值等参数进行赤潮信息的提取。要准确地提取出赤潮信息，合理设置各参数特征阈值是关键。

本文中叶绿素a浓度阈值的选取使用固定阈值法，通过叶绿素a浓度平均值与方差的宏观比较以及赤潮水体和非赤潮水体的叶绿素a浓度对比分析，将叶绿素a浓度阈值设为赤潮水体时间7月26日的均值10mg/m3。

Rrs比值、nLw比值的阈值设定与叶绿素a浓度阈值选取相同，直接把赤潮水体时间7月26日的峰值设为阈值，即Rrs551/Rrs469阈值1.8，nLw531/nLw443阈值2.8。

将叶绿素a浓度、Rrs比值、nLw比值异常图叠加分析，得到异常图的时间序列，如图7所示，其中红色代表赤潮发生区域，蓝色为正常区域。

从图7可以看出：赤潮从7月23日起就出现在厦门大金门岛以南的小范围海域，集中爆发是在7月26日，聚集在厦门岛以东区域，7月31日赤潮基本退去。

图7 水体异常图(图中蓝色区域表示海域，红色区域表示赤潮)

根据新华网对此次赤潮的报道[22]，7月26日，赤潮面积约25km2，主要位于厦门同安湾顶琼头以及鳄鱼屿以南海域；7月28日，赤潮面积扩大到约105km2，其中厦门同安湾海域扩大到约60km2，厦门西海域出现约25km2，东部海域出现约20km2；7月30日，赤潮面积减少到20km2，其中厦门同安湾15km2，西海域减少到5km2，东部海域赤潮消失。本文所提取的不同时间段的赤潮分布与该报道大体吻合。

5 结束语

本文通过分析2011年福建厦门海域赤潮水体与非赤潮水体的特征变化，得到如下主要结论：

①叶绿素a浓度在赤潮发生前变化很小，赤潮发生时迅速升高，赤潮退去后平稳下降；海水表面温度变化平稳；Rrs比值、nLw比值在赤潮最为严重时出现极大值。赤潮水体叶绿素a浓度、海水表面温度、Rrs比值、nLw比值均大于非赤潮水体，赤潮水体在赤潮发生时变化明显，有明显极值点。

②利用时间序列叶绿素a浓度变化提取赤潮的效果最好，海水表面温度可以用于辅助判别赤潮藻类，Rrs比值、nLw比值可以辅助判别出赤潮最为严重的区域。

由于云层遮挡以及影像质量等原因，使得遥感数据的时间序列不是很完整，且MODIS 1km数据在空间分辨率上不具有优势，对更为准确地提取赤潮面积信息有较大的限制，下一步将开展更多赤潮案例的研究，完善基于时间序列遥感数据的赤潮监测方法。

参考文献：

[1] 王其茂，马超飞，唐军武，等.EOS/MODIS遥感资料探测海洋赤潮信息方法[J].遥感技术与应用，2006，21(1)：6-10.

[2] 丛丕福，曲丽梅，王臣立，等.基于MODIS 模拟的辽东湾叶绿素a的遥感反演模型[J].生态环境学报，2009，18(6)：2057-2060.

[3] 李继龙，唐援军，郑嘉淦，等.利用MODIS遥感数据探测长江口及邻近海域赤潮初步研究[J].海洋渔业，2007，29(1)：25-30.

[4] 孙丽雅，陈芸芝，汪小钦.时间序列MODIS数据赤潮信息提取研究[J].遥感信息，2012，27(3)：71-77.

[5] KIM Y，BYUN Y，KIM Y，et al.Detection of Cochlodinium polykrikoides red tide based on two-stage filtering using MODIS data[J].Desalination，2009，249(3)：1171-1179.

[6] CARVALHO G A，MINNETT P J，FLEMING L E，et al.Satellite remote sensing of harmful algal blooms：A new multi-algorithm method for detecting the Florida Red Tide(Karenia brevis)[J].Harmful Algae，2010，9(5)：440-448.

[7] SON Y B，ISHIZAKA J，JEONG J C，et al.Cochlodinium polykrikoides red tide detection in the South Sea of Korea using spectral classification of MODIS data[J].Ocean Science Journal，2011，46(4)：239-263.

[8] MORADI M，KABIRI K.Red tide detection in the Strait of Hormuz(east of the Persian Gulf)using MODIS fluorescence data[J].International Journal of Remote Sensing，2012，33(4)：1015-1028.

[9] 顾德宇，许德伟，陈海颖.赤潮遥感进展与算法研究[J].遥感技术与应用，2003，18(6)：434-440.

[10] 中国海洋信息网.2011年中国海洋环境质量公报[EB/OL].http：//www.coi.gov.cn/gongbao/nrhuanjing/nr2011/201207/t20120710_23200.html.

[11] 张春桂，曾银东，张星，等.海洋叶绿素a浓度反演及其在赤潮监测中的应用[J].应用气象学报，2007，18(6)：821-831.

[12] 刘良明，周军元.MODIS数据的海洋表面温度反演[J].地理空间信息，2006，4(2)：7-9.

[13] 丘仲锋，崔廷伟，何宜军.基于水体光谱特性的赤潮分布信息MODIS遥感提取[J].光谱学与光谱分析，2011，31(8)：2233-2237.

[14] 王林，赵冬至，杨建洪，等.赤潮对近岸水体生物光学特性的影响[J].环境科学，2011，32(10)：2855-2860.

[15] 申力，许惠平，吴萍.基于光谱特征的MODIS影像东海水体类别信息提取[J].遥感信息，2011(3)：71-76.

[16] 张涛，苏奋振，杨晓梅，等.MODIS遥感数据提取赤潮信息方法与应用——以珠江口为例[J].地理信息科学学报，2009，11(2)：244-249.

[17] 吴文谨.基于光谱曲线特性和波谱角分类的赤潮监测方法[J].遥感信息，2009(8)：50-54.

[18] 霍文毅，俞志明，邹景忠，等.胶州湾中肋骨条藻赤潮与环境因子的关系[J].海洋与湖沼，2001，32(3)：311-318.

[19] 夏平，陆斗定，朱德弟，等.浙江近岸海域赤潮发生的趋势与特点[J].海洋学研究，2007，25(2)：47-56.

[20] 陈翰林，吕颂辉，张传松，等.2004年东海原甲藻赤潮爆发的现场调查和分析[J].生态科学，2006，25(3)：226-230.

[21] 张秀芳，刘永健.东海原甲藻Prorocentrum donghaiense Lu 生物学研究进展[J].生态环境，2007，16(3)：1053-1057.

[22] 新华网.厦门海域发生赤潮暂未发现有毒赤潮藻类[EB/OL].http：//news.xinhuanet.com/2011-07/31/c_121749922.htm.