尹航,戚洪帅, ,蔡锋, ,张弛,刘根,赵绍华, ,宋嘉诚, ,赵国润,
( 1. 自然资源部第三海洋研究所 海洋与海岸地质实验室, 福建 厦门 361005;2. 福建省海洋生态保护与修复重点实验室,福建 厦门 361005;3. 河海大学 港口海岸与近海工程学院,江苏 南京 210098)
海岸带是陆地与海洋之间的交界地带,是地球表面最为活跃的自然区域,也是社会经济领域中的“黄金地带”[1]。近年来,随着全球海平面上升[2]、河流入海泥沙锐减以及沿海开发加剧,海岸侵蚀、海岸人工化等导致海岸带快速变化[3-4]。而海岸线的位置变迁是海岸带演化最直观的体现[5],遥感手段能够快速获取大范围海岸线数据,对海岸带研究和保护管理皆具有十分重要的意义[6]。
高分辨率遥感影像是进行海岸线提取和精细化制图的重要数据源[7],依此提取精细海岸线数据是当前海岸带遥感研究的热点。遥感岸线提取的核心是对水陆边缘像素识别与定位,现有的自动提取方法主要包括阈值分割法[8-9]、边缘检测法[10-11]、面向对象法[12-13]、区域生长法[14-15]等,其中以边缘检测法最为成熟且应用最为广泛[16]。目前,国内外学者已使用多种边缘检测算子法开展遥感影像岸线提取研究,包括Sobel算子法[17]、Roberts算子法[18-19]、Canny算子法[20-22]、高斯-拉普拉斯(LoG)算子法[10]等。针对砂质海岸线提取,Kass等[23]提出了一种基于Snakes活动轮廓模型的新型边缘检测算法,可应用于砂质岸线的提取和平滑处理,但由于考虑参数过多致使模型较为复杂[16];De Laurentiis等[24]提出了基于图像纹理特征和模糊连通性的砂质海岸水边线检测定位方法;García-Rubio等[25]使用ISODATA非监督分类法对多光谱影像进行海陆分离,并对水边线进行检测提取和平滑处理;张旭凯等[26]使用改进归一化差异水体指数(Modified Normalized Difference Water Index, MNDWI)[27]分 离水体和陆地,结合Canny算子提取得到水边线,并进行了潮位校正。上述方法尚存在噪声敏感性、阈值不稳定性等问题,尤其在处理空间分辨率高、背景复杂的海岸影像时,易出现噪点多、伪边缘、岸线不连续等问题,严重影响岸线提取的效率和精度[6,16]。
由于遥感影像提取的岸线是卫星过境时刻的瞬时水边线,岸滩坡度较缓的砂质岸线易受潮位影响而出现较大的位置偏移[28]。针对这一问题,Liu和Jezek[11]及马小峰等[29]提出了基于岸滩坡度的线性潮位校正模型,Stockdonf等[30]及刘善伟等[7]则通过获取岸滩数字高程模型(Digital Elevation Model, DEM)数据进行岸线优化与校正。然而现实中,砂质海岸海滩剖面坡度并非线性,通常是呈上凹状的平衡剖面形态[31]。而岸滩DEM数据获取复杂且不同时期存在变化,不具备普适性。
本文利用高分遥感影像开展复杂背景下砂质海岸线精细化提取研究,针对上述岸线提取与校正方法存在的不足,引入一种高精度、强噪声鲁棒性的结构森林边缘检测(Strected Forests Edge Detection, SE)算法[32]识别水边线位置,并提出建立拟合剖面模型进行潮位校正得到海岸线,克服了传统线性潮位校正模型误差偏大的问题。
SE算法是一种基于结构化学习模型的新型边缘检测算法,相比于传统边缘检测算法,其运行速度更快、噪声鲁棒性更强、检测精度更高[32]。
SE算法的核心是对随机森林(Random Forest,RF)模型的训练,RF模型是由n个不相关决策树组成的分类器(图1),每棵决策树由图像数据集S中随机排列的数据训练,最终预测结果通过决策树的综合投票得到[33],该模型具有抗干扰能力强、抗过拟合能力强、预测精度高等突出优点[32-34]。
图1 随机森林模型示意图Fig. 1 Schematic diagram of random forest model
使用单个决策树对给定图像数据集S∈X×Y进行训练时,对x∈X(x表示图像块),依据量化特征,得到相应的分类结果y∈Y(y表示标签)。对于每个节点j,其分离函数h定义为
式中,θj为增益参数;k是x的某一量化特征;τ是该特征的阈值。如果h为0,x归类为节点j左侧,如果h为1,则x归类为节点j右侧。
为使分离函数h性能最大化,定义标准信息增益,依据信息增益最大化原则确定,执行递归训练:
式中,H(S)为香农熵;py为训练数据集中图像块x出现标签y的概率;SL用于训练左节点;SR用于训练右节点,当信息增益达到阈值时,训练停止。
结构化学习是对映射问题的探索,其输出结果可以是图像、语句等复杂表示。SE算法通过将RF模型扩展到结构化输出空间,把边缘检测问题公式化为在给定输入图像块的情况下预测局部分割掩膜[32,35]。对于给定节点j,将该节点上的所有标签y(y∈Y)离散化映射到离散化标签c(c∈C),便可用c代替y计算Ij:
式中,Y为信息增益未良好定义的结构化空间;C为信息增益定义良好的离散空间;Z为使相似性度量易于计算而定义的中间空间。
以上RF模型及SE算法构建步骤详见Breiman[34]及Dollár和Zitnick的研究[32]。
使用SE算法从3个颜色、2个梯度幅值和8个梯度方向,共13个通道扩充输入影像的每个图像块,检测边缘特征。对于32×32结构化图像块,使用16×16的窗口预测分割掩膜,取各个掩膜在单一像素点的平均值作为最终像素值,输出结果是1幅灰度图像(图2),各像素点灰度值代表其为边缘点的概率。
图2 结构森林边缘检测算法边缘概率Fig. 2 Edge probability of strected forests edge detection
为得到数字化水边线,依据输出结果设定阈值,将边缘概率图像转化为0和1的二值图像,在保留边缘检测效果的同时便于进行形态学处理。应用Arc-Scan自动矢量化方法,从栅格数据中提取边缘线,经过后期编辑处理即可得到数字化水边线结果。
潮位校正是依据海滩地形和潮位变化关系,将遥感提取瞬时水边线校正至一般意义海岸线(平均大潮高潮线)的过程,是基于遥感影像进行岸线精细提取的重要环节,对于砂质海岸的岸线提取尤为关键。常用的线性潮位校正模型原理如图3所示:假设海滩剖面为均一的斜坡,基于不同潮时影像水边线c1、c2推算岸滩坡度θ,计算校正距离d[26]。该模型主要缺点有两个:(1)忽略了岸滩地形、地貌因素影响,校正距离与实际偏差过大;(2)坡度数据多来源于不同潮位时刻水边线推算,容易引入遥感预处理误差、解译误差等新误差。
图3 线性潮位校正模型原理(修改自文献[27])Fig. 3 The principle of linear tide level correction model(modified from reference[27])
自然状态下,砂质海岸的海滩通常发育平衡剖面,一般呈上凹形态,坡度随向海距离增加而变缓,可使用数学模型描述其形态[31]。Bruun[36]最早提出海滩平衡剖面模式,Dean[37]将其理论化为
式中,h为当地水深;x为距海岸线距离;m为统计得到的系数;A为剖面尺度参数,与当地泥沙粒径有关。依据Dolan等[38]对美国东海岸及墨西哥湾504条剖面观测结果,可知m值分布符合高斯分布,期望值为2/3,研究中一般可取2/3作为m的值[39]。
本文为进行海岸线精细提取,考虑研究区海滩地形和潮位变化,基于实测海滩剖面数据和Bruun-Dean平衡剖面模式,建立拟合剖面模型进行潮位校正(图4)。由于模型拟合需要卫星过境时刻潮位值及研究区多年平均大潮高潮位值,而验潮站提供的潮位数据为当日高低潮和逐时潮位预报结果,可使用经典潮汐插值模型对卫星过境时刻潮高进行推算:
图4 拟合剖面潮位校正模型原理Fig. 4 The principle of the tide level correction model of the fitted profile
式中,t为卫星过境时刻;H为对应潮高;Hh为当日高潮潮高;th为当日高潮时;Hl为当日低潮潮高;tl为当日低潮时。
为进行数据拟合,以平均大潮高潮线作为校正后海岸线位置,将其对应高程面作为基面,对RTK设置高程和验潮站潮位数据高程作统一调整。而后基于Bruun-Dean平衡剖面模式,使用研究区实测剖面数据进行模型拟合,公式为
式中,x为与海岸线水平距离;h为对应水深;a为剖面尺度参数;n为系数,参数a、n由数据拟合得到。通过卫星过境时刻潮高求得h的数值,代入方程即可得到瞬时水边线相对海岸线的潮位校正距离。
研究区域为海南省海口市西海岸,东起秀英港、西至天尾角,海岸线总长度约12 km(图5)。该区域海岸背景复杂,岸线类型以砂质岸线为主,剖面形态为上凹耗散型,附近建筑物密集,部分岸段在沙滩后方建有海堤。
图5 研究区地理位置Fig. 5 Geographical location of the study area
本文使用了1景WorldView-2卫星于2019年9月23日摄取的多光谱影像,其幅宽为16.4 km,覆盖整个西海岸区域。在进行边缘检测岸线提取之前,按照辐射校正-几何校正-正射校正-影像融合的技术路线对原始影像进行了预处理。辐射校正包含辐射定标和大气校正两个步骤,将卫星传感器测量值转换为地表真实反射率;之后结合地面控制点(Ground Control Point,GCP)对影像进行几何精校正,控制均方根误差(Root Mean Squared Error,RMSE)在0.5个像元以内并使用基于有理函数(Rational Polynomial Coefficient,RPC)模型的正射校正法[40],消除地球曲率、地形起伏引起的影像畸变;最后采用超分辨率贝叶斯法自动融合法[41]对多光谱影像进行数据融合。融合后的影像空间分辨率由1.85 m提升至0.5 m,在保留空间纹理信息的同时,水陆边缘愈为清晰明显。
本文所用的数据还包括与影像获取时间同一季度的野外实测数据(2019年11月),包含海滩剖面数据、海岸线数据等。测量所用仪器设备为思拓力(STONEX)S9Ⅱ专业型实时差分定位移动站(RTK),控制平面精度为 -8~8 mm均方根(Root Mean Square,RMS),高程精度为-15~15 mm RMS;此外,还搜集了位于研究区域内的秀英验潮站(站点潮高基准面在平均海平面下150 cm)2015-2019年潮位数据,用于建立拟合剖面模型,整体技术流程如图6所示。
图6 砂质岸线精细提取技术流程Fig. 6 Sandy shoreline extraction technology process
为评估SE算法岸线识别效果,本文使用3种常用的边缘检测算子法(Roberts算子法、Canny算子法、LoG算子法)与之进行对比。筛选出一处研究区内背景较为复杂的岸段图像单元,分别使用4种算法进行边缘检测并获取结果(图7)。为评估边缘检测质量效果,本文使用帧速率(单位:fps)、连续性(Cdr)和F值3种评价指标,对4种算法进行定量评价(表1)。其中,帧速率用于体现边缘检测算法的运行效率;连续性能够反映海陆边缘提取的完整程度;F值则为精确率与召回率的调和平均数,是定量评估边缘检测算法精度的常用指标[42]。
图7 边缘检测算法效果Fig. 7 Edge detection algorithm effect
依据表1结果,对4种边缘检测算法输出结果进行对比分析。LoG算子法输出结果F值最低且Cdr值最低,尤其在水陆边界处出现大量断点,无法有效定位水边线;Roberts算子法和Canny算子法F值和Cdr值相对较高,可用于水边线的识别,但在噪声水平较高的位置出现明显边缘断点和伪边缘现象;本文所用SE算法输出结果清晰、细腻,较少出现断点和伪边缘现象,其评价指标F值和Cdr值在4种算法中均为最高,说明SE算法在水陆边缘检测精度和连续性方面效果最好。在算法运行速率方面,由于使用随机森林模型,SE算法帧速率达到60 fps,是边缘检测算子法的4倍。由此可见,针对高分影像边缘检测,SE算法相比于传统边缘检测算子法更加精准高效,适用于高分影像砂质岸线的识别和精细提取。
表1 边缘检测算法质量评价结果Table 1 Quality evaluation results edge detection algorithm
4.4.1 瞬时水边线提取
将WorldView-2影像数据输入训练好的SE算法模型,按2.2节和2.3节所述方法步骤提取得到数字化水边线,结果如图8所示。显见,SE算法输出水边线结果具有高信噪比、高精度、强连续性的特点,且在算法运行过程中,不需要手动调节阈值参数,自动化程度高,具有很强的场景自适应性。
图8 水边线检测结果Fig. 8 Detection results of land and water boundary lines
4.4.2 潮位校正与海岸线提取
本文使用10条研究区RTK实测海滩剖面数据,依据3.3节所述方法步骤建立拟合剖面模型:h(x)=axn(图9a),同时建立拟合线性模型:h(x)=bx(图9b)进行对比评估。由拟合结果(图9)可知,拟合剖面模型参数a拟合结果范围为0.162 8~0.206 7,期望值为0.184 7,参数n拟合结果范围为0.653 4~0.715 7,期望值为0.682 5,与Bruun-Dean平衡剖面模式中m期望值的2/3相近;拟合线性模型参数b的拟合结果为0.083 9。
图9 潮位校正模型拟合结果对比Fig. 9 Comparison of fitting results of tide level correction models
计算两种模型拟合评价参数(表2),对比可知,拟合剖面模型的确定系数(R-square)为0.93,高于线性模型的0.78,说明前者的拟合优度更高。在误差平方 和(Sum of the Squared Errors, SSE)和 均 方 根 误 差(Root Mean Squared Error, RMSE)方面,拟合剖面模型结果值相较于线性模型都更小,说明使用拟合剖面模型进行潮位校正,预测精度更高。
表2 拟合模型评价参数Table 2 Fitting model evaluation parameter
使用秀英验潮站潮位数据,依据潮汐插值模型和基面高程计算水深值h。将卫星过境时刻潮位数据对应水深ht代入拟合剖面方程,求得瞬时水边线与平均大潮高潮线(海岸线)之间的垂直距离,并依此将矢量水边线分段沿其垂线向陆移动,最终提取得到真实的海岸线结果(图10)。
图10 研究区海岸线提取结果Fig. 10 Extraction results of coastline in the study area
4.4.3 拟合剖面模型分析与评估
使用拟合剖面模型进行潮位校正,相较于传统线性模型,海岸线结果精度获得提升;与DEM高程模型潮位校正相比,在保障精度的同时,该方法仅使用海滩剖面数据,不需要进行全地形测绘,减少了野外实测工作量。
然而,在进行大范围、长时间尺度的砂质岸线演变研究时,建立模型所需的实测海滩剖面数据量依然较大,需对上述方法进行适当简化。针对方程:h(x)=Axm,Moore[43]通过大量实验数据和实地测量数据检验了所谓的剖面尺度参数A与泥沙粒径D之间的相关性;而沙粒沉降遵循斯托克斯定律,泥沙粒径D与沉降速率之间存在方程关系;Dean[44]随后提出了基于沉降速率w下参数A与粒径D之间的转换关系,公式为
式中,A单位为m1/3;w单位为cm/s。为便于查询,Dean[39]给出了参数A与泥沙粒径D之间的对应简表。在实际应用中,可基于研究区泥沙粒径数据,直接建立拟合剖面模型进行潮位校正。
4.5.1 精度验证结果
为验证结果精度,以RTK实测岸线数据作为参考岸线,使用断面法对提取岸线进行定量分析。研究区岸线总长度约为12 km,以5 m为间隔布置了2 324条断面,并将各断面与参考岸线的2 324个交点作为参考点集(图11)。
图11 断面及参考点布置图Fig. 11 Transects and reference points layout drawing
计算提取岸线相对于参考岸线的偏移量,统计各断面对应结果(图12)。红色实线为相对参考岸线0偏移位置,数据位于该线以上表示提取岸线相对参考岸线向陆偏移,位于该线以下表示提取岸线相对参考岸线向海偏移。经统计分析,提取岸线平均偏移量为2.34 m,说明提取岸线精度较高;标准差为1.22 m,说明数据离散程度较小、稳定性较强。
4.5.2 结果评价与误差分析
根据参考点集,求得提取岸线均方根误差为2.23 m,由于参考岸线测量精度为毫米级(详见4.2节),说明提取岸线定位精度优于2.5 m。依据我国基本比例尺地形图制图标准,本文提取岸线整体满足1∶5 000比例尺海岸带制图需求[1]。但由于本文所用WorldView-2影像空间分辨率优于1 m,岸线定位精度超过了两个像元,分析误差可能来源于以下几个方面:
(1)波浪因素对水边线位置的影响。同天文潮汐一样,近岸波增水、波浪爬高等因素引起的近岸水位变化会引起瞬时水边线的位置平移[45]。由于缺乏卫星过境时刻波浪观测数据,本文未进行波浪校正。
(2)影像背景复杂程度对边缘检测算法的影响。依据图12中数据,提取岸线相对于参考岸线最大向海偏移距离为5.33 m,最大向陆偏移距离为5.16 m。根据断面编号定位海岸线位置,发现偏移距离较大的位置多为建筑物密集位置以及岬角、河口等复杂地理地貌单元,使得该位置岸线提取结果精度较低。
图12 提取岸线偏移量统计Fig. 12 Extract coastline offset statistics
(3)Bruun-Dean平衡剖面模式局限性影响。本文提出的拟合剖面模型潮位校正法相比于传统线性模型精度更高,所用Bruun-Dean模式简明但存在参数A物理意义不明以及在破波带以内适用性较差等局限性,与实际海滩剖面形态依然存在一定偏差。
(4)不确定性误差累积的影响。由于遥感数据从获取、处理、分析等一系列过程中都有不同类型和程度的不确定性引入,并在进一步分析中传播,尤其是对于遥感岸线提取验证等定量分析,传感器误差、影像校正误差、人员误差等不确定性累积会对数据结果精度产生更大的影响。
使用SE算法对高分辨率海岸影像进行了水边线检测与提取。通过与传统边缘检测算子法(Roberts算子法、Canny算子法、LoG算子法等)对比研究,发现SE算法具有更高的边缘定位精度、更强的噪声鲁棒性、更快的运行速率以及针对复杂背景海岸线识别的自适应能力。
针对线性模型误差较大的问题,提出了基于拟合剖面模型的潮位校正方法。通过RTK实测剖面数据结合经典的Bruun-Dean平衡剖面模式,建立拟合剖面模型,提升了遥感提取海岸线优化与校正的精度。
基于同期实测岸线数据,使用断面法对提取岸线进行了精度评估。定量分析了提取岸线相对参考岸线的偏移量,结果显示,提取岸线的定位精度优于2.5 m,满足1∶5 000比例尺海岸带制图需求。
下一步工作中,拟制作高分辨率遥感影像数据集专门用于算法的训练和测试,并研究学习适用性更强的平衡剖面模式以及近岸波浪岸线校正方法,进一步提升岸线提取的精度和效率。