1979—2012年莱州湾南岸海水入侵与区域海岸线变动时空耦合分析

蒙永辉， 王集宁， 张丽霞， 罗 梅

(山东省地质环境监测总站,济南 250014)

0 引言

海水入侵问题是近海岸带地区地质环境灾害研究领域中的一个热点，也是国际学术界一直关心的一个焦点环境问题[1-6]。海水入侵是指在人为和自然因素的干扰下，滨海地区海水和陆地地下淡水的平衡状态遭到破坏，导致海水不断向内陆方向侵入的过程和现象[7-9]。虽然海水的密度(矿化度)高于淡水，但是陆地淡水的水位高于海水，所以滨海地区的地下淡水和联通海水通过压力抗衡，在理论上能够维持一个相对平衡稳定的界面。但是，这种稳定一旦受到外界的干扰，例如： 气候变化、淡水资源的过渡开采、海平面上升等，平衡就会被打破，海水就会不断地侵入内陆淡水层。大量的海水向陆地倒灌，就会严重破坏区域生态环境，导致原本可以利用的淡水资源被污染，地下水水质变咸，土壤发生盐碱化，生产和生活用水困难，对区域社会发展造成诸多负面影响并引起一些环境问题，如淡水水资源匮乏、耕地资源退化、土壤次生盐碱化、原生植被凋亡等[10-14]。

位于山东省潍坊市北部的莱州湾南岸海岸带地区是我国最早发现海水入侵现象的地区之一[9,15]，也是我国目前海水入侵最为严重的地区[8,16]，已引起了当地政府部门和研究人员的高度重视。山东省人民政府早在1990年就出台了山东省莱州湾地区海水侵染灾情分析与综合治理方案，1994年出版了山东省莱州湾地区海水侵染综合治理规划[17]。相关领域科研工作者在海水入侵特征、成因、过程、模拟和评价方面开展了大量研究，取得了丰富的成果[3-5,18-27]。在海水入侵过程与特征研究方面，李新运等[20]对莱州湾南岸的海水入侵过程进行了分析和预测； 孙云华等[23]分析了1979―2008年海岸地貌过程演变与人类活动地下海水入侵的关系； 蒙永辉等[24]研究分析了潍坊北部海咸水入侵的特征和现状。在海水入侵成因与对策研究方面，刘桂仪[25]探讨了莱州湾南岸的海水入侵原因和防治对策； 丰爱平等[21]分析了莱州湾南岸1980―2000年海岸侵蚀过程并分析了海水入侵原因； 孙云华等[23]以1979、1989、2001和2008年4个时相的遥感影像为数据源，采用景观空间分析方法，分析了人类活动对海岸带地貌过程及海水入侵的影响。在海水入侵的模型模拟与评价研究方面，李福林[26]对莱州湾东岸海水入侵进行了动态监测与数值模拟研究； 苏乔等[27]对莱州湾南岸海水入侵现状进行评价； 陈广泉等[22]评价分析了莱州湾地区海水入侵灾害风险。综上所述，对莱州湾地区海水入侵研究的成果相当丰硕，但是对海水入侵特征过程仍缺乏定量系统的研究，对最新的海水入侵动态关注较少； 对海水入侵变化原因多为定性描述，缺少定量探讨； 对海水入侵年际时空变化规律并不明确[8]。

本文在已有莱州湾海水入侵监测资料的基础上，结合现代遥感观测手段和空间分析技术，对1979―2012年莱州湾南岸海水入侵发生与发展动态演变过程进行数字重建，根据相应时间段的海岸线信息，分析其时空演化特征，重点探讨和分析区域海水入侵过程与海岸线变迁的时空耦合关系。

1 研究区与数据源

1.1 研究区概况

本文的研究区位于莱州湾南岸滨海地区(图1)，地处胶东半岛西北部，濒临渤海，地理坐标在118°40′～119°40′E， 36°40′～37°20′N之间，涉及寿光、潍坊(寒亭)和昌邑3个市。气候属暖温带季风型大陆性气候，干旱少雨，年内降雨量不均。地质构造复杂、地貌类型多样，主要为由河流冲洪积与海水作用形成的冲海积平原，广泛分布着各类第四纪沉积物和粉沙质海岸地貌[9]。其中第四系沉积层厚度较大，平均200～300 m，沉积物成因类型以海相、湖沼相、河流冲积物为主； 入海河流形成的古河道带透水性较强，成为咸水入侵的良好通道[24]。该地区资源丰富、人口密集、海洋经济发达，拥有得天独厚的卤水矿资源，是我国沿海最大的卤水矿区，但也是我国最早发现海水入侵现象的地区之一[16]。受气候变化和人类活动的共同影响，该区域已经成为我国目前海水入侵最严重的地区，给当地工、 农业生产和人民生活造成了极大危害。

图1 研究区位置 (影像底图为TM B4(R),B3(G),B2(B)假彩色合成影像)Fig.1 Location of study area

1.2 数据源

本文中使用的海水入侵数据主要根据鲁勘字[2011]14号项目《黄河三角洲高效生态经济区(潍坊)海(咸)水入侵调查与监控预警系统建设报告》和莱水字[2011]39号文件，以及相关文献中的数据综合整理得到[17,21,23]。对海水入侵数据整理、利用中，存在空间和尺度的误差，需要进行空间配准和一致性检验，以去除彼此矛盾的数据。根据这些数据的采集时间，同时考虑遥感数据质量、受云层干扰情况等，去匹配相应时间的卫星影像。最终，本文挑选6期(1979, 1990,1995,2000,2003和2013年)遥感数据用于监测岸线变化(表1)。对遥感数据进行了预处理，将DN值转换为光谱反射率[28-29]。

表1 用于海岸线监测的遥感数据Tab.1 Remote sensing data used for shoreline monitoring

2 研究方法

本文的总体技术流程如图2所示。

图2 技术流程图Fig.2 Technical flow chart

在多时相海水入侵数据整理和相应时段遥感影像选择的基础上，首先对海水入侵专题图进行配准和数字化，整合1979―2012年6期海水入侵锋线； 对遥感影像进行预处理，提取1979―2013年的海岸线位置分布信息。然后，采用纵向剖面分析(垂直断面建模方法)[28]，每隔200 m采样，获取每个空间采样点的定量变化信息； 再根据端点速率法(end point ratio，EPR)分析模型，分别计算和统计分析入侵锋线、海岸线的时间分段变化速率及多年平均变化率; 最后，通过Pearson相关分析与显著性t检验，分析海岸线变化与海水入侵在时间和空间上的耦合关系。

2.1 海岸线提取与锋线数字化

通过遥感影像监测海岸线的变化，在国内外已有大量的算法研究和应用。本文主要利用面向对象分类技术，根据岸线在遥感影像中的空间分布特征和水陆光谱差异，依据归一化差值水体指数(normalized difference water index，NDWI)类间光谱特征差异最大原理，完成水陆分离； 然后通过人工互信息操作，实现海岸线的高精度识别，具体技术流程参考文献[28]。首先采用大津算法(OTSU)完成NDWI水陆分离； 然后依据空间关系特征判断实现陆地水域与海洋判别； 最后通过互信息人工操作对海岸线提取结果进行优化处理和微地形修正，得到研究区的各期海岸线分布图。互信息后处理操作主要包括： 河口处理需保留大型河口港湾特征和小河齐陆地基线； 泻湖需考虑其内环境，当泻湖与海洋有宽阔水域通道时，在泻湖内部采集海岸线。海水入侵锋线是在几何精确配准的基础上，对历史时期的海水入侵锋线进行手工数字化，得到各个历史时期海水入侵锋线的空间分布位置。

需要说明的是，本文中的海岸线数据没有考虑平均大潮位修正，不是真正意义的海岸线； 准确地说是利用遥感影像提取的瞬时水边线。理论上做潮位修正效果可能会更好，但因这个区域地表变化很大，SRTM高程数据在这个区域不是特别准确，潮位修正后的结果和遥感影像套合不上。所以，综合考虑到本文中分析数据是通过EPR 模型计算多个时期时间序列岸线的平均变化速率，尽管用水边线代替海岸线时潮位会对分析结果会产生一定的影响(不可避免)，但基本可以近似地反映区域海岸线的动态变化趋势(与其他文献中描述的趋势基本一致)，具有一定的可比性，对分析海水入侵锋线位置影响不大。

2.2 纵向断面建模方法

垂直断面建模方法是一种对数字岸线线性目标变迁速率计算传统的方法，其中以端点速率法 (EPR)计算模型最为简单常用[28]。EPR模型统计计算岸线的变化速率原理如图3所示。

图3 EPR局部垂直(纵向)断面建模计算示意图Fig.3 Sketch map of EPR model

图中P1，P2，P3，…,Pn基于基线的局部垂直采样点(n的大小由采样密度决定)；timei是获取时间；D1，D2，D3，…，Di是垂线(纵向剖面线)与变迁线的交点；distance(Di-Dj)是Di点与Dj点之间的距离。

首先建立垂直剖面线，本文在提取海水入侵锋线和岸线的基础上，通过缓冲区分析建立基线，设置采样间隔为200 m； 然后根据垂线与所有断面上的岸线交点距离，计算平均变化率。那么，对于2个历史时期的岸线变化速率，可直接用2个历史时期的变迁距离除以时间来计算，EPR1,2的数学表达式为

(1)

式中：D1，D2分别为time1，time2两个历史时期岸线的位置；distance(D1-D2) 为D1，D2两点之间的距离。

对于多个时期的一组时间序列岸线的平均变化速率，EPRave的数学表达式为

(2)

式中：timelatest为最新历史时间；timeoldest为最老历史时间；Dlatest,Doldest分别为timelatest,timeoldest两个历史时期岸线的位置。

3 结果分析

3.1 海水入侵锋线动态变化过程

1979―2012年莱州湾海岸带区域海水入侵锋线位置动态变化见图4。从咸水入侵锋线图可以清楚地看出区域海水入侵的动态变化过程。1979―1989年是区域海水入侵最快的阶段，随后区域海水锋线移动明显变缓。图5为海水入侵EPR时空动态分布图。综合分析海水入侵的动态过程和EPR时空动态分布图，可以从总体上将区域海水入侵分为3个阶段： ①1979―1989年的快速入侵阶段。此阶段是莱州湾南岸海水入侵发展最快的时期，每年入侵面积增加29.22 km2； 海咸水平均每年入侵速率为380 m/a，最大入侵距离在昌邑市青乡南部，入侵速率达到690 m/a； ②1989―2000年的慢速发展阶段。此阶段莱州湾南岸海水入侵的速度明显变缓，无论是入侵速度，还是入侵面积，都明显下降。海水入侵锋线基本上在1995年的位置附近波动； ③2000―2012 年的相对稳定与徘徊阶段。此阶段海咸水入侵动态锋线基本处于徘徊状态，并且局部出现了退缩(图5(e))。

图4 1979―2012年莱州湾南岸海水入侵锋线动态过程 (影像底图为TM B3(R)，B2(G)，B1(B)假彩色合成影像)Fig.4 Processes of seawater intrusion strike in south coast of Laizhou Bay from 1979 to 2012

(a) 1978—1989年海水入侵速率 (b) 1989—1995年海水入侵速率(c) 1995—2000年海水入侵速率

(d) 2000—2008年海水入侵速率(e) 2008—2012年海水入侵速率(f) 1979—2012年海水入侵速率

图5EPR空间动态分布图

Fig.5SpatialdynamicdistributionofEPR

3.2 海岸线动态变化遥感制图与分析

从解译的结果可以看出，莱州湾南部岸线的时空变化并不是一致的(图6)。

图6 1979―2013年区域海岸线变迁图 (影像底图为MSS B3(R),B2(G),B1(B)假彩色合成影像)Fig.6 Regional coastline transition from 1979 to 2013

从空间上看，区域海岸线总体上呈现向陆地蚀退的趋势(局部区域受到人类活动的干扰(如人工修建海上设施)除外)； 从时间上看，不同地点区域海岸线的蚀退速率也不一样，相同地点不同时段的蚀退速率也不同。

3.3 海水入侵与岸线变化耦合关系分析

1979―2012年莱州湾南岸的海水入侵与岸线变化速率见图7。从图7可以看出2条EPR曲线的变化特征和时空耦合关系。在空间分布格局上，海水入侵的强度和海岸线的蚀退强度相关性具有很强的一致性； 在时间变化速率上，海水的入侵速率变化曲线和海岸线变化的变化速率曲线形状相似。这一变化特征说明，莱州湾南岸地区海水入侵和海岸线蚀退在空间分布上和时间变化上均具有一定的相关性。

(图中Y轴表示基准岸线，正值表示向海方向，负值表示向陆地方向)图7 1979―2012年莱州湾南岸岸线与海水入侵锋线变化速率图Fig.7 Change rates of coastline and seawater intrusion front line of Laizhou Bay from 1979 to 2012

表2从定量的角度统计了区域岸线蚀退和海水入侵的变化特征以及二者之间的相关关系。表中海岸蚀退正值表示向海方向淤进，负值表示向陆地方向蚀退； 海水入侵负值表示向陆地方向侵入； ** 表示在0.01水平(双侧)显著相关。

表2 岸线蚀退与海水入侵变化量统计与耦合关系分析Tab.2 Relationship analysis between seawater intrusion and coastline changes (m/a)

1979―2012年莱州湾南岸海水入侵速率均值达到了 177.23 m/a； 最大值分布在昌邑市青乡南部区域，入侵速率达到 435.28 m/a； 最小值分布在潍坊北部寒亭,为 124.02 m/a。1979―2012年区域岸线的蚀退速率约为10.44 m/a； 蚀退的最大速率为124.02 m/a，空间位置上分布在昌邑市北部下营镇； 蚀退的最小速率分布在寒亭区以北的岸线段。此处由于人工建造海上游乐场和港口，区域岸线呈现向海洋快速淤进，最大速度达到191.14 m/a。海岸线EPR与入侵锋线EPR相关系数为0.407，显著性水平P值为0.00，通过99%的信度水平(双侧)检验。由此可见，区域海水入侵与岸线蚀退在时空上存在强耦合关系。

4 结论与讨论

4.1 结论

本文对1979―2012年莱州湾南岸海水入侵锋线演化的时空过程进行数字重建和区域海岸线的多期遥感监测； 通过EPR定量模型探索性地分析了二者之间的时空耦合关系机制。研究得到以下结论：

1)区域海水入侵经历了由快到慢的变化过程，1990年以后速率明显减缓，入侵锋线基本稳定在1995年锋线附近； 2008―2012年入侵锋线局部发生后退，目前保持相对稳定。这一研究结论与之前相关学者的研究结果基本一致[22]。

2)除了局部人工造陆导致岸线向海扩张外，1979―2012年区域海岸线以蚀退型海岸为主，但是蚀退速度在空间上差异明显。

3)海水入侵锋线变化与海岸线进退二者之间在时空上存在强耦合关系，相关系数达到0.407，显著性水平P<0.01(双侧)，说明海岸线工程对区域海水入侵速率有显著的影响。以上研究结论可为区域海水入侵的预防和治理提供数据支撑和科学依据。

4.2 讨论

1)在变化原因上，影响咸淡水界面迁移速率的因素可分为自然因素和人为因素，自然因素主要是地质条件、大气降水减少和风暴潮等； 人为因素主要是地下水开采及卤水的开采。特殊的地质条件(广泛分布的第四纪海相、湖相沉积以及河流冲积物，透水性强)为研究区内海咸水入侵提供了有利条件； 而气候变暖及海平面升高是该区海咸水入侵的大环境背景，人类活动则改变了海水入侵时空变化进程和入侵强度[9,24]。

2)在变化特征上，1979―1989年海咸水快速入侵阶段，主要是由于地下淡水资源的过渡开采，地下水位大幅度下降； 1989―2000年海咸水慢速发展阶段，由于当地逐渐认识到海咸水入侵的危害，对地下水开采进行了适当限制，同时近海卤水的开采量在逐年增大； 2000―2012 年海咸水相对稳定与后退阶段，主要是得益于海水入侵治理工程得实施，通过协调和控制地下水与卤水的开采，较好地控制了咸淡水界面。因此，按照研究区咸淡水界面变化情况，大量的研究分析共同表明： 地下水资源开发和卤水资源开采起了决定性作用，是区域海水入侵的直接原因[3,19, 22-23]。

3)海岸线的进退是影响区域海水入侵的又一叠加因素。海岸线的淤进，尤其是人工修筑的岸线，在一定程度上可以阻止海水入侵。可能的原因是人工修筑的岸线有利于阻挡风暴潮，堵塞海水入侵的部分通道。

以上研究结论建立在文中数据和方法的基础上。由于EPR模型方法计算线性目标的变化率是通过剖面上两点距离与时间的比值得到的，对于短期的加速、缓慢或者逆转变化不够敏感。尤其是在计算多期平均变化率时，短期的逆转性变化对整体变化趋势会产生一定的影响。此外，多源数据的空间精度误差问题、时间尺度问题不可回避。在空间尺度上，遥感数据源空间分辨率为30～60 m，在提取的海岸线精度上难免有些误差； 海水入侵历史数据资料的几何变形，导致局部精度难以控制。在时间尺度上，受到观测数据的限制，通过6期数据分析变化过程，时间分辨率精度不够高，对演化进程分析尤其是时间拐点的确定不够准确。这些可能的误差会对数据分析以及部分研究结论产生一定的影响。在未来的研究中有些方面将有待于进一步改进和提高，例如加大海水入侵锋线动态监测时间密度，改进算法提高对短期动态变化的敏感性等。