海岸线(Shoreline)是多年平均大潮高潮所形成的海水和陆地分界线的痕迹线^[1-2]，是划分海洋与陆地行政管理区域的基准线，也是国际地理数据委员会(International Geographic Data Committee, IGDC)认定的27个地表要素之一^[3]。海岸线的位置及其变化是安全导航、城镇扩建和沿海资源规划和管理的基础^[4-5]，而且随着卫星遥感技术的进步，中高光谱数据成本会大大下降并广泛应用到滨海地物的提取分离中。因此，及时的提取和监测海岸线的变化对海岸领域管理具有重要意义。

遥感在海岸线提取中占有重要的作用，遥感数据是提取海岸线的重要数据源。本文使用遥感数据提取的是卫星过顶时拍摄的瞬时水陆分界线亦即海岸瞬时水边线^[6-7]，常用于海岸瞬时水边线提取的遥感图像主要分为SAR图像、LiDAR图像和中高光谱图像^[8-9]。SAR图像提取海岸瞬时水边线具有宽覆盖、高分辨率以及全天候的特点，但是由于SAR图像容易受到散斑效应的影响，导致成像时会有众多的噪点，因此对海岸瞬时水边线的提取精度造成影响^[10-11]。对于LiDAR可以快速获得大范围的图像，但是这种数据一般受到成本的限制，在对多幅图像进行处理时，成本较高^[8]。对于多光谱图像可以获得覆盖面积较大的区域，但是由于受到波段数目的限制而产生光谱信息较少的离散信号，水、植被、不透水层和土壤的光谱特征难以精确地体现出来，会造成海岸地物误分的现象，影响海岸瞬时水边线提取的精度^[12-13]。然而，考虑到海岸地带的环境复杂，目标种类繁多，高光谱图像包含连续的光谱信息和丰富的海岸地表特征信息，能够以丰富的光谱信息采集捕捉到水、植被、不透水层和土壤的特征信息，具有较高的能力区分海水和其他海岸地物。但是，若采用多光谱图像(如Landsat)，由于包含较少的海岸地物的光谱信息，会对水、植被、不透水层和土壤地物的真实光谱造成干扰，影响海岸瞬时水边线的提取精度。因此使用高光谱图像获得丰富的海岸地物光谱信息能够提高海岸瞬时水边线提取精度。

从中高光谱图像中提取海岸瞬时水边线的方法有基于像元级和亚像元级的方法。像元级方法是以像元作为地物类别划分的基本单位，将像元中所占丰度最大的地物类别作为像元类别，由于该方法忽略了像元中丰度较小的地物，使得该类端元没有识别出来，影响了边界线的位置，限制了海岸瞬时水边线提取的精度^[14-15]。然而，基于亚像元级方法在一定程度上可以解决海岸边缘混合像元问题，其允许一个像元含有多个端元类别，可以给出精确可靠的水、植被、不透水层和土壤端元代表，获得比像元级方法精度较高的海岸瞬时水边线。但是已存的亚像元级方法对于复杂海岸地带的混合像元的提取和端元的确定并不适合^[16]；而且，在光谱分解时通常只使用一个混合像元的端元丰度值之和为1或者是端元丰度值为非负的约束条件，由于低估或高估混合像元的水体丰度值而导致提取海岸瞬时水边线精度较差^[17]。除此之外，很少有学者使用亚像元制图方法获得不同水体丰度具体的空间位置，因此在一定程度上难免会对海岸瞬时水边线提取精度造成影响。而且这些方法一般应用在湖面边线和城市河流边线的提取^[18-19]，但将其用于泥沙较多的复杂海岸提取水边线的研究较少。

本文提出一种改进的亚像元分解方法对上海长江入海口实验区域进行海岸瞬时水边线提取，主要目标为：(1)使用了一种水体-植被-不透水层-土壤模型(W-V-I-S)用于检测和确定海岸地区的混合像元和纯净端元光谱；(2)使用全约束最小二乘法(FCLS)估计混合像元中水体丰度值；(3)使用空间吸引力模型将丰度图像转换成地物分类图像获得最终的海岸瞬时水边线。

1 研究区域及数据介绍

我国海岸类型大致可分为基岩海岸、砂质海岸、淤泥海岸以及生物海岸^[20]，而长江入海口属于淤泥海岸。但是近几年长江口地带出现将大量的滩涂湿地转为城市建设的现象，不但不利于长江口湿地生态系统的保护，而且使得长江口海岸每年都发生着明显的变化。因此，为了提供更好的海岸管理，本次研究选取该地段进行研究。

由于在长江入海口海岸一般是淤泥海岸和人工海岸，根据该海岸地区的海水特性，由于图像的光谱信息包含较少，难以区分海岸水、植被、不透水层和土壤地物，影响海岸瞬时水边线的提取精度，因此采用高光谱图像对海岸地物进行处理。如图 1所示，其中红色矩形框区域是本文研究区域，采用了2006年11月8日的EO-1卫星搭载的Hyperion传感器获取的高光谱遥感影像作为实验数据，中心经度为121.393 2°E，中心纬度为30.657 4°N，该数据是经过几何校正、投影校正和地形校正的L1Gst级数据。该数据是总共包含242个波段，光谱范围为357~2 576 nm，空间分辨率为30 m，但是由于高光谱数据波段1~7，58~76和225~242未经过辐射校正，被记录为坏波段，其值设为0 ^[21]，因此在高光谱数据处理时将不考虑这些波段。

2 改进的亚像元分解高光谱海岸瞬时水边线提取方法

图 2是基于ISPCW方法提取海岸瞬时水边线的总体框架图，该方法主要研究了以下3个问题：首先，使用W-V-I-S模型获得研究区域的混合像元和端元光谱；其次，使用线性光谱混合模型估计混合像元中的端元水体丰度图；最后，使用空间吸引力模型获得海岸瞬时水边线。该方法将在下面几节中详细讨论。

2.1 基于W-V-I-S模型的混合像元提取和端元确定

长江入海口的地表特征主要包括水、植被、不透水层和土壤等^[22]，根据这些地物特征，水体-植被-不透水层-土壤模型(Water-Vegetation-Impervious-Soil, W-V-I-S)使用归一化水体指数(Normalized Difference Water Index, NDWI)、归一化植被指数(Normalized Difference Vegetation Index, NDVI)、归一化建筑物指数(Normalized Difference Built-Up Index, NDBI)和归一化土壤指数(Normalized Difference Soil Index, NDSI)四种指数被建立。该模型可以获得W-V-I-S混合像元和地表水、植被、不透水层和土壤特征的参考光谱，可以为海岸地带的混合像元的分解提供有效的指导^[23]。表 1是NDWI、NDVI、NDBI和NDSI四种指数的计算公式^{[15, 24-25]}，其中公式中ρ(Green)、ρ(NIR)、ρ(R)、ρ(SWIR)和ρ(Yellow)分别是绿波段、近红外波段、红波段、短波红外波段和黄色波段的反射率。

为了自动提取W-V-I-S混合像元及其对应的纯净像元，首先在NDWI、NDVI、NDBI和NDSI直方图中，使用局部加权散点图平滑法对四种指数的直方图进行平滑，然后W-V-I-S模型以直方图谷底阈值作为起始点分别向两边搜索W-V-I-S混合像元。本文以直方图曲线斜率作为判别标准，根据经验值通常被设置为30°、45°和60°^[19]，然后根据不同的曲线斜率分别搜索W-V-I-S混合像元，判断在30°、45°和60°曲线斜率下混合像元的提取情况，选取能够最大程度包含海岸边缘的W-V-I-S混合像元的斜率作为该地物的判别标准。由于在本研究区域中，相对于纯净像元而言，混合像元更容易与复杂海岸地表特性相混淆，根据这一特性，纯净水、植被、不透水层和土壤的阈值分别被选为60°、60°、30°和30° ^[19]。当首次满足这个条件时对应横坐标就是纯净水、植被、不透水层和土壤的阈值，对应的W-V-I-S混合像元即可获得。

本文使用像元纯净指数(Pixel Purity Index, PPI)算法^[26]对海岸地带的W-V-I-S混合像元进行纯净端元提取，同时为了确定PPI算法提取的端元类别，考虑到海岸光谱的易变性，使用光谱库中的光谱作为参考光谱很难得到较好的检测效果。因此，为了获得更真实可靠的参考光谱，根据海岸地带的地物特征，使用NDWI、NDVI、NDBI和NDSI灰度值较大的若干个像素点取其平均值作为参考光谱，能够更精确的确定本研究区域中提取的端元光谱的类型。然后采用基于光谱角匹配方法对端元光谱和W-V-I-S模型获得的参考光谱进行匹配，光谱角越小，表明与该地物越匹配^[27]。

2.2 基于空间吸引力模型的高光谱海岸瞬时水边线提取

全约束最小二乘法(Fully Constrained Least Squares，FCLS)用来估计海岸水体丰度值，该方法包含两个约束条件，即一个像元内端元的丰度值之和为1且每个端元的丰度值必须在0和1之间。W-V-I-S模型假设海岸地表特征，包括水、植被、不透水层和土壤是线性光谱混合模型，即图像中的每个像元光谱都可以由所有端元光谱线性组合近似地表示。公式(1)代表的是在两个约束条件下通过多次迭代之后的最小残差，公式(2)代表的是两个约束条件。基于公式(1)和公式(2)，海岸地表的水、植被、不透水层和土壤在混合像元中的丰度值分别得到，则线性光谱混合模型可以描述为^[19]：

式中: p_w、p_v、p_i、p_s、和p_m分别是纯净水、植被、不透水层、土壤和混合像元的反射率，f_w、f_v、f_i、和f_s分别代表的端元水、植被、不透水层和土壤的丰度值。

为了获得不同丰度值在混合像元的空间分布情况，亚像元制图方法可以精确的定位到端元在混合像元中实际的空间位置。由于海岸瞬时水边线是连续的，空间吸引力模型能够更好的表达海岸水、植被、不透水层和土壤之间的空间相关性，获得端元在混合像元中具体的空间位置，从而提高海岸瞬时水边线提取精度。

在本研究中，空间吸引力模型的比例尺度S通过水体丰度直方图确定。为了方便计算像元和亚像元之间的距离，将其中心坐标作为该像元的起始点，一个单元就是一个亚像元大小。为了得到不同类别的吸引力值，将中心像元的邻域像元中不同类别的丰度值用于计算，本研究中选用八邻域类型。公式(4)计算得到亚像元与像元之间的距离并进而通过公式(3)计算得到各个类别对每个亚像元的吸引力值，其具体公式如下：

式中：p_{a, b}表示位置为(a, b)的亚像元，p_{a, b} (c)表示类别c对亚像元p_{a, b}的吸引力值，P_{i, j}表示位置为(i, j)的像元，P_{i, j} (c)表示像元P_{i, j}中类别c的丰度值，d(p_{a, b}, P_{i, j})表示亚像元p_{a, b}与像元P_{i, j}之间的距离，N_t (p_{a, b})表示选择的领域类型，即八邻域^[28]。

为了验证本文提出的ISPCW方法提取海岸瞬时水边线的精度，考虑到实验数据采用的是2006年11月8日数据，为了获取更真实可靠的基线，选用了2006年11月12日的高空间分辨率(优于1 m)Google Earth^TM图像获得的数字化水边线作为基线，由于实验区域的空间分辨率是30 m，而Google Earth^TM的高分辨率数据优于1 m，因此可以作为基线对提取的海岸瞬时水边线的精度进行评估。在计算使用多种方法提取的海岸瞬时水边线与数字化水边线之间的距离中，通过使用垂直断面法对海岸瞬时水边线的变化进行统计，该方法是由美国地质调查局开发的ArcGIS的扩展模块DSAS(Digital Shoreline Analysis System)^[29]。由于高光谱图像的空间分辨率是30 m，在进行断面间隔设置时，分别采用了30 m、60 m、90 m和120 m对本文提出的方法和对比方法提取的海岸瞬时水边线进行断面，并分别计算出提取的海岸瞬时水边线到数字化水边线之间的垂直距离，计算出均方根误差(Root Mean Square Error，RMSE)^[20]。RMSE值越小，表明该方法提取的海岸瞬时水边线与实际岸线越接近，效果越好。

3 结果与分析 3.1 基于W-V-I-S模型的端元提取和确定

根据2.1节介绍的W-V-I-S模型提取的混合像元，首先使用最小噪声分离变换(Minimum Noise Fraction, MNF)对原始图像进行数据降维，去除图像噪声，然后获得MNF变换后的前两个波段W-V-I-S混合像元的主要成分散点图，如图 3 (a)所示，横坐标是第一个波段的主要成分，纵坐标是第二个波段的主要成分；考虑到PPI是利用端元是遥感图像在特征空间中形成的单形体的端点的特点及单形体的向量投影性质进行端元提取^[30]。因此，根据图 3 (a)可以获得4个明显的凸集端点A~D，每个端点代表着W-V-I-S混合像元中的一种纯净像元，即分别对应着4个海岸边缘的水、植被、不透水层和土壤纯净端元。图 3 (b)是提取的纯净像元A~D对应的端元光谱结果，以及根据W-V-I-S模型获得的水、植被、不透水层和土壤的参考光谱曲线，从图中可以看出，使用该模型建立的光谱库与端元光谱具有较高的匹配结果。

为了进一步确定端元光谱的类别，使用光谱角匹配方法计算端元光谱和参考光谱之间的夹角，光谱角越小，表明与该地物越匹配。如表 2显示了各个端元光谱与参考光谱之间的夹角，与参考水、植被、不透水层和土壤进行相比，端元类别A~D的夹角最小值分别是0.0660 rad、0.0220 rad、0.0585 rad和0.0350 rad，即端元A~D的地物类别分别对应的是水、植被、不透水层和土壤。从实验结果来看，尽管不透水层和土壤的光谱特征是相似的，但是具体到本实验研究的长江入海口区域，不透水层和土壤还是很明显的从复杂的海岸地带区分开。

3.2 基于ISPCW方法的高光谱海岸瞬时水边线提取

根据全约束最小二乘法(FCLS)线性混合模型，可以获得水体、植被、不透水层和土壤丰度图，为了获得各个端元的不同丰度值在混合像元中的实际空间分布情况，亚像元制图是必要的。根据2.2节中介绍的空间吸引力模型可以获得像元中不同地物端元之间的空间相关性。根据水体丰度值的直方图，可以确定空间吸引力模型中的比例尺度S是3，即可以将一个像元划分为9个亚像元，通过实验结果表明，该比例尺度可以精确的获得不同丰度值的海岸地表特征的空间分布情况。

基于ISPCW方法提取海岸瞬时水边线的结果如图 4所示，其中图 4(a)是原始图像，根据公式(3)计算出亚像元和邻近像元之间的吸引力值，获得最终的亚像元制图结果，如图 4(b)所示。与原始图像相比，从实验结果中可以明显的看出提取的含有不同丰度值的端元水、植被、不透水层和土壤在混合像元中具体的空间分布情况，特别在海岸边缘地带，水和土壤的混合程度比较严重，表明在研究区域内海岸边缘淤泥比较严重。图 4(c)是本文提出的ISPCW方法在亚像元制图的基础上提取的海岸瞬时水边线结果，从图中可以看出，ISPCW方法提取的海岸瞬时水边线能够更好的体现出淤泥海岸真实的水陆分界线。

3.3 基于多种方法的高光谱海岸瞬时水边线提取

为了验证本文提出的ISPCW方法提取海岸瞬时水边线的精度，采用了多端元光谱混合分析(Multiple Endmember Spectral Mixture Analysis, MESMA)、混合调谐匹配滤波法(Mixture Tuned Matched Filtering, MTMF)、连续最大角凸锥(Sequential Maximum Angle C0onvex Cone, SMACC)、能量约束最小化(Constrained Energy Minimization, CEM)4种混合像元分解方法和归一化水体指数(Normalized Difference Water Index, NDWI)的像元级方法作为对比方法，以数字化水边线为基线，对本文提出的ISPCW方法提取海岸瞬时水边线的精度进行评估。

如图 5所示，图 5 (a)是原始数据，图 5 (b)~(h)分别是数字化水边线和ISPCW、MESMA、MTMF、SMACC、CEM、NDWI方法提取的海岸瞬时水边线结果图。从图中可以看出，本文提出的ISPCW方法提取海岸瞬时水边线的效果优于其他亚像元级的方法和像元级的方法。在对比方法中，SMACC、CEM和NDWI方法提取的海岸瞬时水边线明显偏离数字化水边线，而MESMA和MTMF方法与以上三种方法相比在提取海岸瞬时水边线上具有较好的性能，但是具体到更小的细节区域，这两种方法提取的海岸瞬时水边线结果含有相应的误差。对于长江入海口海岸地带，本文使用的ISPCW方法充分利用了水体光谱信息，可以精确的提取在混合像元中的水体端元，有效的避免由于漏分和错分像元产生的误差。因此，相比较而言，本文使用的ISPCW方法可以提供较准确的海岸瞬时水边线的位置。

3.4 精度验证

本次实验利用2006年11月12日的高空间分辨率(优于1 m)Google Earth^TM图像获得的数字化水边线作为基线，对海岸瞬时水边线提取结果进行精度评估。图 6是基于ISPCW、MESMA、MTMF、SMACC、CEM和NDWI方法的在断面间距分别为30 m、60 m、90 m和120 m时计算出的RMSE值，RMSE值越小，表明该方法提取海岸瞬时水边线的效果较好。从图中可以看出，本文提出的ISPCW方法提取海岸瞬时水边线的效果明显优于其他对比方法。对于不同的断面间距而言，ISPCW方法提取海岸瞬时水边线效果最好的是在30 m，即RMSE值为11.4 m，效果最差的是在120 m，即RMSE值为12.0 m，同样，在不同的断面间距提取海岸瞬时水边线效果中，对比方法MESMA、MTMF、SMACC、CEM和NDWI方法有着相同的趋势，即随着断面间距的增大RMSE值也会增大。因此，断面间距的选取也会影响海岸瞬时水边线提取精度，但是本文提出的方法在各个断面间距上提取海岸瞬时水边线的精度都优于其他方法。而且与数字化水边线相比，断面间距越小可以更清晰的表达出提取海岸瞬时水边线与实际岸线之间的差距，而且考虑到使用的图像的空间分辨率是30 m，因此，与其他断面间距提取的海岸瞬时水边线相比，使用30 m的断面间距提取海岸瞬时水边线可以获得较高精度。

在30 m的断面间距中，ISPCW、MESMA、MTMF、SMACC、CEM和NDWI方法的RMSE值如图 7所示，其中MESMA、MTMF、SMACC、CEM和NDWI方法提取海岸瞬时水边线的精度分别是0.49个像元、0.57个像元、0.66个像元、0.67个像元和0.78个像元，而本文提出的ISPCW方法提取海岸瞬时水边线的RMSE值是11.4 m，即精度为0.38个像元，与对比方法相比，精度分别提高了22.4%、33.3%、42.4%、43.2%和51.3%。虽然本文采用的高光谱图像的空间分辨率有所限制，但是长江入海口实验区域包含大范围的水，其水的纯净像元很容易找到，与原始数据相对照，选取水的精度会很高，可以充分利用丰富的水体光谱信息获得比较好的海岸瞬时水边线提取精度。因此，实验结果表明，本文提出的ISPCW方法使用高光谱图像提取海岸瞬时水边线获得较高的精度。

这是因为CEM方法对阴影部分比较敏感，由于海岸地带含有大量的阴影部分，而且海岸地带的水和阴影的光谱信息是相似的，因此使用该方法难以将这两种地物区分开^[31]。对于像元级的NDWI方法，图像中的每个像元被硬性的分为在像元中所占比例比较大的地物，限制了水边缘的提取；虽然NDWI阈值分割法可以快速的实现水体提取，但是NDWI最佳阈值一般很难确定，而且只使用到水体的2个特征光谱波段，精度上没有充分保证^[32]。因此，相对于ISPCW方法，这些方法提取海岸瞬时水边线的效果较差。尽管MESMA、MTMF和SMACC在提取海岸瞬时水边线上比CEM和NDWI方法性能更好，但是由于在光谱分解时只满足一个约束条件，即端元在混合像元中的丰度值之和为1或者每个端元的丰度值在0和1之间，因此在混合像元分解时由于低估或者高估端元的丰度值而导致提取海岸瞬时水边线的精度较差^[17]。

在对比方法中，MESMA获得较好的海岸瞬时水边线提取精度，但是在复杂的海岸地带，当端元光谱较多时端元组合情况较为复杂，像元数目庞大导致耗时较长，不利于大区域多端元光谱库的分解。因此，根据以上精度统计分析来看，与MESMA、MTMF、SMACC、CEM和NDWI方法相比，本文提取的ISPCW方法在海岸瞬时水边线提取上性能最好，达到亚像元级的精度。

4 结论

本文提出了一种改进的亚像元分解方法对上海长江入海口区域进行海岸瞬时水边线提取，并与传统的MESMA、MTMF、SMACC、CEM和NDWI方法在提取海岸瞬时水边线精度方面进行了对比分析。

实验结果表明，断面间距选取的不同会影响海岸瞬时水边线提取精度，但是本文提出的ISPCW方法在各个断面间距提取海岸瞬时水边线的精度都优于其他方法，其精度达到0.38个像元，与MESMA、MTMF、SMACC、CEM和NDWI方法相比，精度分别提高了22.4%、33.3%、42.4%、43.2%和51.3%。

因此，从定性角度和定量角度来看，本文提出的ISPCW方法在总体上获得了较优的海岸瞬时水边线提取性能，可以精确的分解复杂海岸地带的混合像元，与传统的像元级和亚像元级方法相比，可以获得较高的海岸瞬时水边线提取精度。

感谢德国宇航中心(DLR)提供的免费数据；感谢上海海洋大学信息学院通信导航实验室老师及同学对本文中实验的帮助。