基于改进黑暗像元法的环境一号CCD数据海岸带区域大气校正研究＊

李 微，郭锡杰，刘 远，2，李媛媛，牟 蒙，刘长发，2

（1.大连海洋大学 海洋科技与环境学院，辽宁 大连116023；2.辽宁省高校近岸带海洋环境科学与技术重点实验室，辽宁 大连116023）

环境一号（HJ－1）卫星是我国自主研发的专门用于环境与灾害监测的新型卫星，具有高时间分辨率、高空间分辨率和覆盖范围广等优势，其不仅为环境与减灾业务运行系统提供重要保障，还成为很多部门日常业务的重要数据源。由于受大气的影响，卫星数据不能真实地反映地物的情况，制约了利用卫星数据进行地面参数反演和地面监测的应用。因此，对卫星数据进行高精度的大气校正是十分必要的，大气校正的精度决定了后续遥感卫星数据地面参数反演的精度。

目前，针对卫星遥感数据的大气校正算法分为2类：1）基于辐射传输原理建立起来的大气校正辐射传输模型，该算法精度较高，包括6S模型、MODTRAN模型、FLAASH模型等，但该类算法模型计算量大，需要很多不易获得的地面同步参数；2）不需要大气和地面的实测数据的黑暗像元法，该方法主要通过卫星数据本身的信息就可完成，校正精度能满足一般的遥感预处理和研究，应用于只有卫星遥感数据，而缺乏必要大气和地面参数数据的情况。

近几年，一些学者针对HJ－1卫星CCD数据的大气校正算法进行了研究。金鑫等［1］利用HJ－1A卫星CCD数据提出一种基于辐射传输模型的大气校正方法，对太湖水域大气校正结果较好。孙长奎等［2］利用查找表方法来实现HJ－1ACCD数据的大气校正，得到了研究区域的高空间分辨率的地表反射率数据。方莉等［3］运用像元纯净指数来提取CCD数据上的纯像元，并由 HJ－1A星和B星的多时相CCD观测数据结合地表双向反射率模型，确定了纯像元的地表反射特性，进而较大程度地提高了HJ－1ACCD数据的精度。同时，部分学者也深入研究了适用于HJ－1卫星的黑暗像元大气校正算法。针对HJ－1ACCD数据，王中挺等［4］将密集植被作为暗像元，提出一种基于密集植被的暗目标法，利用辐射传输模型构建查找表，获取可见光范围的地物反射率进而求算气溶胶光学厚度。刘其悦［5］以HJ－1ACCD数据为例，深入探讨了基于辐射传输模型大气校正和基于黑暗像元法的大气校正在不同下垫面条件下应用的可行性，发现不同气候和土地类型对模型的选择有较大影响。综上所述，不同专家学者针对特定区域利用HJ－1CCD数据分别进行了不同大气校正算法的研究，但针对海岸带区域的研究尚不多见。

因此，本文以海岸带为研究区域，针对HJ－1CCD数据，提出了一种改进的黑暗像元大气校正算法。笔者以双台子河口区域为例，利用黑暗像元法大气校正的原理，首先结合归一化植被指数（Normalized Difference Vegetation Index，NDVI）和改进的归一化差异水体指数（Improved Normal Differential Water Index，INDWI）确定初始黑暗像元区域，再利用区域增长法对该区域进行优化并确定最终黑暗像元，然后计算大气的影响因素程辐射，进而去除该影响获得大气校正后的地表反射率R。最后，将该算法与真实地物光谱以及ENVI软件下的Dark－subtract模块和FLAASH大气校正模块的校正结果进行对比分析。

1 研究区域概况

我国海域面积大，海岸带长，近海岸带有着丰富的自然资源以及生物多样性，因此海岸带的监测和保护显得尤为重要。近海岸区域存在大量水体，而且有的区域存在浓密的植被，如红树林、芦苇或翅碱蓬等，另外，部分区域还会存在小岛阴影，而水体、浓密的植被以及阴影在卫星遥感数据上都是理想的黑暗像元。

本文以辽宁双台子河口为研究区域，该区域位于辽宁省盘锦市境内双台子河口国家级自然保护区河口湿地，双台子河入海口，是中国高纬度地区面积最大的滨海湿地。该地地貌为冲海积平原，海岸带地势低洼且平坦，同时分布着大面积的芦苇沼泽区、翅碱蓬滩涂和浅海海域，海岸带较长。因此该区域存在大片水体和浓密的植被，这为黑暗像元的选取提供了条件。

2 大气校正过程

遥感成像过程中，由于大气的影响遥感影像存在一定的辐射量失真现象，降低了遥感数据的精度，因此获得卫星遥感数据后必须进行大气校正。大气校正的目的是消除大气和光照等因素对地物反射率的影响，以获得地物的真实反射率。

2.1 辐射定标

大气校正前，为了建立遥感传感器的数字量化输出值DN与其所对应视场中辐射亮度值之间的定量关系，首先需要对数据进行辐射定标，将DN值转化为表观辐射度L。本文采用中国资源卫星应用中心2013－02－28公布的方法和定标系数（表1）对HJ－1卫星进行辐射定标［6］，其公式如下：

式中，A为定标系数增益，单位为DN／（W·m－2·sr－1·μm－1）；L0为定标系数偏移量，单位为 W·m－2·sr－1·μm－1。

表1 HJ－1卫星CCD数据的定标系数Table 1 Radiometric calibration coefficient of HJ－1CCD data

2.2 大气校正

本文提出的算法有如下假定：1）待校正区域存在黑暗像元；2）地表为朗伯面反射和大气性质均一；3）大气多次散射辐照作用和邻近像元漫反射作用可以忽略。因此反射率小的黑暗像元在大气的影响下，反射率会相对增加［7］。

地表反射率R可通过式（2）计算：

式中，L为真实地物辐亮度，可由公式（1）得到；LP为程辐射，E0为大气层外太阳光谱辐照度，参见表1；ED为天空光漫射到地表面的光谱辐照度，可利用Moran［8］的方法，即用πLP估算；θ为太阳天顶角；φ为传感器天顶角，均可由遥感图像头文件中获取；Tθ和Tφ为大气透过率，从可见光到短波红外过程中，当θ，φ小于70°时，图像受到大气散射和弱吸收的影响，大气透过率Tθ和Tφ可由以下两式估算［9］：

式中，τ为光学厚度，针对近海岸带大气环境，依据Kaufman［10］提出的简化计算方法，利用式（5）估计得出τ：

式中，λ为波长，单位为μm。

由上可知，去除大气的影响获得真实的地表反射率，关键就是获得程辐射，该值可以通过黑暗像元值来确定。因此准确、快速地提取黑暗像元值是黑暗像元法大气校正的核心内容。本文提出的改进黑暗像元大气校正算法具体步骤如下。

2.2.1 初始黑暗像元区域的确定

笔者利用归一化植被指数（NDVI）、改进的归一化差异水体指数（INDWI）对初始黑暗像元区域进行提取。

1）浓密植被备选区域确定

利用归一化植被指数来确定植被备选区域。NDVI是植被生长状态及植被覆盖的最佳指示因子，NDVI值越大，植被覆盖或长势越好。NDVI表达式如下：

式中，Lband4和Lband3分别为HJ－1CCD数据的近红外波段和红光波段的表观辐射度，其计算公式参见式（1）。

研究区域没有浓密的林地，但是存在成片、长势茂密的苇田和生长密集的翅碱蓬群落。通过大量试验最终确定，当NDVI＞0.5时可以提取出浓密的芦苇和翅碱蓬，并将其作为黑暗像元的备选区域。

2）纯净水体区域确定

通常认为，纯净的水体可以作为黑暗像元。但近岸水体中含有大量泥沙以及叶绿素，属于二类水体，因此如何提取到纯净水体是研究重点。相关研究证明，采用归一化的比值指数运算可以最大程度地抑制植被、城市和山地等造成的噪声，从而达到突出水体信息的目的［11］。本文通过研究提出了改进的归一化差异水体指数（INDWI）组合，计算公式如下：

式中，Lband1，Lband3和Lband4分别为HJ－1CCD数据的蓝波段、红波段和近红外波段的表观辐射度，其计算公式参见式（1）。

通过大量试验最终确定，当INDWI1＞0.1且INDWI2＞0.2时即可认为提取到较纯净的水体，作为黑暗像元的备选区域。

2.2.2 黑暗像元区域优化与提取

区域增长法的基本思想是将具有相似性质的像元整合成为一个区域。本文利用区域增长法对黑暗像元区域进行优化提取，具体算法参见李微等［12］。具体实现过程如下：1）确定初始黑暗像元生长点：即将黑暗像元备选区域中第一个使LP大于零的非零像元值所对应的像元作为生长点，可以是一个或者多个生长点；2）确定生长规则：将每个生长点8邻域中值和标准差以及整幅图像平均值作为特征值，可降低随机噪声的影响，避免中值过大的情况；3）进行多阈值区域增长：采用8邻域增长，当符合生长规则时，将其作为新的生长点并更新中值和标准差，直至无法生长为止。如果区域生长范围过小，则认定原生长点是噪声点，同时为了节省运算时间，当区域增长范围过大时终止运算，并将该范围作为黑暗像元区。

统计每一个波段所有增长区域的像素平均值和增长区域个数，将所有增长区域像素平均值相加，除以增长区域个数，即获得黑暗像元值DW′。

2.2.3 程辐射LP的计算

由于程辐射的影响，在传感器上黑暗像元反射率很小但并不为零，由此可以计算出程辐射。一些研究学者例如 Kaufman等［13］假设黑暗像元的反射率为（0.02±0.01），还有一些学者［8，14］假设黑暗像元的反射率为0.01。本文结合研究区域现状，假设黑暗像元具有1%的地表反射率，因此程辐射LP计算方法如下：

式中，L′为黑暗像元区域的辐亮度，计算公式如式（1）所示。将程辐射LP带入到式（2）中，即可得到真实的地表反射率R。

3 大气校正结果与分析

本文利用成像日期为2013－09－01的HJ－1ACCD2级数据为例，验证算法的有效性。研究区域原始数据真彩色合成影像颜色较暗，图1为2%线性拉伸后研究区域真彩色合成影像，由图可见图像对比度差。

图1 2013－09－01研究区域 HJ－1ACCD 3，2，1波段真彩色合成影像Fig.1 True color composite image of HJ－1ACCD with bands 3，2and 1in study area on September 1，2013

现以该数据的第三波段为例说明黑暗像元区域优化提取的过程（图2）。其中：图2a为利用NDVI和INDWI确定的黑暗像元备选区域，图中白色表示黑暗像元备选区域；图2b为第三波段数据初始黑暗像元生长点的位置，用红色标示；图2c为区域增长优化后黑暗像元最终分布图，用红色标示。

图2 黑暗像元区域提取过程图Fig.2 The schematic drawing of extracted dark pixel areas

图3为3种大气校正算法校正后的真彩色合成图方法，其中图3a和3b分别为ENVI软件下的黑暗像元法（Dark－Subtract）和FLAASH方法的校正结果，图3c为改进的黑暗像元大气校正算法（IMDS）的校正结果。三者与图1对比可见图像对比度均明显提高，其中，IMDS方法目视效果更好。

图3 大气校正以后3，2，1波段合成真彩色影像图Fig.3 The true color image with bands 3，2，1after atmospheric correction

为了验证和定量评价改进的黑暗像元大气校正算法，本文对不同大气校正方法校正前后的研究区域内4种典型地物（芦苇、翅碱蓬、水体和滩涂）光谱曲线以及地面实测光谱进行分析对比，实测光谱利用美国ASD地物光谱仪测得的光谱，按照CCD波段带宽取平均得到，如图4所示。从图中可以看出，相比大气校正前，IMDS、Dark－Subtract法和FLAASH法在一定程度上均可以消除大气的影响，相同地物的光谱曲线趋势大致相同，并且相比FLAASH方法，IMDS方法与地面实测值更接近，效果较好。

由图4a可知，IMDS和FLAASH法得到的芦苇光谱曲线在绿光波段（0.55μm）有一小反射峰，红光波段（0.68μm）有一光合作用吸收谷［15］，曲线趋势相近，而Dark－subtract法在绿光波段没有出现明显的反射峰，效果较差；影像获取时值9月，翅碱蓬群落呈现紫红色，从图4b可见，3种方法获得光谱曲线在红光波段反射率都较高，在绿光波段也没有出现反射峰［16］，但IMDS与FLAASH法更接近；图4c为大气校正前后水体的光谱曲线，由图中可见，IMDS和FLAASH法得到水体光谱曲线在绿光波段出现小的反射峰，之后反射率逐渐下降［17］，Dark－subtract法没有出现明显的反射峰；图4d中IMDS和FLAASH法得到的土壤的光谱曲线均优于Dark－subtract法。总的来说，IMDS更接近基于辐射传输模型的FLAASH法，与真实地物实测光谱特征拟合较好，精度较高。

图4 不同方法大气校正前后典型地物光谱与实测光谱值对比Fig.4 Comparison between in situ－measured spectra and satellite－derived ones by different atmospheric correction methods

地物反射率的正常范围为0～1，而大气校正后会出现负值，我们利用“数据可利用率”表示非负值像元个数所占的比例。IMDS与FLAASH法对HJ－1CCD数据可利用率对比见表2。由表2可知，IMDS能有效的改善这种现象，提高了HJ－1CCD数据的可利用率。

表2 大气校正后IMDS与FLAASH数据可利用率（%）对比Table 2 The usable percentage（%）comparison between IMDS and FLAASH after atmospheric correction

本文在计算大气透过率时，IMDS算法只是考虑了大气分子的衰减（式（5）），没有考虑气溶胶的衰减作用；在估算地表下行天空光辐照度时，也只假设了天空光为完全漫射光，没有考虑天空光存在显著的各向异性。这是该算法的不足之处，也是今后需要改进的重点。

4 结论

本文基于黑暗像元大气校正模型，针对海岸带的地表状况，提出了一种改进的黑暗像元大气校正模型。首先结合归一化植被指数（NDVI）、改进的归一化差异水体指数（INDWI）来确定黑暗像元，同时利用区域增长法对黑暗像元区域进行优化，计算出程辐射，实现了大气校正过程。结果表明，IMDS得到地物反射率曲线相近，与真实地物光谱特征拟合较好，且进行大气校正之后的图像地物边缘清晰度明显增强，图像亮度增加，图像对比度较高，目视效果好。同时与FLAASH法相比，IMDS能改善FLAASH法校正过程中出现负值的现象，可以有效的保证了HJ－1CCD数据的可利用率。

致谢：文中HJ－1A数据来源于国家环境保护部卫星环境应用中心（http：∥www.secmep.cn），在此表示感谢！

（References）：

［1］JIN X，LI Y M，WANG Q，et al.Atmospheric correction of satellite HJ－1CCD image based on aerosol distribution in Lake Taihu［J］.Journal of Lake Sciences，2010，22（4）：504－512.金鑫，李云梅，王桥，等.基于太湖气溶胶类型分区的环境一号卫星 CCD大气校正［J］.湖泊科学，2010，22（4）：504－512.

［2］SUN C K，SUN L，MA S F，et al.Atmospheric correction method based on HJ－1CCD data［J］.Journal of Remote Sensing，2012，16（4）：831－836.孙长奎，孙林，麻盛芳，等.HJ－1CCD数据大气校正方法研究［J］.遥感学报，2012，16（4）：831－836.

［3］FANG L，YU T，GU X F，et al.Aerosol retrieval and atmospheric correction of HJ－1CCD data over Beijing［J］.Journal of Remote Sensing，2013，17（1）：158－164.方莉，余涛，顾行发，等.北京地区 HJ－1卫星CCD数据的气溶胶反演及在大气校正中的应用［J］.遥感学报，2013，17（1）：158－164.

［4］WANG Z T，LI Q，TAO J H，et al.Monitoring of aerosol optical depth over land surface using CCD camera on HJ－1satellite［J］.China Environmental Science，2009，29（9）：902－907.王中挺，厉青，陶金花，等.环境一号卫星 CCD 相机应用于陆地气溶胶的监测［J］.中国环境科学，2009，29（9）：902－907.

［5］LIU Q Y.The atmospheric correction for visible－neainfrared region of HJ－1constellation data［D］.Henan：Henan Polytechnic University，2010.刘其悦.基于环境星数据可见光－近红外波段的大气校正方法研究［D］.河南：河南理工大学，2010.

［6］Chinese Center for Resources Satellite Data and Application.HJ－1AB，ZY－3and 02Cabsolute radiometric calibration coefficient in 2013.［EB／OL］.［2013－09－17］.http：∥www.cresda.com／n16／n1115／n1522／n2103／index.html.中国资源卫星应用中心.2013年 HJ－1AB，ZY－3和02C绝对辐射定标系数［EB／OL］.［2013－09－17］.http：∥www.cresda.com／n16／n1115／n1522／n2103／index.html.

［7］DENG S B.ENVI RS image processing methods［M］.Beijing：Science Press，2010.邓书斌.ENVI遥感图像处理方法［M］.北京：科学出版社，2010.

［8］MORAN M S，JACKSON R D，SLATER P.Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output［J］.Remote Sensing of Environment，1992，41（2）：169－184.

［9］ZHENG W，ZENG Z Y.Dark－object methods for atmospheric correction of remote sensing image［J］.Remote Sensing for Land ＆ Resources，2005，63（3）：8－11.郑伟，曾志远.遥感图像大气校正的黑暗像元法［J］.国土资源遥感，2005，63（3）：8－11.

［10］KAUFMAN Y J，WALD A E，REMER L A，et al.The MODIS 2.1μm channel－correlation with visible reflectance for use in remote sensing of aerosol［J］.IEEE Transactions on Geoscience and Remote Sensing，1997，35（5）：1286－1296.

［11］XU H Q.A study on information extraction of water body with the modified normalized difference water index（MNDWI）［J］.Journal of Remote Sensing，2005，9（5）：589－595.徐涵秋.利用改进的归一化差异水体指数（GNDWI）提取水体信息的研究［J］.遥感学报，2005，9（5）：589－595.

［12］LI W，TIAN Y，LIU Y，et al.The algorithm to extract automatically dark objects in coastal zone based on region growing［J］.Journal of Dalian Ocean University，2013，28（5）：502－505.李微，田彦，刘远，等.基于区域增长的海岸带黑暗像元自动提取算法研究［J］.大连海洋大学学报，2013，28（5）：502－505.

［13］KAUFMAN Y J，SENDRA C.Algorithm for automatic atmospheric corrections to visible and near－IR satellite imagery［J］.International Journal of Remote Sensing，1988，9（8）：1357－1381.

［14］CHAVEZ Jr P S.Image－based atmospheric correction revisited and improved［J］.Photogrammetric Engineering and Remote Sensing，1996，62：1025－1035.

［15］MEI A X，PENG W L，QIN Q M，et al.An introduction to remote sensing［M］.Beijing：High Education Press，2001.梅安新，彭望琭，秦其明，等.遥感导论［M］.北京：高等教育出版社，2001.

［16］WU T，ZHAO D Z，KANG J C，et al.Research on remote sensing inversion biomass method based on the Suaeda Salsa＇s Measured Spectrum［J］.Spectroscopy and Spectral Analysis，2010，30（5）：1336－1341.吴涛，赵冬至，康建成，等.基于现场光谱的翅碱蓬生物量遥感反演方法研究［J］.光谱学与光谱分析，2010，30（5）：1336－1341.

［17］ZHANG Y，ZHANG Y，WANG J J.Correlations between the spectral reflectance of water mass containing suspended sediments and the concentration and particle size of the suspended sediments［J］.Advances in Marine Science，2008，26（3）：340－346.张芸，张鹰，王晶晶.悬沙水体光谱反射率与质量浓度、粒径的相关关系［J］.海洋科学进展，2008，26（3）：340－346.