1979—2012年萊州灣南岸海水入侵與區(qū)域海岸線變動時空耦合分析

蒙永輝， 王集寧， 張麗霞， 羅 梅

(山東省地質環(huán)境監(jiān)測總站,濟南 250014)

0 引言

海水入侵問題是近海岸帶地區(qū)地質環(huán)境災害研究領域中的一個熱點，也是國際學術界一直關心的一個焦點環(huán)境問題[1-6]。海水入侵是指在人為和自然因素的干擾下，濱海地區(qū)海水和陸地地下淡水的平衡狀態(tài)遭到破壞，導致海水不斷向內陸方向侵入的過程和現(xiàn)象[7-9]。雖然海水的密度(礦化度)高于淡水，但是陸地淡水的水位高于海水，所以濱海地區(qū)的地下淡水和聯(lián)通海水通過壓力抗衡，在理論上能夠維持一個相對平衡穩(wěn)定的界面。但是，這種穩(wěn)定一旦受到外界的干擾，例如： 氣候變化、淡水資源的過渡開采、海平面上升等，平衡就會被打破，海水就會不斷地侵入內陸淡水層。大量的海水向陸地倒灌，就會嚴重破壞區(qū)域生態(tài)環(huán)境，導致原本可以利用的淡水資源被污染，地下水水質變咸，土壤發(fā)生鹽堿化，生產和生活用水困難，對區(qū)域社會發(fā)展造成諸多負面影響并引起一些環(huán)境問題，如淡水水資源匱乏、耕地資源退化、土壤次生鹽堿化、原生植被凋亡等[10-14]。

位于山東省濰坊市北部的萊州灣南岸海岸帶地區(qū)是我國最早發(fā)現(xiàn)海水入侵現(xiàn)象的地區(qū)之一[9,15]，也是我國目前海水入侵最為嚴重的地區(qū)[8,16]，已引起了當?shù)卣块T和研究人員的高度重視。山東省人民政府早在1990年就出臺了山東省萊州灣地區(qū)海水侵染災情分析與綜合治理方案，1994年出版了山東省萊州灣地區(qū)海水侵染綜合治理規(guī)劃[17]。相關領域科研工作者在海水入侵特征、成因、過程、模擬和評價方面開展了大量研究，取得了豐富的成果[3-5,18-27]。在海水入侵過程與特征研究方面，李新運等[20]對萊州灣南岸的海水入侵過程進行了分析和預測； 孫云華等[23]分析了1979―2008年海岸地貌過程演變與人類活動地下海水入侵的關系； 蒙永輝等[24]研究分析了濰坊北部海咸水入侵的特征和現(xiàn)狀。在海水入侵成因與對策研究方面，劉桂儀[25]探討了萊州灣南岸的海水入侵原因和防治對策； 豐愛平等[21]分析了萊州灣南岸1980―2000年海岸侵蝕過程并分析了海水入侵原因； 孫云華等[23]以1979、1989、2001和2008年4個時相的遙感影像為數(shù)據(jù)源，采用景觀空間分析方法，分析了人類活動對海岸帶地貌過程及海水入侵的影響。在海水入侵的模型模擬與評價研究方面，李福林[26]對萊州灣東岸海水入侵進行了動態(tài)監(jiān)測與數(shù)值模擬研究； 蘇喬等[27]對萊州灣南岸海水入侵現(xiàn)狀進行評價； 陳廣泉等[22]評價分析了萊州灣地區(qū)海水入侵災害風險。綜上所述，對萊州灣地區(qū)海水入侵研究的成果相當豐碩，但是對海水入侵特征過程仍缺乏定量系統(tǒng)的研究，對最新的海水入侵動態(tài)關注較少； 對海水入侵變化原因多為定性描述，缺少定量探討； 對海水入侵年際時空變化規(guī)律并不明確[8]。

本文在已有萊州灣海水入侵監(jiān)測資料的基礎上，結合現(xiàn)代遙感觀測手段和空間分析技術，對1979―2012年萊州灣南岸海水入侵發(fā)生與發(fā)展動態(tài)演變過程進行數(shù)字重建，根據(jù)相應時間段的海岸線信息，分析其時空演化特征，重點探討和分析區(qū)域海水入侵過程與海岸線變遷的時空耦合關系。

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

本文的研究區(qū)位于萊州灣南岸濱海地區(qū)(圖1)，地處膠東半島西北部，瀕臨渤海，地理坐標在118°40′～119°40′E， 36°40′～37°20′N之間，涉及壽光、濰坊(寒亭)和昌邑3個市。氣候屬暖溫帶季風型大陸性氣候，干旱少雨，年內降雨量不均。地質構造復雜、地貌類型多樣，主要為由河流沖洪積與海水作用形成的沖海積平原，廣泛分布著各類第四紀沉積物和粉沙質海岸地貌[9]。其中第四系沉積層厚度較大，平均200～300 m，沉積物成因類型以海相、湖沼相、河流沖積物為主； 入海河流形成的古河道帶透水性較強，成為咸水入侵的良好通道[24]。該地區(qū)資源豐富、人口密集、海洋經(jīng)濟發(fā)達，擁有得天獨厚的鹵水礦資源，是我國沿海最大的鹵水礦區(qū)，但也是我國最早發(fā)現(xiàn)海水入侵現(xiàn)象的地區(qū)之一[16]。受氣候變化和人類活動的共同影響，該區(qū)域已經(jīng)成為我國目前海水入侵最嚴重的地區(qū)，給當?shù)毓ぁ?農業(yè)生產和人民生活造成了極大危害。

圖1 研究區(qū)位置 (影像底圖為TM B4(R),B3(G),B2(B)假彩色合成影像)Fig.1 Location of study area

1.2 數(shù)據(jù)源

本文中使用的海水入侵數(shù)據(jù)主要根據(jù)魯勘字[2011]14號項目《黃河三角洲高效生態(tài)經(jīng)濟區(qū)(濰坊)海(咸)水入侵調查與監(jiān)控預警系統(tǒng)建設報告》和萊水字[2011]39號文件，以及相關文獻中的數(shù)據(jù)綜合整理得到[17,21,23]。對海水入侵數(shù)據(jù)整理、利用中，存在空間和尺度的誤差，需要進行空間配準和一致性檢驗，以去除彼此矛盾的數(shù)據(jù)。根據(jù)這些數(shù)據(jù)的采集時間，同時考慮遙感數(shù)據(jù)質量、受云層干擾情況等，去匹配相應時間的衛(wèi)星影像。最終，本文挑選6期(1979, 1990,1995,2000,2003和2013年)遙感數(shù)據(jù)用于監(jiān)測岸線變化(表1)。對遙感數(shù)據(jù)進行了預處理，將DN值轉換為光譜反射率[28-29]。

表1 用于海岸線監(jiān)測的遙感數(shù)據(jù)Tab.1 Remote sensing data used for shoreline monitoring

2 研究方法

本文的總體技術流程如圖2所示。

圖2 技術流程圖Fig.2 Technical flow chart

在多時相海水入侵數(shù)據(jù)整理和相應時段遙感影像選擇的基礎上，首先對海水入侵專題圖進行配準和數(shù)字化，整合1979―2012年6期海水入侵鋒線； 對遙感影像進行預處理，提取1979―2013年的海岸線位置分布信息。然后，采用縱向剖面分析(垂直斷面建模方法)[28]，每隔200 m采樣，獲取每個空間采樣點的定量變化信息； 再根據(jù)端點速率法(end point ratio，EPR)分析模型，分別計算和統(tǒng)計分析入侵鋒線、海岸線的時間分段變化速率及多年平均變化率; 最后，通過Pearson相關分析與顯著性t檢驗，分析海岸線變化與海水入侵在時間和空間上的耦合關系。

2.1 海岸線提取與鋒線數(shù)字化

通過遙感影像監(jiān)測海岸線的變化，在國內外已有大量的算法研究和應用。本文主要利用面向對象分類技術，根據(jù)岸線在遙感影像中的空間分布特征和水陸光譜差異，依據(jù)歸一化差值水體指數(shù)(normalized difference water index，NDWI)類間光譜特征差異最大原理，完成水陸分離； 然后通過人工互信息操作，實現(xiàn)海岸線的高精度識別，具體技術流程參考文獻[28]。首先采用大津算法(OTSU)完成NDWI水陸分離； 然后依據(jù)空間關系特征判斷實現(xiàn)陸地水域與海洋判別； 最后通過互信息人工操作對海岸線提取結果進行優(yōu)化處理和微地形修正，得到研究區(qū)的各期海岸線分布圖。互信息后處理操作主要包括： 河口處理需保留大型河口港灣特征和小河齊陸地基線； 瀉湖需考慮其內環(huán)境，當瀉湖與海洋有寬闊水域通道時，在瀉湖內部采集海岸線。海水入侵鋒線是在幾何精確配準的基礎上，對歷史時期的海水入侵鋒線進行手工數(shù)字化，得到各個歷史時期海水入侵鋒線的空間分布位置。

需要說明的是，本文中的海岸線數(shù)據(jù)沒有考慮平均大潮位修正，不是真正意義的海岸線； 準確地說是利用遙感影像提取的瞬時水邊線。理論上做潮位修正效果可能會更好，但因這個區(qū)域地表變化很大，SRTM高程數(shù)據(jù)在這個區(qū)域不是特別準確，潮位修正后的結果和遙感影像套合不上。所以，綜合考慮到本文中分析數(shù)據(jù)是通過EPR 模型計算多個時期時間序列岸線的平均變化速率，盡管用水邊線代替海岸線時潮位會對分析結果會產生一定的影響(不可避免)，但基本可以近似地反映區(qū)域海岸線的動態(tài)變化趨勢(與其他文獻中描述的趨勢基本一致)，具有一定的可比性，對分析海水入侵鋒線位置影響不大。

2.2 縱向斷面建模方法

垂直斷面建模方法是一種對數(shù)字岸線線性目標變遷速率計算傳統(tǒng)的方法，其中以端點速率法 (EPR)計算模型最為簡單常用[28]。EPR模型統(tǒng)計計算岸線的變化速率原理如圖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點之間的距離。

首先建立垂直剖面線，本文在提取海水入侵鋒線和岸線的基礎上，通過緩沖區(qū)分析建立基線，設置采樣間隔為200 m； 然后根據(jù)垂線與所有斷面上的岸線交點距離，計算平均變化率。那么，對于2個歷史時期的岸線變化速率，可直接用2個歷史時期的變遷距離除以時間來計算，EPR1,2的數(shù)學表達式為

(1)

式中：D1，D2分別為time1，time2兩個歷史時期岸線的位置；distance(D1-D2) 為D1，D2兩點之間的距離。

對于多個時期的一組時間序列岸線的平均變化速率，EPRave的數(shù)學表達式為

(2)

式中：timelatest為最新歷史時間；timeoldest為最老歷史時間；Dlatest,Doldest分別為timelatest,timeoldest兩個歷史時期岸線的位置。

3 結果分析

3.1 海水入侵鋒線動態(tài)變化過程

1979―2012年萊州灣海岸帶區(qū)域海水入侵鋒線位置動態(tài)變化見圖4。從咸水入侵鋒線圖可以清楚地看出區(qū)域海水入侵的動態(tài)變化過程。1979―1989年是區(qū)域海水入侵最快的階段，隨后區(qū)域海水鋒線移動明顯變緩。圖5為海水入侵EPR時空動態(tài)分布圖。綜合分析海水入侵的動態(tài)過程和EPR時空動態(tài)分布圖，可以從總體上將區(qū)域海水入侵分為3個階段： ①1979―1989年的快速入侵階段。此階段是萊州灣南岸海水入侵發(fā)展最快的時期，每年入侵面積增加29.22 km2； 海咸水平均每年入侵速率為380 m/a，最大入侵距離在昌邑市青鄉(xiāng)南部，入侵速率達到690 m/a； ②1989―2000年的慢速發(fā)展階段。此階段萊州灣南岸海水入侵的速度明顯變緩，無論是入侵速度，還是入侵面積，都明顯下降。海水入侵鋒線基本上在1995年的位置附近波動； ③2000―2012 年的相對穩(wěn)定與徘徊階段。此階段海咸水入侵動態(tài)鋒線基本處于徘徊狀態(tài)，并且局部出現(xiàn)了退縮(圖5(e))。

圖4 1979―2012年萊州灣南岸海水入侵鋒線動態(tài)過程 (影像底圖為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空間動態(tài)分布圖

Fig.5SpatialdynamicdistributionofEPR

3.2 海岸線動態(tài)變化遙感制圖與分析

從解譯的結果可以看出，萊州灣南部岸線的時空變化并不是一致的(圖6)。

圖6 1979―2013年區(qū)域海岸線變遷圖 (影像底圖為MSS B3(R),B2(G),B1(B)假彩色合成影像)Fig.6 Regional coastline transition from 1979 to 2013

從空間上看，區(qū)域海岸線總體上呈現(xiàn)向陸地蝕退的趨勢(局部區(qū)域受到人類活動的干擾(如人工修建海上設施)除外)； 從時間上看，不同地點區(qū)域海岸線的蝕退速率也不一樣，相同地點不同時段的蝕退速率也不同。

3.3 海水入侵與岸線變化耦合關系分析

1979―2012年萊州灣南岸的海水入侵與岸線變化速率見圖7。從圖7可以看出2條EPR曲線的變化特征和時空耦合關系。在空間分布格局上，海水入侵的強度和海岸線的蝕退強度相關性具有很強的一致性； 在時間變化速率上，海水的入侵速率變化曲線和海岸線變化的變化速率曲線形狀相似。這一變化特征說明，萊州灣南岸地區(qū)海水入侵和海岸線蝕退在空間分布上和時間變化上均具有一定的相關性。

(圖中Y軸表示基準岸線，正值表示向海方向，負值表示向陸地方向)圖7 1979―2012年萊州灣南岸岸線與海水入侵鋒線變化速率圖Fig.7 Change rates of coastline and seawater intrusion front line of Laizhou Bay from 1979 to 2012

表2從定量的角度統(tǒng)計了區(qū)域岸線蝕退和海水入侵的變化特征以及二者之間的相關關系。表中海岸蝕退正值表示向海方向淤進，負值表示向陸地方向蝕退； 海水入侵負值表示向陸地方向侵入； ** 表示在0.01水平(雙側)顯著相關。

表2 岸線蝕退與海水入侵變化量統(tǒng)計與耦合關系分析Tab.2 Relationship analysis between seawater intrusion and coastline changes (m/a)

1979―2012年萊州灣南岸海水入侵速率均值達到了 177.23 m/a； 最大值分布在昌邑市青鄉(xiāng)南部區(qū)域，入侵速率達到 435.28 m/a； 最小值分布在濰坊北部寒亭,為 124.02 m/a。1979―2012年區(qū)域岸線的蝕退速率約為10.44 m/a； 蝕退的最大速率為124.02 m/a，空間位置上分布在昌邑市北部下營鎮(zhèn)； 蝕退的最小速率分布在寒亭區(qū)以北的岸線段。此處由于人工建造海上游樂場和港口，區(qū)域岸線呈現(xiàn)向海洋快速淤進，最大速度達到191.14 m/a。海岸線EPR與入侵鋒線EPR相關系數(shù)為0.407，顯著性水平P值為0.00，通過99%的信度水平(雙側)檢驗。由此可見，區(qū)域海水入侵與岸線蝕退在時空上存在強耦合關系。

4 結論與討論

4.1 結論

本文對1979―2012年萊州灣南岸海水入侵鋒線演化的時空過程進行數(shù)字重建和區(qū)域海岸線的多期遙感監(jiān)測； 通過EPR定量模型探索性地分析了二者之間的時空耦合關系機制。研究得到以下結論：

1)區(qū)域海水入侵經(jīng)歷了由快到慢的變化過程，1990年以后速率明顯減緩，入侵鋒線基本穩(wěn)定在1995年鋒線附近； 2008―2012年入侵鋒線局部發(fā)生后退，目前保持相對穩(wěn)定。這一研究結論與之前相關學者的研究結果基本一致[22]。

2)除了局部人工造陸導致岸線向海擴張外，1979―2012年區(qū)域海岸線以蝕退型海岸為主，但是蝕退速度在空間上差異明顯。

3)海水入侵鋒線變化與海岸線進退二者之間在時空上存在強耦合關系，相關系數(shù)達到0.407，顯著性水平P<0.01(雙側)，說明海岸線工程對區(qū)域海水入侵速率有顯著的影響。以上研究結論可為區(qū)域海水入侵的預防和治理提供數(shù)據(jù)支撐和科學依據(jù)。

4.2 討論

1)在變化原因上，影響咸淡水界面遷移速率的因素可分為自然因素和人為因素，自然因素主要是地質條件、大氣降水減少和風暴潮等； 人為因素主要是地下水開采及鹵水的開采。特殊的地質條件(廣泛分布的第四紀海相、湖相沉積以及河流沖積物，透水性強)為研究區(qū)內海咸水入侵提供了有利條件； 而氣候變暖及海平面升高是該區(qū)海咸水入侵的大環(huán)境背景，人類活動則改變了海水入侵時空變化進程和入侵強度[9,24]。

2)在變化特征上，1979―1989年海咸水快速入侵階段，主要是由于地下淡水資源的過渡開采，地下水位大幅度下降； 1989―2000年海咸水慢速發(fā)展階段，由于當?shù)刂饾u認識到海咸水入侵的危害，對地下水開采進行了適當限制，同時近海鹵水的開采量在逐年增大； 2000―2012 年海咸水相對穩(wěn)定與后退階段，主要是得益于海水入侵治理工程得實施，通過協(xié)調和控制地下水與鹵水的開采，較好地控制了咸淡水界面。因此，按照研究區(qū)咸淡水界面變化情況，大量的研究分析共同表明： 地下水資源開發(fā)和鹵水資源開采起了決定性作用，是區(qū)域海水入侵的直接原因[3,19, 22-23]。

3)海岸線的進退是影響區(qū)域海水入侵的又一疊加因素。海岸線的淤進，尤其是人工修筑的岸線，在一定程度上可以阻止海水入侵。可能的原因是人工修筑的岸線有利于阻擋風暴潮，堵塞海水入侵的部分通道。

以上研究結論建立在文中數(shù)據(jù)和方法的基礎上。由于EPR模型方法計算線性目標的變化率是通過剖面上兩點距離與時間的比值得到的，對于短期的加速、緩慢或者逆轉變化不夠敏感。尤其是在計算多期平均變化率時，短期的逆轉性變化對整體變化趨勢會產生一定的影響。此外，多源數(shù)據(jù)的空間精度誤差問題、時間尺度問題不可回避。在空間尺度上，遙感數(shù)據(jù)源空間分辨率為30～60 m，在提取的海岸線精度上難免有些誤差； 海水入侵歷史數(shù)據(jù)資料的幾何變形，導致局部精度難以控制。在時間尺度上，受到觀測數(shù)據(jù)的限制，通過6期數(shù)據(jù)分析變化過程，時間分辨率精度不夠高，對演化進程分析尤其是時間拐點的確定不夠準確。這些可能的誤差會對數(shù)據(jù)分析以及部分研究結論產生一定的影響。在未來的研究中有些方面將有待于進一步改進和提高，例如加大海水入侵鋒線動態(tài)監(jiān)測時間密度，改進算法提高對短期動態(tài)變化的敏感性等。