渤海湾大规模围填海导致的岸线变化及潮流场响应分析

白玉川，史丰硕，徐海珏，黄哲

（1.天津大学 水利工程仿真与安全国家重点实验室，天津 300350；2.天津大学 河流海岸工程泥沙研究所，天津 300350）

渤海湾区域具有地形坡度较小、泥沙资源丰富、海岸线弧度较大等自然条件特点，有利于进行大规模的填海造陆活动。20 世纪70 年代以来，此区域在经济和政策的双重驱动下，对土地资源的需求更加强烈，进一步加快了填海造陆的进程（国家海洋局，2016；朱高儒等，2012）。渤海湾岸线人工化程度较高，截至2015 年，渤海湾沿岸除部分入海河口附近的海岸线为河口型自然岸线外，其他海域沿岸几乎均为人工岸线（孙晓宇等，2014），人类活动造成了岸线的永久改变，进而引起附近海域潮流场的变化（张立奎，2012）。潮流的流速、流向决定了研究海域盐度、污染物、排海高温水体等的扩散过程和泥沙运动沉积过程，从而影响该海域的生态环境和底床形态变化（Lu et al，2009）。因此，对于渤海湾岸线及潮流场变化的研究具有很重要的现实意义和工程意义。

对于渤海湾岸线及潮流场的研究主要有卫星遥感解译和数值模拟分析两种方法。卫星遥感是一种将地学分析、数学方法和物理手段综合在一起的应用技术（庄悦，2019），具有观测范围广、精度较高、获取信息量大、可反演等优点（侯庆志等，2013；Song et al，2004），与传统方法相比，在时空尺度上具有显而易见的优势，在岸线、土地属性等方面的时空分析中发挥着越来越重要的作用。众多学者利用此方法分析了渤海湾地区岸线的演变规律。李建国等（2010） 利用多源遥感数据提取了2000—2009 年天津滨海新区海岸线，得到岸线的变化特征并对驱动力进行了分析。朱高儒等（2012） 基于1974—2010 年MSS/TM/ETM 影像资料，利用RS/GIS 技术分析了渤海湾西北岸填海造陆的动态空间分布和数量变化特征。侯庆志等（2013） 通过分析1993—2011 年TM 数据与HJ-1A/1B 卫星影像，结合2003 年水下地形图，利用ArcGIS 分析了近20 年来曹妃甸海域围海工程造成的岸线及水动力环境变化。叶小敏等（2016）利用Landsat-5/TM、Landsat-7/ETM 及HJ-1A、B/CCD遥感数据分析了渤海湾近30 年来水域面积、岸线长度和岸线分形维数的变化，发现养殖场、港口及围填海工程建设是渤海湾岸线变迁的主要因素。

随着计算机和模拟技术的不断发展，数值模拟逐渐成为研究潮流场、泥沙运动及污染物迁移核扩散的重要手段，众多学者利用此方法在渤海湾区域取得了许多研究成果。王永刚等（2014） 基于ROMS 模式根据不同水深及岸线数据，建立了渤海海域潮波数值模型，提高了计算精度，研究表明1972—2002 年期间M2分潮无潮点向东北方向迁移了约30 km。宋军等（2016）建立了高分辨率数据同化模型系统，在ROMS 模式最优解的基础上进一步研究了黄渤海M2分潮的潮汐特征、潮余流、潮能及其扩散、潮混合和潮动能平衡。韩树宗等（2012）基于FVCOM 三维水动力模型针对近岸海域围填建设研究了天津近海海域流速、流向、冲淤环境、余流等方面的变化。蒋秋飙（2011）通过建立二维数学模型模拟工程区周围潮汐环流、海浪及泥沙运动规律，分析了码头区建设对海洋环境带来的具体影响，并为码头区各方面的科学决策提供了重要依据。陆永军等（2007）建立二维水沙数学模型分析了波流共同作用下曹妃甸港区的泥沙运动特征，发现甸头岬角效应是维持深槽水深的主要因素。夏波等（2006） 利用SWAN 和ADCIRC 模型创新性地建立了渤海湾西南岸风浪、潮汐、水流联合作用耦合数学模型，分析了近岸区水位、流场对风浪模拟结果的影响，发现水位变化对结果有明显影响。李希彬等（2018）利用FVCOM 模型模拟了未来大规模海岸工程对渤海潮动力环境的影响，并建议在工程建设之前，应充分考虑工程对水动力环境及水体交换能力的影响。聂红涛等（2008）应用数学模型模拟分析了海岸带不同开发活动对渤海湾近岸海域水环境的影响，研究结果表明，一次性排放大量污水会增大河口及近岸海域污染面积及污染程度。

虽然部分学者利用遥感影像技术与数值模型相结合的方法来分析围填海对渤海湾区域的影响，但大都是针对某一工程单独进行的，计算范围较小（陈静等，2016；Cui et al，2011），大范围区域的研究成果较少，而从整体角度可以同时研究区域内多个工程，方便在同一条件下相互对比分析，有利于宏观上把控演变规律。因此，本文利用ENVI 提取了渤海湾2000—2015 年卫星图像数据中典型年份岸线，分析得到渤海湾岸线时、空演变规律，在此基础上建立二维数学模型重点研究渤海湾大范围围填海工程前、后不同岸线及地形条件下潮流场的变化，并为后续的海洋环境、海底泥沙运动、污染物扩散等方面的研究提供理论基础。

1 渤海湾2000—2015 年岸线演变分析

1.1 区域概况

渤海湾地处我国大陆东部北端，渤海以西，三面环陆（宋媛 等，2021），位于117毅35忆E—118毅51忆E，38毅00忆N—39毅14忆N 之间，海域面积约为1.59 万km2，是典型的半封闭缓坡淤泥质海湾，包括隶属于天津的滨海新区、隶属于河北省的曹妃甸工业区和黄骅港以及老黄河口，东南部为莱州湾（图1）。

图1 研究区域示意图

1.2 数据来源与处理方法

为了研究渤海湾区域2000—2015 年岸线演变规律，本文使用美国陆地卫星Landsat 系列遥感数据，其分辨率为30 m，选取覆盖研究区域的卫星影片，轨道号为p121—r34 及p122—r33，条带号分别为122、033，利用每个典型年份（2000 年、2004 年、2011 年和2015 年）当年一整年的图像数据进行岸线的提取。

首先，对于6、5、4 或5、4、3 波段利用ENVI（邓书斌，2014）对所有卫星图像进行假彩色合成，用以突出海岸信息；接着对假彩色合成图像进行分类，提取陆地与水体的分界线，此分界线是提取年份内所有时期高潮位的平均位置，即为平均高潮线。对于淤泥质海岸，利用这种方法能得到较为清晰的岸线，而对于部分凸堤或挡沙堤等港口岸线，也在提取时算作有效岸线；最后用ArcGIS对ENVI 识别出来的陆-水分界线进行人工修整，得到最终岸线并进行对比分析。由于岸线的走向是曲折蜿蜒的，在进行人工修整的时候不一定严格按照网格来确定，有时会取对角线连线作为最终岸线，精确到个位数以后会导致岸线长度不是30 m的倍数。不同年份岸线具体信息见表1，不同年份岸线时空变化见图2。

表1 不同年份岸线具体信息

1.3 渤海湾岸线演变分析

分析表1 和图2 可以发现，渤海湾岸线长度整体呈现增长趋势，按增长强度可以分为三个阶段：2000—2004 年为岸线缓慢增长阶段；2004—2011年为岸线快速增长阶段；2011—2015 年为岸线增长减缓阶段。多期影像及岸线长度数据表明，整体演变呈现岸线向海洋扩张的趋势（Zhu et al，2012）。

图2 渤海湾岸线时空变化

2000—2004 年期间的围填海规模仍处于较低水平，年均填海速率为4.78 km2/a，渤海湾岸线共增长了46.45 km。填海工程主要集中在滩涂区域，但也有不少填海区延伸到近海。这一时期岸线及填海造地面积变化最大的区域主要位于天津滨海新区的北疆、南疆和临港工业区以及曹妃甸工业园区，这是因为天津滨海新区东疆港区的围埝工程和曹妃甸连岛公路在这一时期基本完工，这也为下一阶段的大规模开发奠定了基础。

2004 年之后天津滨海新区、曹妃甸工业园区以及黄骅港等重点建设区的岸线长度增速明显，岸线向海扩张的速度也明显加快。2004—2011 年期间年均围填海面积高达96.96 km2/a，岸线总长度增加了197.93 km。岸线变动及填海面积变化范围主要集中在曹妃甸工业区和天津港北港区；另外在曹妃甸新区的曹妃甸新城、南堡油田填海区，滨海新区的北疆电厂、中心渔港、南港工业区等均进行了一定规模的填海工程。填海位置开始向更深的海域推进，最远处达到了6 m 等深线。

随着国家可持续发展战略的提出和生态环境保护意识的增强，围填海工程速度明显减缓，2011—2015 年期间渤海湾陆域年均增长面积为6.645 km2/a，岸线总长度增加了58.88 km，变化范围集中在黄骅港区域，黄骅港区岸线向海推进，同时其东南方岸线也迅速向海推进，这主要是由滩涂围垦用作渔业养殖所致。

总之，岸线长度与围填海面积的不断增长与近年来防波堤、护堤等海岸工程以及围填海工程陆续的建设密不可分。

2 渤海湾数值模型

2.1 数学模型

近年围绕渤海湾进行的围填海工程已致使岸线整体向外海推进，工程区域逐渐由高滩转为低滩，近岸滩涂向海淤涨，地形发生显著变化，导致了海湾内潮流动力场的改变。本文基于2000 年和2015年渤海湾岸线和地形资料，利用Mike 21 建立了渤海湾潮流场动力模型（衣秀勇，2014），模拟了渤海湾潮流动力场，在此基础上分析了大规模围填海工程实施前后动力场变化的趋势和原因。

平面二维水流连续方程：

动量方程：

2.2 模型设置

模拟开边界条件指定为水位过程，陆地边界为固边界。建立渤海整体区域大模型（包括辽东湾、渤海湾、莱州湾和渤海部分海域），模型的开边界取小长山岛—威海连线，小模型为渤海湾数学模型。

大模型水深资料取自中国人民解放军海军司令部航海保证部出版的黄渤海海域海图，小模型水深资料取自渤海湾海域海图，并将各地水深统一到1985 国家高程基准面，不同时期的模型分别采用当年的水深资料来计算。渤海大模型开边界水位条件由ChinaTide 潮汐预报软件提供（李孟国等，2007），考虑Q1、P1、O1、K1、N2、M2、S2、K2、Sa 9 个主要分潮，调和常数由数值计算结果进行调和分析得到。模拟时间为2000 年2 月10 日00 颐00 颐00—2001 年3 月30 日00 颐00 颐00 和2015 年2 月10 日00 颐00 颐00—2016 年3 月30 日00 颐00 颐00，时间步长为30 s。采用三角形无结构化网格来划分，网格尺度为3000 m，2000 年模型单元数为33 187，结点数为16 991，2015 年模型单元数为32 864，结点数为16 917。渤海湾小模型开边界条件由渤海大模型插值计算得到。模拟时间为2000 年3 月10日00 颐00 颐00—2001 年3 月10 日00 颐00 颐00 和2015 年3 月10 日00 颐00 颐00—2016 年3 月10 日00 颐00 颐00，时间步长为30 s。采用三角形无结构化网格来划分，网格尺度为1000m，2000 年模型单元数为31 347，结点数为16 039，2015 年模型单元数为88 096，结点数为45 026。大、小模型的网格分布如图3 所示。

图3 模型网格分布图

渤海大模型和渤海湾小模型的底床摩擦力均选择曼宁系数，曼宁系数取值为0.0125 s/m1/3，该值的选取是基于王万战等（2007）对渤海流场基本特性的率定结果。

2.3 模型验证

为了验证潮流数学模型的准确性，引入模型技能得分（SkillScore）来量化模型误差（Song，2013a，2013b），SS 定义为均方根（RMS）误差与观测值标准差之间的比值（Murphy，1988），当SS 跃0.65时，认为模拟得很好，当0.65逸SS逸0.2 时，认为模拟得较好，当SS臆0.2 时，认为模拟得较差（Allen et al，2007）；同时引入相关性系数（Correlation Coefficient）来量化模型相关性，当CC 跃0.8 时属于强相关，当0.8逸CC逸0.5 时属于中等相关，当0.5逸CC逸0.3 时属于弱相关，当0.3 跃CC 时属于不相关。

分别利用环渤海各地验潮站的实测资料（国家海洋信息中心，2004）与计算结果进行比较，选取塘沽、黄骅港和东风港三个潮位观测站2000 年4月5—8 日4 天的观测值与计算潮位对比分析，结果如图4。

从图4 可以看出，各验潮站潮位观测值与模型计算值吻合良好，高、低潮时刻及潮位过程曲线趋势都基本一致，模型技能得分均在0.918 以上，相关性系数均在95.9%以上，因此可以认为计算结果是准确的，能够为渤海湾小模型提供精确的开边界水位条件。

图4 各验潮站潮位过程对比

为验证小模型的准确性，同样利用上文提到的模型技能得分（SS）来量化模型误差，利用相关性系数（CC）来量化模型相关性。选取天津港水文全潮观测站2000 年5 月18 日12 时至5 月19 日14时大潮观测值与小模型计算得到的流速、流向进行对比验证。天津港海区内8 个验潮测点的布设及位置见图5。将不同测点的实测资料同计算结果绘制在同一幅图像中进行对比，结果见图6。

图5 天津港区水文测站布置示意图

图6 实测与计算流速、流向对比

从实测和计算流速对比图中可以看出，1 号测点实测流速在研究时段内均大于计算流速，而其他测点的实测流速与计算流速曲线几乎重合，相差很小，这是由于1 号测站位于南疆港防波堤的西南侧，该处受到东疆港防波堤的阻力作用，潮流动力减弱，并且在模型中未考虑波浪、风等因素，因此，计算结果和实测数据存在一定的误差。流速整体技能得分为0.838、整体相关性系数为94.2%。8个潮流验证点的计算流向与实测流向相比仅有个别点存在微小差异，其余各点基本一致。流向整体技能得分为0.87、整体相关性系数为93.9%。总体而言，小模型计算准确性高，可以应用于后续计算。

3 结果分析

3.1 大规模围填海前后潮流动力场分析

由遥感影像岸线分析图（图2） 可以明确看出，2000—2015 年渤海湾岸线变化主要发生在曹妃甸工业区、天津港和黄骅港等海域，在此期间渤海湾海域面积减少了1178 km2，约占渤海湾海域总体面积的12%，其中围填海工程是造成渤海湾海域面积减少的主要因素。本文以天津港涨急落急时刻为基准，选取2000 年与2015 年两种地形条件下渤海湾大潮期间涨急时刻和落急时刻的水动力场进行分析。

由表2 可以看出，大规模围填海实施前后涨急最大流速和落急最大流速均呈现出先增大后减小的趋势，而涨急平均流速与落急平均流速变化都不大。为了更加直观地展示出渤海湾水文全潮流场的改变，分别绘制了渤海湾2000 年和2015 年大规模围填海前后涨落急流速差值等值线图进行对比，如图7 所示。

表2 2000—2015 年涨落急时刻流速对比单位：m·s-1

图7 大规模围填海前后涨落急流速差值等值线图

由图分析得到，2000 年渤海湾内最大涨落潮流速差值范围为-0.475 耀0.683 m/s，2015 年渤海湾内最大涨落潮流速差值范围为-0.492 耀0.717 m/s，即大规模围填海实施后流速差值较实施前有小幅增加，可见工程后研究区部分海域（工程区附近海域）潮流运动更加剧烈，流速变化更大。近岸海域和开边界处差值均有负值，渤海湾中部海域的流速差值多为正值，范围在0.1 耀0.15 m/s 之间。对比图6（a）和图6（b）可以发现，2000 年曹妃甸工业区、天津港和黄骅港附近海域流速差值等值线由近岸向外海方向分布均匀：曹妃甸工业区附近海域差值在0.1 耀0.15 m/s 之间，天津港和黄骅港附近海域差值在-0.1 耀0.15 m/s 之间。2015 年三个主要工程区附近海域的流速差值等值线由近岸向外海方向分布较为散乱：曹妃甸工业区附近海域流速差值增大，最大值达到0.3 m/s，天津港附近海域流速差值范围增大至-0.15 耀0.15 m/s，黄骅港附近海域最大差值同样有所增大，在港区岬角处流速差值达到了0.3 m/s。由此可见，渤海湾大规模围填海工程对三个主要工程区的流场影响都较大，这是由于防波堤对水流起着阻碍作用，导致流速减小，而在深水航道内和港区岬角处由于挑流的作用，流速明显增加，导致岬角处有被冲刷侵蚀的风险。

为了更直观地了解2000 年和2015 年渤海湾流速场变化，分别绘制了涨、落急时刻大规模围填海工程实施前后水文全潮流速差值图，如图8 所示。

图8 工程前（2000 年）、后（2015 年）差值流速场

从图8 可以看出，围填海工程实施后近岸海域流速有所减小，变化幅度约为0.15 m/s。与落急时刻不同的是黄骅港航道南部海域流速在涨急时刻有所增加，且涨急时刻流速变化幅度较大，工程区附近流速变化更为明显；离岸海域流速变化较小，变化幅度小于0.05 m/s，曹妃甸和滨州港之间海域、黄骅港附近海域流速变大，最大增幅能达到0.5 m/s以上。天津港附近海域流速变化复杂，近岸海域的潮流流向变化较大。

3.2 大规模围填海前后潮流调和分析

潮流通常是指天文潮汐涨、落而导致的海水流动。潮流流向一般呈现旋转式变化，一个潮周期内左旋或右旋360毅（徐辉奋等，2011），潮流性质依下式划分：

式中WK1、WO1、WM2分别是K1、O1和M2分潮潮流椭圆长半轴。如果F臆0.5，则为正规半日潮流；如果0.5 约F臆2.0，则为不正规半日潮流；如果F 跃2.0，则视情况分为不正规半日潮流或正规日潮流。本文在3 个主要工程区域共设置9 个特征点，如图9所示。利用小模型计算的2000 年和2015 年5 月1日00 颐00 颐00—6 月1 日00 颐00 颐00 潮流模拟结果，对深度平均流进行调和分析得到9 个特征点大部分分潮的特征值。本文主要针对M2分潮进行分析，并绘制M2分潮潮流椭圆。在表3 中给出9 个特征点的F 值，可以看出F 值在0.5耀2.0 之间，即渤海湾潮流属于不正规半日潮流性质（娄安刚 等，2002）。

表3 大规模围填海前后特征点的F 值

图9 渤海湾特征点位置示意图

潮流的运动形式可用潮流椭圆旋转率K 来判定，K 为潮流椭圆短轴与长轴的比值。当K 绝对值大于0.25 时，潮流表现为较强的旋转性；当K绝对值小于0.25 时以往复流为主，当K 值为正时表示潮流沿逆时针旋转，负值表示沿顺时针旋转（乔方利，2016），选取M2分潮来分析研究区域的潮流运动形式。

由表4 可知，对于M2分潮，黄骅港附近特征点g 在填海工程实施前潮流运动形式为旋转流，实施后转为往复流，这是由于在黄骅港附近海域修建大量防波堤等水工建筑物后，潮流运动受地形约束作用强烈，出现较明显的往复流特征。其他特征点所在位置潮流运动均为往复流，工程前研究区域潮流主要为逆时针旋转，但在工程后研究区域内潮流顺时针旋转范围增加，多发生在深水航道附近和受防波堤掩蔽的海域中。

表4 大规模围填海前后特征点的K 值

根据调和分析结果，绘制了M2分潮潮流椭圆。潮流椭圆的长轴方向为最大流速方向，短轴方向为最小流速方向（图10、图11） （Xu et al，2016）。在研究海域内围绕主要工程区域共提取了9 个特征点，以此分析工程建设区附近海域潮流的变化，见图12；根据前文涨急、落急分析，对于远离工程区的渤海湾区域，其流态受大规模围填海的影响较小，因此不再选取特征点进行分析。

图11 2015 年渤海湾整体M2 分潮椭圆

由图10 可以看出，大规模围海造陆工程前，渤海湾M2分潮潮流场分布均匀：由外海向近岸，M2分潮最小流速呈现出先增大后减小的趋势，最大流速则逐渐减小，在临近海岸时，最小流速和最大流速均有明显的减小。在渤海湾西北角附近海域M2分潮主要为西北向往复流；在曹妃甸附近海域潮流为旋转流，随着水深的增大，逐渐转为东西向的往复流；渤海湾中部大部分海域为东西向的往复流，在黄骅港岸线变动较大的区域附近转为旋转流，这是由于黄骅港区建设的挡沙堤凸出海岸，阻碍了水流的往复运动，产生了局部漩涡而导致的；渤海湾南部海域运动方式以旋转流为主，但是在岸线附近最大流速和最小流速急剧减小，主要以南北向往复流的方式运动。

图10 2000 年渤海湾整体M2 分潮椭圆

随着围海造陆工程的开展，渤海湾海域M2分潮运动方式也随之发生改变。分析图11 可以得到，大规模围填海工程实施后，由外海向近岸，渤海湾大部分海域M2分潮的最小流速有所增大，最大流速逐渐减小；在工程区附近的M2分潮运动由于受地形变化影响，流向发生变化，呈现出与深水航道轴线方向和防波堤堤线方向平行的往复流动。

为了更直观地进行工程前后渤海湾M2分潮潮流运动特征的对比，现将三个主要工程区附近9 个特征点不同年份的潮流椭圆绘制在同一幅图中，如图12 所示。

图12 不同年份工程区特征点潮流椭圆对比

通过对比分析可以发现，工程前后曹妃甸附近海域M2分潮的运动变化较小，较为明显的是最小流速的减小，工程后该海域主要以往复流为主，往复流运动方向主要由防波堤方向决定：最大流速方向与防波堤方向近似平行；天津港位于渤海湾西部，工程后该海域M2分潮的最小流速明显减小，致使往复流范围增大，旋转流范围减小；黄骅港海域在工程后潮流变化较大，主要表现在最小流速的减小和运动方向逐渐转向与深水航道轴线或防波堤堤线平行方向，往复流范围增大。围填海工程后，原有的自然岸线被码头、船坞、防波堤等人工岸线所代替，这不仅改变了原有的岸线属性，工程区所在海域的地形也随之发生变化，潮流运动受到地形的约束而发生改变：在深水航道内，M2分潮的最大流速有明显增大的趋势，最小流速减小，运动方式多为往复流；在港区岬角处，由于受到挑流的作用，最大流速和最小流速增大，往复流和旋转流同时存在，导致岬角附近存在被冲刷侵蚀的风险，降低了港区的整体稳定性；而受防波堤掩蔽的海域以及防波堤后方海域，M2分潮运动方式以往复流为主，受到建筑物的阻碍作用，最大、最小流速均有所减小，导致泥沙淤积的风险加剧。

4 结论

通过分析渤海湾不同时期的卫星遥感影像发现，渤海湾海岸线发生了巨大变化。2000—2015年海岸线不断地向外海推进，增长速度先快后慢，基本形成了以北部曹妃甸、西部天津港和南部黄骅-滨州港为中心的三大人工岸线聚集区。渤海湾北部向海推进最大距离约为20 024 m；渤海湾西部向海推进最大距离约为14 670 m；渤海湾南部向海推进最大距离约为22 011 m。

大规模围填海工程后，在岸线与地形变化较为明显的工程区域附近，潮流场变化较大，2000 年和2015 年渤海湾内最大涨落潮流速差值范围分别为-0.475耀0.683 m/s 与-0.492耀0.717 m/s，流速变化更大，潮流运动更剧烈；近岸海域流速有所减小，变化幅度约为0.15 m/s，但在工程区附近流速增大，增幅可达0.5 m/s 以上。

工程后渤海湾M2分潮流速变化明显：深水航道内最大流速增大、最小流速减小，潮流运动形式多为往复流，但在港区附近往复流和旋转流同时存在，黄骅港区由于挡沙堤的建设产生了局地涡旋，导致其附近潮流运动形式变为旋转流；岬角处由于受到挑流的作用，最大、最小流速都有所增大，导致岬角附近存在被冲刷侵蚀的风险，降低了港区的稳定性；受防波堤掩蔽的海区及防波堤后方海域，运动形式以往复流为主，最大、最小流速均有所减小，泥沙淤积的风险加剧。M2分潮潮流椭圆分布工程前、后也发生较为明显的变化：工程前以逆时针旋转为主，工程后潮流顺时针旋转范围增加，多发生在深水航道附近和受防波堤掩蔽的海域中；在工程建设区域，围海工程的修建改变了其最大流速运动方向，多为与建筑物布置方向、深水航道轴线方向一致，且往复流的范围增大。