基于计算机视觉的平水桥海滩风暴潮响应分析

摘要

为探究海滩对风暴潮的响应过程，采用视频监测设备和计算机视觉方法，通过设立在岸边高处的Argus相机获取海岸带彩色图像，利用改进后的Canny算法结合双阈值处理法自动提取水边线，使用相机标定获得水边线地理坐标，对多时相的水边线进行高程插值，得到海岸线位置及海滩地形。通过对秦皇岛平水桥海滩的连续监测分析，发现平水桥海滩对风暴潮过程短时响应剧烈，泥沙向西、向海输移，滩肩明显冲刷，砂质粗化，潮间带整体淤积，海滩坡度变缓，岸线向海推进4.34 m；风暴潮后的恢复期，平水桥海滩岬角掩护岸段高潮带微冲而中低潮带淤积，泥沙离岸输移； 平直岸段演变主要受中小恢复浪作用，呈现上淤下冲的态势，泥沙由中低潮带向高潮带输移。

关键词

计算机视觉; 视频监测; 水动力模型; 风暴潮

风暴潮是由强风或气压骤变引起的海平面突变现象，其在短期内向海滩输入大量的能量，这种极端天气条件下产生的强风-浪-流相互作用对海滩地貌的塑造比正常天气高出1~2个量级^［1-2］，进而造成海滩地貌发生显著变化^［3］。风暴潮过程常伴随有风暴增水，这往往会使风暴潮对海滩作用的上界向陆推移并造成沿岸侵蚀灾害；另外，波高增大会造成碎波带显著变宽，可能导致风暴潮巨浪的能量在到达海滩之前已经耗散殆尽，这种情况不造成强烈的海滩侵蚀^［4-5］。风暴潮作为海滩地貌短期演变的重要动力形式，一直是研究的热点和重点。对海滩进行监测是评估风暴潮影响的基本手段之一，监测方法包括人工实地测量和摄影测量^［6-7］，其中摄影测量可以有效解决风暴潮期间人工测量困难的问题。

随着计算机视觉技术的发展，摄影测量的精度和自动化程度都得到了极大提高。计算机视觉技术为海岸线的精确提取指明了方向。近些年来，计算机视觉技术在海岸工程中的应用逐渐增多。郭俊丽等^［8］基于Argus记录的舟山东沙海滩图像，研究了东沙海滩对两次连续风暴的响应，结果表明两次风暴使得海滩滨线较风暴前后退10.91 m。Cao等^［9］将Argus系统监测与传统监测相结合，利用计算机视觉技术对日照海岸的监测图片进行分析，研究结果表明养滩工程两个月后海拔-6 m以上的沉积物体积大致保持稳定，通过视频图像可以量化波浪的翻越特征。

研究海滩对风暴潮响应的方法通常为视频图像数据或数值模拟，本文将视频图像技术与数值模拟相结合，以提高风暴潮反演精度、更多元化地丰富海滩动力地貌过程研究。首先利用Delft 3D模型对2016年7月20日的“720”风暴潮过程的动力环境进行复演，然后基于计算机视觉技术处理“720”风暴潮前后的视频图像数据，反演风暴潮前后地形，结合模拟的动力环境与反演的地形，分析极端天气过程影响下秦皇岛平水桥海滩演变短时响应特征及恢复期过程。

1 数据来源与方法

平水桥海滩位于秦皇岛北戴河西海滩东部，为砂质海滩，岸线呈SW~NE走向，海岸沉积以中细沙为主，东侧分布有老虎石天然岬角和人工潜堤岬头，西侧为平直岸段。平水桥属于弱潮海岸，多年平均潮差约为0.74 m，最大潮差1.50 m。波浪以风浪为主，常浪向为东南向。2016年7月19~21日，受温带气旋影响，渤海沿海出现一次较强的风暴潮过程，最大风暴增水发生在河北省黄骅站，为115 cm；秦皇岛站达到了黄色警戒水位，风暴增水为59 cm ^［10］。

2014年10月，河北省地矿局第八地质大队在秦皇岛平水桥海滩安装了Argus视频监测系统，其完整地记录下了“720”风暴潮的全过程。Argus系统由现场视频测站、可编程控制器、现场高程控制点系统、专用通讯网、中心服务器系统等硬件和Deltares公司开发的后处理软件2个部分组成^［11］。如图1所示，图1a为Argus视频测站安装位置，图1b为摄像终端，图1c为中心服务器系统，图1d为6个摄像头监测图像融合后的海滩全貌。Argus系统对海滩进行24 h连续监测，图像采集频率为30 min，记录2560×1920 pixel²的瞬时图像和10 min时均图像。

图1  研究区域现场及Argus系统

Fig.1  Study area and Argus system

基于Argus监测系统采集图像可以反演得到地形等信息，其工作技术路线如图2所示。首先对所采集的海滩图像进行处理，结合潮流模型提供的潮位信息；再根据摄影测量原理进行坐标转换，即可得到水边线坐标；最后由不同时刻的水边线反演地形的变化。

图2  地形反演流程

Fig.2  Flow chart of topography extraction

1.1　数学模型及参数设置

基于图2的技术流程，需要建立精确的潮流模型，以提供拍摄时刻的潮位信息。采用Delft 3D-Flow模块建立二维潮流模型模拟研究区域的潮位过程，使用渤海-秦皇岛双重嵌套网格^［12］，通过潮位过程来控制模型边界，其中渤海大模型以大连潮位站至烟台潮位站的连线为边界，对大连至烟台沿线各点的潮位数据插值，得到大模型的海域开边界；秦皇岛小模型的边界由三条海域开边界和一条岸线闭边界组成，其中三条海域开边界的潮位条件由渤海大模型的计算结果给出，岸线闭边界流速为零。渤海大模型的曼宁系数取0.012~0.018，秦皇岛小模型的曼宁系数取均值0.014 5，计算时间步长均设置为1 min，水平涡黏系数取15 m²·s^-1。

通过布设在平水桥南面700 m处的波浪监测浮标提供每小时的波浪数据，浮标处水深约6 m。在“720”风暴潮期间，浮标遭到破坏，造成部分监测数据的缺失，因而采用第三代SWAN模型模拟研究区域波浪场，以补充缺失时间段内的波浪数据。其中风应力计算如下：

$$\overrightarrow{{\tau }_{\mathrm{a}}}={\rho }_{\mathrm{a}}{C}_{\mathrm{d}}\overrightarrow{W_{\mathrm{10}}}\left|\overrightarrow{W_{\mathrm{10}}}\right|$$

（1）

式中：ρ_a为空气的密度；W₁₀表示海洋表面上方10 m处的风速；C_d代表风拖曳力系数，是精确模拟风暴潮的关键参数，宋竑霖^［13］基于前人建立的风速与风拖曳力系数之间的关系，归纳得

$$\overrightarrow{C_{\mathrm{d}}}=(\mathrm{0.8}+\mathrm{0.123}{W}_{\mathrm{10}})\times {\mathrm{10}}^{-\mathrm{3}}$$

（2）

采用欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts， ECMWF）提供的风场数据作为模型中的风作用过程。图3为模型输入的2016年7月风速和风向。可以看出，风暴潮整个生长衰落过程自7月19日18时~7月21日7时，历时约1.5 d。风速从7月19日18时起开始逐渐增大，7月20日9时~7月21日5时风速均在10 m·s^-1以上，风向为80°~180°，其后至7月21日7时风速逐渐衰减。

图3  2016年7月ECMWF风速及风向

Fig.3  Wind speed and direction in July 2016 from ECMWF

1.2　数字高程模型

结合已建立的潮流模型，对摄得的图像经过预处理、水边线提取、坐标转换等步骤，得到拍摄时刻的水边线实地三维坐标，对多时相的水边线进行插值，得到研究区域的数字高程模型（Digital Elevation Model， DEM）。

1.2.1　预处理图像

现场视频测站测得的图片包含大量的噪声，对海滩进行水边线提取时，需先通过空间滤波图像处理方法对噪声进行预处理。空间滤波用于模糊处理和降低噪声。典型的随机噪声由灰度级的急剧变化产生，通过使用滤波器模板将图像每个像素的灰度值设置为其邻域内像素的平均值，从而降低噪声灰度的急剧变化，有效防止其在边缘检测识别中被误检。高斯滤波以二维高斯分布为基础生成滤波模板。滤波分布函数为

$$G(x,y)=\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }{\sigma }^{\mathrm{2}}}{e}^{-\frac{x^{\mathrm{2}}+y^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}$$

（3）

式中：x、y为像素点的坐标；σ为分布函数标准差；

G（x，y） 为高斯函数。

用高斯函数和输入图像的卷积形成一幅平滑后的图像如下：

$$f_{\mathrm{s}}(x,y)=G(x,y)\times f(x,y)$$

（4）

式中：f（x， y） 为输入图像；f_s（x， y） 为输出图像。

图4为瞬时图像与高斯滤波后图像的对比。经高斯滤波降噪后，图像噪点被平滑，更利于后续水边线的提取。

图4  高斯滤波效果

Fig.4  Effect of Gaussian filter

1.2.2　提取水边线

经空间滤波图像预处理后，对图像中的水边线进行提取，主要步骤为梯度计算、非最大抑制、双阈值处理、像素坐标导出。

计算图像中每个像素点的梯度幅度M（x，y）和方向a（x，y），则

$$M(x,y)=\sqrt[]{g_x^{\mathrm{2}}+g_y^{\mathrm{2}}}$$

（5）

a （x， y） = arctan（gx/gy）

（6）

式中：$g_{\mathrm{x}}=\frac{\partial f_{\mathrm{s}}}{\partial x},g_{\mathrm{y}}=\frac{\partial f_{\mathrm{s}}}{\partial y}$。

首先，定义0°、45°、90°和135°四个方向，使用梯度幅度和方向来估计每一点处的边缘强度和方向，初步检测出边缘像素。然后，对边缘像素应用非极大值抑制，以去除边缘检测过程中出现的杂散边缘信息。对每个像素进行非极大值抑制的算法是：①将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较；②如果当前像素的梯度强度与另外两个像素相比最大，则该像素点保留为边缘点，否则该像素点将被抑制。最后，使用双阈值处理和连接分析，从边缘像素中锁定待提取的水边线，并使用提取工具导出水边线边缘的像素坐标，提取效果见图5。

图5  水边线提取

Fig.5  Shoreline extraction

1.2.3　转换坐标

图像拍摄过程中可能会产生畸变，通过相机标定可消除畸变^［14］，并建立起图像坐标与真实坐标的对应关系。根据解析摄影测量中的经典共线方程，可得到像点坐标与海滩坐标的对应关系如下：

$$x=-f\frac{a_{\mathrm{1}}(X-X_{\mathrm{S}})+b_{\mathrm{1}}(Y-Y_{\mathrm{S}})+c_{\mathrm{1}}(Z-Z_{\mathrm{S}})}{a_{\mathrm{3}}(X-X_{\mathrm{S}})+b_{\mathrm{3}}(Y-Y_{\mathrm{S}})+c_{\mathrm{3}}(Z-Z_{\mathrm{S}})}$$

$$y=-f\frac{a_{\mathrm{2}}(X-X_{\mathrm{S}})+b_{\mathrm{2}}(Y-Y_{\mathrm{S}})+c_{\mathrm{2}}(Z-Z_{\mathrm{S}})}{a_{\mathrm{3}}(X-X_{\mathrm{S}})+b_{\mathrm{3}}(Y-Y_{\mathrm{S}})+c_{\mathrm{3}}(Z-Z_{\mathrm{S}})}$$

（7）

式中：（x， y）为像点的像素坐标；（X， Y， Z）像点对应的地面点的实地坐标；（X_S， Y_S， Z_S）为拍摄终端的实地坐标；f为拍摄终端主距；a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为外方位元素。

使用6个高程控制点的实测数据，采用最小二乘法平差计算求解得到外方位元素，完成相机标定。水边线的Z坐标即为图像拍摄时刻的潮位，将模型计算的潮位数据结合水边线像素坐标代入已标定的共线方程中，即可输出水边线的（X， Y， Z）坐标。

1.2.4　反演地形

将图片提取出的多个时刻的水边线视作等高线^［15］，然后使用克里金插值算法对多时相的水边线进行插值，即可得到海滩潮间带的DEM。

2 模型验证

采用实测资料对数学模型进行了验证，潮位验证资料为2016年7月19~23日秦皇岛海洋站实测数据^［16］，波浪验证资料为2016年7月19~21日平水桥波浪监测浮标实测波浪数据，地形验证资料为RTK-DGPS测量的2016年7月19日地形数据，其中波浪和地形资料由河北省地矿局第八地质大队提供。经验证，潮位和波浪的模型计算结果在大小和相位上与实测值均比较吻合，如图6所示。7月20日14：00~7月21日3：00期间增水均在0.4m以上，风暴潮增水峰值达0.55m，与实测值0.59m ^［10］相差0.04m，相对误差为6.78%；波高误差均值为0.15 m，波向相对误差为20.93%^［17］。

图6  潮位及波浪验证

Fig.6  Verification of tide level and wave

在潮间带内测量了36个高程控制点以验证反演地形的精度，部分验证结果如图7所示，高程绝对误差集中在0~0.30 m，平均值为0.12 m，与其他已有研究的误差相当^［8， 18］。反演基于同一套转坐标参数，在评价潮间带地形变化时候，部分误差会相互抵消，因此将该地形反演结果用于评价海滩潮间带的地貌演变时，具有较好的可靠性。

图7  反演地形的误差

Fig.7  Deviation of topography extraction

3 结果与讨论

3.1　动力变化

图8为Argus系统3号摄像机记录下的2016年7月20日海滩图像，此区域为海滩的平直岸段，且坡度较缓，水位增加对岸滩影响较大。可以看出，在风暴潮的袭击下，3号摄像机监视的区域显示出不同的状态。13：30水位明显较平时高；17：00海水漫上岸滩；18：30海水进一步蔓延，对附近房屋造成破坏；19：00，水位继续增高。

图8  Argus3号摄像头下的2016年“720”风暴潮过程

Fig.8  Process of 2016 "720" storm surge captured by the Argus camera No.3

图9为风暴潮期间波浪和潮位过程，有效波高从7月19日18时的0.63 m开始增加，至20日17时达到峰值2.22 m，21日7时衰减到0.70 m，整个风暴潮过程的平均有效波高为1.47 m。岸线与正北方向顺时针成65°~80°角，故波向在65°~280°之间的波浪均可以直接入射海滩。由图9a可知，风暴潮期间的波浪均可以直接入射海滩，风暴潮对海滩存在持续且较强的影响。由图9b可知，风暴潮发生时平水桥海滩正处于大潮潮情，最大潮差为1.45 m，高于全年平均潮差0.74 m，最大潮位出现在20日21时，与有效波高峰值出现时刻接近，增强了对海滩的影响，最大风暴增水达0.55 m。

图9  2016年“720”风暴潮期间波况、潮位及风暴增水

Fig.9  Wave, tidal level and storm surge elevation during the "720" in 2016

3.2　海滩地形演变

3.2.1　风暴期海滩冲淤演变

图10a、图10b分别为视频反演的2016年7月19日、7月22日的潮间带数字高程模型DEM，风暴潮对潮间带地形产生了重大影响，促使潮间带坡度趋于平缓，泥沙自滩肩向水下输移。进一步地，将两天的地形相减，得到风暴潮前后海滩地形冲淤变化，其中正值为淤积，负值为冲刷，如图10c所示。风暴潮期间，海滩沉积物整体呈离岸输运的趋势。海岸线（平均大潮高潮线）高程以上，海滩处于冲刷状态，潮间带处于淤积状态。大风浪导致岸滩及高潮带泥沙掀扬侵蚀，而在中、低潮带区域落淤，且自岸向海方向淤积厚度增加，低潮带最大淤积厚度达0.6 m。

图10  2016年7月19、22日潮间带DEM及7月19~22日地形变化

Fig.10  Topography of intertidal zone and its difference between July 19 and 22, 2016

潮间带冲淤分布存在明显的空间差异。海滩中部受风暴潮影响最为剧烈，其坡度比两侧坡度更为陡峭，且无岬角等地形掩护，根据地形动力学的观点，更陡的坡度对波能耗散性更弱，海滩受到风暴潮的影响也更为强烈^［19］。海滩西侧处于平直岸段，地形变化最为剧烈；海岸线以上出现强烈的侵蚀，最大冲刷深度达0.6 m，潮间带淤积明显。风暴潮期间波浪破碎，波生流方向自东向西，主导泥沙向西输移。平水桥海滩东侧，由于天然岬角及NE-SW向的人工潜堤岬头的遮蔽，海滩对风暴潮响应较弱，潮间带地形变化不明显，但冲淤变化情况较为复杂，整体以淤积为主，自东向西淤积厚度增加，靠近岬角侧存在部分区域侵蚀，天然岬角和人工潜堤岬头对岸段起到一定保护作用，提高了抵御风暴能力。

3.2.2　风暴期岸线变化

“720”风暴潮在短期内对海滩输入了大量能量，造成海滩地形的剧烈变化。图11为平水桥海滩2016年7月18日、22日的同一潮高下4、5号摄像头的实景照片。通过对比风暴潮前后的实景照片，发现滩肩上存在较为明显的冲刷，砂质粗化，如方框区内的滩面沉积物明显变粗；在同一潮高下的海滩出露面积明显减小，表明潮间带淤积，印证了风暴潮造成的滩肩及高潮带侵蚀，海滩整体坡度趋于平缓，岸滩及高潮带泥沙向海输移。对图10海滩2016年7月19日和22日地形进行处理，得到对应的海岸线位置，如图12所示。风暴潮发生前（7月19日），滩肩宽度为130.80 m，风暴潮之后（7月22日）为135.14 m，风暴潮导致滩肩冲刷，滩肩高程降低，冲刷下的泥沙落淤在潮间带，导致海岸线向海推进4.34m。

图11  2016年7月18日和7月22日同一潮高下海滩图像

Fig.11  Images of beach at the same tidal level on July 18 and July 22, 2016

图12  2016年7月19日和7月22日高潮线位置

Fig.12  Location of high tide line on July 19 and July 22, 2016

3.2.3　恢复期海滩演变

风暴潮过后的一段时间里，海滩处于恢复阶段。为更详细地了解恢复期的海滩地形演变规律，对风暴潮后14d（8月3日）的20张Argus图片进行解译，提取水边线，进而插值反演得到风暴过后8月3日的地形，如图13a所示；将8月3日与7月22日的地形相减，得到7月22日~8月3日地形变化，即风暴潮过后14d的海滩地形冲淤变化，如图13b所示。

图13  潮间带地形及地形变化

Fig.13  Topography of the intertidal zone and its change

由图13可知，自7月22日~8月3日，潮间带整体以淤积状态为主，在空间上将海滩大致分为东、中、西三个区域，海滩东部受人工潜堤岬头和天然岬角影响，地形变化较为复杂，高潮带轻微冲刷，中低潮带淤积，淤积厚度小于0.2 m；海滩中部为平直岸段，整体呈上淤下冲的状态，高潮带区域淤积厚度小于0.1 m，中低潮带冲刷，冲刷深度小于0.2 m；海滩西侧也为平直岸段，地形演变规律与海滩中部较为一致，即上淤下冲，但变化程度更为剧烈，高潮带最大淤积厚度达0.3 m，低潮带最大冲刷深度达0.5 m。

综观海滩整体，发现在7月22日~8月3日的潮间带演变中，平直岸段与岬角掩护岸段变化呈现出不同的规律，岬角掩护岸段高潮带微冲而中低潮带淤积，泥沙离岸输移；平直岸段演变主要受中小恢复浪的作用，呈现出上淤下冲的态势，泥沙由中低潮带向高潮带输移。

4 结论

本文首先建立了“720”风暴潮期间的水动力数学模型，复演了风暴潮过程的动力环境，然后结合视频监测，运用计算机视觉技术，反演出平水桥海滩2016年“720”风暴潮发生前后的潮间带地形，结合动力数据分析风暴潮对海滩的影响，得出以下结论：

（1）“720”风暴潮发生时平水桥海滩正处于大潮潮情，整个过程自7月19日18时起~7月21日7时止，共历时约1.5 d，风暴潮期间平均有效波高为1.47 m，最大有效波高为2.22 m，波浪均为入射海滩，最大风暴增水0.55 m。

（2）平水桥海滩对风暴潮过程短时响应剧烈。在强浪的影响下，泥沙向西、向海输移，造成滩肩明显冲刷，砂质粗化，潮间带整体淤积，从而导致海滩坡度变缓，岸线向海推进4.34 m。受波向、地形等因素影响，海滩对风暴潮的响应存在明显的空间差异，其中以平直岸段的海滩中部对风暴潮响应最为剧烈，海滩东部受人工潜堤岬头和天然岬角影响，地形变化幅度最小。因此可认为人工潜堤岬头对其影响岸段起到了一定的保护作用，提高了岸段的抗风暴潮能力。

（3）风暴潮后的恢复期，平水桥海滩平直岸段与岬角影响岸段呈现出不同的变化规律。岬角掩护岸段高潮带微冲而中低潮带淤积，泥沙离岸输移； 平直岸段演变主要受中小恢复浪的作用，呈现出上淤下冲的态势，泥沙由中低潮带向高潮带输移。

作者贡献声明

匡翠萍：论文撰写及修改。

刘 旭：数据处理、论文撰写。

夏子龙：数据处理、图像绘制。

朱 磊：现场工作。

丛 新：论文修改。