基于高光谱的射阳河口悬浮泥沙浓度定量反演研究

潘 洁,张 鹰,谭子辉

(1.南京林业大学 森林资源与环境学院,江苏 南京 210037;2.南京师范大学 地理科学学院,江苏 南京210042;3.山东省临沂市蔬菜办公室,山东 临沂276000)

水体悬浮物含量是重要的水质参数,国内外诸多研究建立了多种基于遥感技术反演水体悬浮物含量的理论和经验模型,常用的主要有对数模式、Gordon模式和指数模式[1]。此外,我国学者还提出了一些改进算法,李炎等[2]提出以NOAAAVHRRCH1和CH2波段反射率差的最大值条带作为大气校正参考点,提出了以海面-遥感器的光谱反射率斜率传递现象为基础的斜率算法;汪小钦等[3]利用两个时相的 TM 影像,用线性光谱混合分析法提取福建闽江口的悬浮物浓度;李四海等[4]在 2000年根据多时相NOAA/AVHRR卫星遥感数据和准同步实测表层含沙量资料,利用斜率法、灰度法和泥沙指数法三种模式分别建立泥沙遥感定量模式,并对其效果和适用性进行了比较。近年来,高光谱影像也开始应用到悬浮泥沙监测领域中来,Hyperion影像兼备有波段资源丰富和空间分辨率高的特点,在水体悬浮泥沙浓度遥感提取方面具有很大的应用潜力[5]。

对河口而言,悬浮泥沙浓度不仅是河口地区港口与航道工程十分关心的问题,同时细颗粒泥沙又是各种营养盐和污染物的重要载体,对河口水质环境研究也相当重要。河口区的悬浮泥沙在径流潮流相互作用、盐淡水混合以及风浪等多种动力因素的作用下,其运动规律极其复杂,时空分布变化快。本研究利用射阳河口实测光谱及模拟Hyperion光谱数据与表层泥沙浓度进行相关性分析,构建定量模型,从而实现了射阳河口水体悬浮泥沙浓度的定量反演。

1 数据获取与处理

1.1 研究区概况

射阳河口位于江苏沿海侵蚀性岸段与淤长性岸段之间,其河口北侧属于侵蚀性岸段,南侧属于淤长性岸段。由于水动力条件复杂,悬浮泥沙质量浓度最大可达到11.08 kg/m3,最小仅为0.035 kg/m3,拦门沙和航道的位置也经常发生变化[6](图1)。

射阳河口作为射阳河的入海口,射阳河排污量大,由 1996年全省海洋调查结果可知,与全省其他入海河口相比,射阳河口水质最差,主要以油类、无机氮磷及重金属污染为主[7]。

1.2 水体光谱测量数据及悬沙质量浓度获取

水体光谱测量主要是为获取离水辐射率、归一化离水辐射率和遥感反射率等遥感参数,进而通过这些参数反演得到水体信息。对于二类水体,表面以上测量法是目前唯一有效的方法[8]。

本研究对于射阳河口的水体光谱测量共进行了两次。试验1测量时间是2007年7月1日,如图1所示,本次试验设置5个测量点(S1,S2,S3,S4,S5),通过GPS定位,S1,S2,S3三个测点以及S3,S4,S5三个测点分别位于一条直线。其中,S1,S2,S3三个测点从9:00至15:00每隔1h连续测量水体光谱值。试验2光谱测量时间是2008年5月29～31日,在射阳河口附近海域进行了三个航次的测量和采样,共在60个测点测量了水体的光谱(图1)。每天测量的时间为9:30～14:30。为了使所建立的模型在射阳河口及其附近海域更具有代表性,第三个航次的取样点选择在离河口较远的海域,以确保样本点不仅仅代表河口,还能代表河口附近海域。

图1 射阳河口研究区2007年、2008年采样点位置Fig.1 Sampling sites in 2007 and 2008

两次试验光谱测量所用光谱仪是ASD公司生产的FieldSpec地物光谱仪,该仪器测定的光谱范围为282～1 090 nm,标准板是经严格定标的反射率为0.3的灰板。在每个样点上按上述观测几何条件分别测取水体、天空光和标准板的DN值各10条,并同时采取水样。同时用 GPS定位并记录当时的风速,以便于对影像进行大气校正。

在测量水体光谱的同时,每个样点(包括试验 1中的 S4,S5样点)同时采取水样,采用烘干法获取各样品悬沙质量浓度。两次试验所获取的悬沙质量浓度的量值范围分别为 0.035～2.15 kg/m3(试验 1)和0.04～3.25 kg/m3(试验 2)。

1.3 Hyperion影像数据获取

本研究选用的Hyperion影像是由美国地质调查局(USGS)处理后生成的L1R产品,其VNIR波段和SWIR波段间的空间错位已经经过纠正[9],成像时间为北京时间2007年7月11日。Hyperion数据采用HDF(Hierarchical Data Format)数据集的形式存储,波段存储格式为 BIL格式,由于后面的预处理程序需要在ENVI软件环境下运行,因此在对数据进行处理前需要将其数据存储格式转化为 ENVI软件能够直接读取的标准格式。

对 Hyperion影像的预处理包括:波段剔除、坏线去除、条纹去除、大气校正和几何校正等。考虑到辐射定标、波段噪声与水体信息提取,Hyperion影像的 242个波段中,我们最后保留了 8～57波段(430～930 nm)共50个波段用于研究。利用像元灰度斜率阈值法可以实现影像各波段的坏线判别。在进行坏线修复以后,影像上还存在另外一种像元值异常情况,即垂直条纹。针对条纹的垂直分布特征,采样全局归一化法GNM(Global Normalization Method),通过像元的列平均值、标准差与波段平均值、标准差之间的差异对像元进行分波段线性化修正,可以消除垂直条纹的影响,利用ENVI IDL编制程序,实现对Hyperion影像的坏线修复和垂直条纹去除。采用ENVI软件自带的基于MODTRAN模型的大气校正模块FLAASH软件来实现Hyperion影像的大气校正。影像的几何精校正采用经过几何精校正的的影像底图进行影像到影像的校正[10]。

2 结果与分析

2.1 射阳河口含沙水体光谱特征分析

通过对射阳河口不同悬沙质量浓度水体的光谱特征分析(图2)可以看出,随着悬沙质量浓度增大,水体光谱反射率有升高的趋势。在350～500 nm之间,含沙水体反射率相对较低;在560～720 nm之间有一个反射峰,当泥沙浓度较小时,其峰值主要在560～610 nm 之间,且峰值反射率较低,当泥沙浓度增大时,其峰移在690～720 nm之间,且峰值反射率较高。此外,在790～820 nm之间,还有一个反射峰,当泥沙浓度较低时,该峰值明显低于第一个反射峰,随着泥沙浓度的上升,该峰值迅速上升,并接近第一个反射峰的值;850 nm之后反射率趋于下降。

图2 不同悬沙质量浓度水体的反射光谱曲线Fig.2 Spectra of water samples with different suspended sediment contents

图3 光谱反射率与悬浮泥沙含量相关性分析Fig.3 Correlation between the ASD-based reflectance and suspended sediment content

对试验2中60个测点的悬浮泥沙浓度与对应的实测光谱反射率值的进行相关分析,结果如图3所示,由于在小于350 nm和大于910 nm的波段受噪声干扰较大,所以此分析只取350～910 nm之间的数据进行。结果表明,在波长350～500 nm,悬沙质量浓度与光谱反射率呈现负相关,且相关系数小于 0.2;波长大于500 nm开始,悬沙质量浓度与光谱反射率的正相关性逐渐增大,最大相关系性出现在波长 898～904 nm间,最大相关系数值为0.8858。

2.2 基于Hyperion高光谱数据的射阳河口悬浮泥沙浓度定量模型构建

模型构建前,有必要对实测的光谱数据按光谱响应函数进行光谱重采样,以模拟传感器的光谱响应特征。由于未能从NASA获得Hyperion的光谱响应函数,而 Hyperion是高光谱影像,各波段半高宽(FWHM)很窄,可以用高斯函数模拟它的光谱响应函数[11]。将模拟得到的Hyperion各波段反射率及波段组合因子分别与悬沙质量浓度作相关性分析,通过相关性分析寻找敏感波段或波段组合来构建射阳河口悬沙质量浓度的定量反演模型。已有研究表明,一阶或二阶微分因子虽能体现高光谱数据的优势,但这一因子对传感器的灵敏度和信噪比的要求很高,对大气校正等影像预处理的要求也十分苛刻,如果通过影像不能获得与水面实测光谱完全一致的光谱曲线,则用微分因子构建的模型很难在影像反演中获得令人满意的结果[11]。因此,本研究不再考虑用微分因子作为敏感波段构建反演模型。

由图4所示,Hyperion光谱反射率与悬沙质量浓度的相关性和实测光谱与悬沙质量浓度的相关性变化趋势基本相同,当波长大于 500 nm后,Hyperion光谱反射率与悬沙质量浓度呈现正相关,且最大相关系数位于中心波长 896 nm(最大相关系数值达到0.8916),这与实测光谱预测悬沙质量浓度的敏感波段吻合。

前人的研究表明,单个波段的反射率难以全面地反映出不同泥沙浓度的光谱信息[4]。因此,对Hyperion的各波段进行了9种波段组合(表1)。通过9种波段组合分别与60个测点的悬浮泥沙浓度进行相关性分析,得出相关系数均不超过 0.8,并没有单波段的最大相关系数高。由此可见,对泥沙质量浓度为0.035～3.25 kg/m3范围的河口水体,896 nm虽然仅代表了单波段信息,但其反射率对水体悬沙质量浓度变化具有更灵敏的反应,更能反映含沙水体的光谱特征。

图4 Hyperion光谱反射率与悬浮泥沙含量相关性分析Fig.4 Correlation between Hyperion reflectance and the suspended sediment content

图5 Hyperion 中心波长 896nm波段光谱反射率与悬浮泥沙含量回归分析Fig.5 Regression of the suspended sediment contents based on the Hyperion reflectance at 896 nm

表1 单波段因子和波段组合因子Tab.1 Factors of single band and bands combination

同时,通过对水体氮磷浓度、重金属浓度与Hyperion光谱反射率的相关性分析,896 nm波段都表现出显著的相关特征[12]。因此,以Hyperion 896 nm波段反射率与悬浮泥沙含量的相关性构建射阳河口悬沙质量浓度的定量模型(图5),如式(1)所示:

其中,CSS指射阳河口悬沙质量浓度(kg/m3)。模型R2值达到 0.799,显著高于各种波段组合构建的射阳河口悬沙质量浓度反演模型。

2.3 模型检验

以2007年进行的试验1中S1,S2和S3三个监测位点的实测悬沙质量浓度及光谱测量数据对上述模型进行检验。3个位点光谱测量数据自中午12:00至下午15:00,每隔1 h测量光谱,并在对应时间取水样,利用烘干法获取悬浮浓度,共获取检验数据12组。实测光谱同样进行重采样以模拟Hyperion影像各波段的反射率。模型检验采用常用的相对根均方差法(relative root mean square error,以E表示)对模拟值与观测值之间的符合度进行统计分析,E值的计算见方程(2)。

其中,Pi和Oi分别为预测值和观测值,为观测值的平均值[13]。如表2所示,以模拟Hyperion光谱数据模拟射阳河口悬沙质量浓度模型检验的E为38.7757%,表明模型具有较好的预测精度。

表2 射阳河口表层悬沙质量浓度预测误差Tab.2 Predicted errors of suspended sediment contents for Sheyang estuary

绘制观测值与预测值之间的1 :1关系图,以直观的展示模拟值与观察值的拟合度和可靠性[14],结果显示模拟值与观测值间具有较好的符合度(图6)。

图6 射阳河口悬浮泥沙质量浓度预测值与实测值比较Fig.6 Modeled and observed values of the suspended sediment contents for Sheyang estuary

2.4 模型反演

将上述射阳河口悬沙质量浓度回归模型应用于悬浮物浓度反演制图,结果如图7所示。

图7 2007年7月11日Hyperion影像射阳河口悬沙质量浓度反演结果Fig.7 Result of model simulation on suspended sediment content by Hyperion image

由图7所示,射阳河口悬浮泥沙浓度由近岸向远海呈现递减的规律,在河口南侧的悬沙质量浓度高于北侧,主要是由于河口南侧属于淤长性海岸。

3 总结

本文利用射阳河口实测光谱及模拟Hyperion光谱数据与表层泥沙浓度进行相关性分析,优选了896 nm波段作为研究水体悬浮物浓度的敏感波段并据此建立了定量反演模型,实现了射阳河口水体悬浮泥沙浓度的遥感制图。尽管已有的研究表明单个波段的反射率难以全面地反映出不同泥沙浓度的光谱信息,但针对射阳河口水域,其悬沙类型为淤泥质粉沙,颗粒较为细腻均一,基本没有明显的分级。悬浮泥沙水体对实测光谱反射率的敏感波段为898～904 nm,对光谱分辨率为10 nm的Hyperion影像敏感波段正好位于896 nm,以此构建的指数模型,相关指数达到 0.89。模型的检验结果表明,相对RMSE值为38.78%,表明单波段896 nm对预测射阳河口的悬沙质量浓度具有良好的精度。

由于模型构建与检验的实测悬沙质量浓度数据均来源于上午落潮后的4～5 h后取样,该时间段与Hyperion影像的成像时间相匹配,此时水流场较稳定,河口水体中的泥沙运移与扩散规律与影像解译结果相似,即:近岸悬沙质量浓度高,远海浓度低;河口南侧由于属于淤长性海岸,悬沙质量浓度比北侧高。

[1]Richard L M,Brent A M.Using MODIS terra 250m imagery to map concentrations of total suspended mattering coastal water[J].Remote Sensing Environment,2004,93:259-266.

[2]李炎,李京.基于海面-遥感器光谱反射率斜率传递现象的悬浮泥沙遥感算法[J].科学通报,1994,44(17):1892-1897.

[3]汪小钦,王钦敏,邬群勇,等.遥感在悬浮物质浓度提取中的应用——以福建闽江口为例[J].遥感学报,2003,7(1):54-57.

[4]李四海,恽才兴.河口表层悬浮泥沙气象卫星遥感定量模式研究[J].遥感学报,2001,5(2):154-160.

[5]Bowers D,Gaffney S,White M.Turbidity in the southern Irish Sea[J].Continental Shelf Research,1998,18(5):487-500.

[6]张忍顺,陆丽云,王艳红.江苏海岸侵蚀过程及其趋势[J].地理研究,2002,21(4):469-478.

[7]江苏省海洋污染基线调查队.江苏省海洋污染基线调查报告[M].南京:河海大学出版社,2001:10-56.

[8]唐军武,田国良,汪小勇,等.水体光谱测量与分析Ⅰ:水面以上测量法[J].遥感学报,2004,8(1):37-44.

[9]谭炳香,李增元,陈尔学,等.EO-1 Hyperion高光谱数据的预处理[J],遥感信息,2005,6:36-41.

[10]张东,张鹰,李欢.海岸带星载高光谱遥感影像预处理方法[J].海洋科学进展,2009,27(1):92-97.

[11]赵祥,梁顺林,刘素红,等.高光谱遥感数据的改正暗目标大气校正方法研究[J].中国科学 D辑,2007,37(12):1653-1659.

[12]Pan Jie,Zhang Ying,Xu Yong.Assessment of water quality in Sheyang Estuary (China) using hyperspectral data[C]//Tan Zhenghua.2010 3rdInternational Congress on Image and Signal Processing Symposium.Yantai:Yantai University,2009:1-10.

[13]Michele R,Nicola L,Zina F.Evaluation and application of the OILCROP –SUN model for sunflower in southern Italy[J].Agricultural Systems,2003,78:17-30.

[14]Snyder R L,Spano D,Cesaraccio C,et al.Determining degree-day thresholds from field observations[J].International Journal of Biometeorology,1999,42(4):177-182.