袁德奎，李 廣，王道生，楊志斌

(天津大學機械工程學院，天津 300072)

渤海灣位于渤海西部，是一個典型的半封閉緩坡淤泥底質(zhì)淺水海灣，波浪潮流作用下污染物沿岸輸移趨勢明顯[1]，水體交換能力較弱，陸源排放的污染物難以運移到渤海中部或外海，從而造成灘涂和近岸海域的嚴重污染．隨著環(huán)渤海經(jīng)濟的快速發(fā)展，一方面，頻繁的人類活動改變了近岸海域的動力條件，進而影響到渤海灣的水交換能力；另一方面，陸源污染物排放量的增加，使渤海灣的生態(tài)環(huán)境承受著巨大的壓力．為了保障該區(qū)域經(jīng)濟與生態(tài)環(huán)境的協(xié)調(diào)發(fā)展，對渤海灣的水交換開展深入的研究，認識渤海灣的水交換特性是非常必要的．許多學者對渤海(含渤海灣)的水交換能力進行了研究[2-4]．

縱觀前人的研究，多采用水質(zhì)模型來獲得區(qū)域的整體水交換特性，得到的整體半交換時間各不相同，原因可能是采用的地形各不相同，或是渤海灣開邊界的選取不同及考慮的分潮有所差異．隨著近年來渤海灣圍填海工程項目的相繼展開，渤海灣岸線和地形發(fā)生顯著變化，有必要進一步研究渤海灣的水交換特性．對于渤海灣這樣的大型海灣，灣內(nèi)各處的水體難以在短時間內(nèi)混合均勻，因此，不僅需要分析其整體的水交換能力，也需要分析其局部水交換特性．此外，對于像渤海灣這樣以往復型潮汐為主要驅(qū)動力的海灣，回流對水交換的影響可能較大，但目前關(guān)于此方面的研究較少．目前存在多種表征水交換能力的時間尺度，如半交換時間、駐留時間、曝光時間和生命等．不同時間尺度的計算基于的假設不盡相同，采用不合理的尺度可能導致低估或高估系統(tǒng)的水交換能力[5-6]，選擇一個合適的時間尺度來評價系統(tǒng)的水交換能力是非常重要的．

本文采用拉格朗日粒子跟蹤方法，從整體和局部層面分析了渤海灣水交換特性：根據(jù)數(shù)值模擬結(jié)果分析了渤海灣的整體水交換能力；通過統(tǒng)計各計算節(jié)點的曝光時間、駐留時間[7]以及返回系數(shù)，探討了水體返回對渤海灣水交換能力的影響；進一步通過情景模擬分析了圍填海工程對渤海灣水交換能力的影響．

1 水交換時間尺度及計算方法

在研究水體的水交換能力時，需要考慮兩個問題：①該水體和外界水體的整體交換能力；②該水體各個子區(qū)域的水交換能力．對于前者，可用具有整體特性的時間尺度來表征，如半交換時間；對于后者，可用具有局部特性的時間尺度來表征，如駐留時間和曝光時間．

1.1 駐留時間和曝光時間

半交換時間作為綜合指標，可以在一定程度上表征所研究區(qū)域的整體水交換能力．但對于水交換能力存在明顯的空間差異的水體，采用這樣的方法來評估系統(tǒng)的水交換能力并不充分，可以采用駐留時間或曝光時間來表征各點的水交換能力，同時也可以通過駐留時間和曝光時間相結(jié)合來評價當?shù)厮w的回流情況．

駐留時間通常定義為自關(guān)心時刻起，水團第1次離開感興趣的區(qū)域所需要的時間，每個水團的駐留時間都是唯一的．由定義可知，駐留時間能體現(xiàn)出空間的差異性，但沒有考慮到水團的回流對水體交換的影響．而曝光時間包括了回流．若水團離開感興趣的區(qū)域后性質(zhì)發(fā)生了顯著的變化，顯然駐留時間更適合來描述區(qū)域內(nèi)的水交換特性；但是，若區(qū)域的邊界是人為劃分的，并沒有顯著的理化特性的區(qū)別，則曝光時間更為可取，因為這時駐留時間會低估污染物在區(qū)域內(nèi)的影響．

1.2 返回系數(shù)

返回系數(shù)[8]b（x, y ）可表示為

式中：τ（ x , y） 為曝光時間；τ0（ x, y ）為駐留時間．

由式(1)可以看出，返回系數(shù)為 0表示曝光時間等于駐留時間，即粒子一旦離開感興趣的區(qū)域就不再返回；返回系數(shù)大于 0表示曝光時間大于駐留時間，即粒子在第1次離開感興趣區(qū)域后曾返回；返回系數(shù)接近于1表示曝光時間遠大于駐留時間，即粒子曾多次返回感興趣區(qū)域．可見，返回系數(shù)可以用來表征水團離開區(qū)域后返回能力的強弱，進而可以說明區(qū)域內(nèi)各點對整體水交換的貢獻．返回系數(shù)越大，表示該區(qū)域?qū)φw水交換能力貢獻越?。?/p>

2 數(shù)學模型

2.1 水動力學及隨機游動模型

渤海灣水交換時間的計算模型由兩部分組成，分別是水動力學模型和隨機游動模型，前者為后者提供基本的流場信息，后者用于計算粒子的運動軌跡，從而進一步統(tǒng)計出渤海灣的水交換時間尺度．

對于渤海這樣的淺水系統(tǒng)，垂向運動的尺度遠小于水平尺度，采用沿水深積分的二維模型不僅可以獲得較滿意的結(jié)果，還可以提高計算效率．沿水深積分的二維水動力數(shù)學模型的控制方程[4]如下所述．

連續(xù)性方程為

動量方程為

式中：t為時間；ζ為水位；ui為沿水深平均的流速在i方向的分量；pi為流體在i方向的單寬通量；H為總水深，H=ζ+h，h為相對于基準面的水深；Sm、ui分別為單位水平面積上的源項強度和該源項初始速度在i方向的分量；β為動量修正系數(shù)；f為科氏力系數(shù)；g為重力加速度；ρ為水體密度；τi為自由水面剪切應力在i方向的分量；σi為水體底部剪切應力在i方向的分量；ε為水深平均渦黏系數(shù)．

在空間交錯網(wǎng)格系統(tǒng)上用 ADI方法對控制方程進行差分求解，其中對流項用二階迎風格式處理．

在每個計算步，用隨機游動模型[9]求解粒子的運動軌跡．

2.2 模型設置

為了分析粒子回流的影響，模型的開邊界取在遠離渤海灣的大連和煙臺連線．模型的計算區(qū)域覆蓋117.5°E～123°E，37°N～41°N，如圖 1所示．

圖1 渤海地形、水位開邊界及渤海灣邊界Fig.1 Topography of Bohai Sea with open boundary and boundary of Bohai Bay

在經(jīng)度和緯度方向的空間步長分別為 0.025°和0.020°，時間步長為 60,s．用大連和煙臺兩個驗潮站的 8 個主要分潮(M2、S2、K1、O1、N2、P1、K2 和Q1)的調(diào)和常數(shù)來插值生成開邊界上各計算節(jié)點上的水位，以驅(qū)動整個研究區(qū)域內(nèi)的潮流運動．模型計算采用冷啟動，初始水位和流速均為零．

初始時刻在渤海灣內(nèi)根據(jù)當?shù)厮畈荚O示蹤粒子，每1,m深度布置1層粒子，小于1,m的地方只布置1層，新地形下共布置19,941個粒子，舊地形下布置 19,773個粒子．文獻[10]的研究表明，在不同的時刻(潮相位)釋放示蹤粒子，統(tǒng)計得到的水交換特征時間存在一定差異．為了分析潮相位對水交換時間的影響，設計了表1所示的8個算例，數(shù)值模擬中分別在4個特征時刻釋放示蹤粒子，粒子跟蹤的計算時間為 3年．為了進一步說明渤海灣內(nèi)新舊地形下各處的水交換路徑，選取10個代表點(選取的理由在后文有具體分析)來描繪粒子在 1年內(nèi)的輸移軌跡，各代表點的初始位置見圖 2(新舊岸線條件下粒子初始位置相同)．

表1 算例描述Tab.1 Description of simulation cases

圖 2為本文所對比的新舊地形(分別是根據(jù)2010年和 2003年海圖數(shù)字化得到)，其中黑色區(qū)域為 2003年的岸線邊界，灰色區(qū)域為 2010年的岸線邊界．

圖2 岸線變化及代表點初始位置分布Fig.2 Changes of coastal line and initial locations of representative points

3 結(jié)果及討論

3.1 模型驗證

本文所用水動力學模型已多次用于渤海流場的模擬，結(jié)果與實測數(shù)據(jù)符合良好，模型的具體驗證可參見文獻[3，11-12]．為保證數(shù)值模擬結(jié)果的可靠性，本文采用2003年7月13日14時至16日12時的71,h對B2站位的實測資料以及2012年11月15日16時至16日16時一晝夜A1站位的實測潮位和流速資料對模型在新舊地形下的模擬結(jié)果進行了進一步驗證，測點位置如圖 2所示，對比結(jié)果如圖 3所示．總體來看，計算值和實測值吻合良好，模型能較好地模擬出渤海灣的流場．

圖3 計算值和實測值的對比Fig.3 Comparison between calculated results and measured data

3.2 圍填海工程導致的渤海灣余流場變化

為了分析渤海灣水交換特性的內(nèi)在機制，統(tǒng)計了新舊地形下渤海灣的余流場，即

式中：Q為各點的水深積分流量；rQ為各點的余流流量；aT為統(tǒng)計時間(1年)．

圖4給出了新舊地形下渤海灣的余流場，可以看出地形變化對余流的影響主要體現(xiàn)在：天津港北部的逆時針環(huán)流消失；天津港南部圍海造地和黃驊港北部工程的擴建，導致該地區(qū)的沿岸流減??；由于黃驊港的阻隔作用，流線發(fā)生變形向深水擴張；曹妃甸淺灘處的余流基本消失，只有西南角的順時針環(huán)流依然存在．以上結(jié)果與秦延文等[13]的結(jié)論基本一致．無論新舊地形，大沽口以北、曹妃甸以南為順時針環(huán)流，大沽口以南存在一個逆時針環(huán)流；灣口北部水域有一支流入渤海灣的余流，但是深入渤海灣不遠即左轉(zhuǎn)向東南偏轉(zhuǎn)，與渤海灣南岸的順時針環(huán)流匯集從渤海灣中部流出；即存在一個雙渦旋結(jié)構(gòu)，這與已有的研究結(jié)果[4]基本一致．

圖4 渤海灣余流場Fig.4 Residual current field of Bohai Bay

3.3 岸線變化對渤海灣整體水交換的影響

圖5 (算例描述見表1)給出了渤海灣內(nèi)剩余的粒子數(shù)與初始時刻釋放的粒子數(shù)的比例(本文將該比例稱為無量綱剩余平均濃度)隨時間的變化過程，實線和虛線分別為新舊岸線條件下的情況．可以看出不同時刻釋放的粒子的無量綱剩余平均濃度都隨著時間的推移而減?。?2給出了根據(jù)無量綱剩余平均濃度統(tǒng)計結(jié)果計算出的渤海灣水體的半交換時間．顯然，粒子釋放時刻的不同對半交換時間的計算結(jié)果有一定影響．在舊岸線條件下，渤海灣整體的半交換時間在 300,d左右，這與已有的研究結(jié)果[2-3]吻合；在新岸線條件下，渤海灣整體的半交換時間超過600,d．

圖5 新舊岸線條件下不同時刻釋放的粒子的無量綱剩余平均濃度隨時間的變化Fig.5 Variation of dimensionless average remaining concentration of particles released at different tide phases with time under new and old coastal line respectively

表2 渤海灣水體半交換時間及交換率的計算結(jié)果Tab.2 Numerical results of semi-exchange time and exchange rate of Bohai Bay

圖6為新舊岸線條件下高潮時刻(表1中的算例1和算例5)釋放的示蹤粒子運動1年后的分布，其余算例的整體分布形式與之類似．粒子的分布與余流場有很強的相關(guān)性：渤海灣中部和西北部的環(huán)流是造成這兩處大量粒子駐留在灣內(nèi)的主要原因；有兩股粒子流出了渤海灣，其中一支沿黃驊港南部流出渤海灣的余流流向膠州灣方向，另一支由于與渤海中部相連的灣口處環(huán)流作用則徘徊在渤海中部．渤海中部徘徊的粒子在渤海灣外平均停留了1年以上，水體已得到凈化，因此不會對渤海灣整體的自凈能力產(chǎn)生明顯的負面影響．

對比新舊岸線條件下粒子的遷移擴散范圍以及對新舊岸線條件下的余流場，可以推斷出，由于曹妃甸圍填海工程的建設使得渤海灣北部沿岸余流消失，加之渤海灣灣口中部環(huán)流的作用，新岸線條件下灣口北部粒子更多地向南輸移，沿南部余流流出渤海灣，而不是直接通過北部流出，降低了該處水體的水交換能力．而黃驊港的阻隔作用，使渤海灣西岸的沿岸余流有所減弱，造成南部沿岸流出渤海灣的粒子遷移擴散能力變?nèi)?，從而使新岸線條件下渤海灣整體的半交換時間變長．

圖6 高潮時刻釋放的粒子1年后的分布Fig.6 Distributions of particles released at high tide after one year

3.4 地形變化對渤海灣局部水交換的影響

為進一步討論岸線變化對渤海灣局部水交換的影響，引入無量綱駐留時間和曝光時間來定量評價渤海灣各處水體對整體水交換的貢獻．定義無量綱駐留時間R和無量綱曝光時間E為

式中T為總的計算時間，1年．

圖 7給出了在新舊地形下渤海灣內(nèi)各點高潮時刻釋放的粒子的無量綱駐留時間、無量綱曝光時間和返回系數(shù)的分布．與粒子的分布類似，粒子釋放時刻的不同對駐留時間、曝光時間和返回系數(shù)分布的細節(jié)雖有一定影響，但對其整體分布形式影響不大，故僅給出了高潮時刻釋放粒子的計算結(jié)果．

在原岸線條件下，渤海灣中部、西部地區(qū)以及黃驊港東側(cè)的近岸區(qū)域?qū)Σ澈痴w的水交換貢獻較小．從圖 4的余流場中可以看出，渤海灣中部環(huán)流一方面導致大量粒子在此處徘徊，另一方面阻礙西部粒子離開渤海灣；南部余流較弱，加之有一股支流沿岸從萊州灣流向渤海灣，阻擋粒子沿岸東移路徑，所以粒子難以流出渤海灣，造成以上地區(qū)駐留時間和曝光時間都較長，且二者相近，返回系數(shù)較低．而黃驊港東北部近岸區(qū)域?qū)Σ澈痴w水交換的貢獻較大，其原因是該區(qū)域存在一股向東方向的余流，使得該區(qū)域的粒子可以比較順利地沿灣口南部流出渤海灣，并且很少返回，因此，該區(qū)域的駐留時間和曝光時間并不太長，且返回系數(shù)較?。澈碁晨谀喜狂v留時間和曝光時間都短，表明這里的粒子在釋放后很快會流出渤海灣且基本不再返回，對渤海灣整體的水交換貢獻大．在灣口北部，雖然駐留時間也短，但曝光時間很長，表明粒子在流出渤海灣后很快被帶回渤海灣，并在灣內(nèi)長期徘徊，因此，灣口北部的水體對渤海灣整體水交換的貢獻小，這主要是曹妃甸南岸的回流所致.

對比圖4中新舊岸線條件下的余流場可以看出，岸線變化對渤海灣中部環(huán)流以及南部余流的影響較小，因此，在新地形下，渤海灣中部、西岸以及南岸水體對渤海灣整體水交換的貢獻依然較小．但明顯不同的是：①在天津港和曹妃甸建設以前，渤海灣西北部及曹妃甸沿岸處水體的駐留時間短，在1年時間內(nèi)至少有 1次離開過渤海灣，與外界水體交換通暢；在新岸線條件下，雖然天津港北部的逆時針環(huán)流消失，但是西北部的環(huán)流依然存在，加之曹妃甸甸頭的阻隔作用以及曹妃甸東向余流的消失，導致駐留時間和曝光時間都較長，返回系數(shù)小，表明該處水體不易流出渤海灣，對渤海灣整體的水交換貢獻??；②灣口北部水體的駐留時間和曝光時間變長、返回系數(shù)變小，說明該處水體對渤海灣整體水交換的貢獻減弱，原因是北部沿岸淺灘處余流因為曹妃甸的建設而消失，使該處的水體難以直接從灣口北部離開渤海灣；在灣口南部，駐留時間和曝光時間變短，返回系數(shù)變小，而且作用范圍變大，說明岸線變化加強了灣口南部水體的水交換能力，該處水體更快地離開渤海灣，且不再返回．

對返回系數(shù)做進一步的討論．從圖 7可以看出：在灣口南部的粒子，離開渤海灣后基本不再返回，駐留時間和曝光時間相近，返回系數(shù)較低，對區(qū)域的整體水交換貢獻較大；而對于灣口北部的粒子，其返回系數(shù)也較低，但粒子并未離開過渤海灣，對區(qū)域的整體水交換貢獻?。@然同樣的返回系數(shù)對應著兩種截然相反的情況，不能單獨依靠返回系數(shù)正確地評價各處粒子對整個區(qū)域水交換的貢獻，而需結(jié)合駐留時間和曝光時間進行分析．若返回系數(shù)為 0，則用駐留時間就足以描述水交換時間；但是若返回系數(shù)不為0，則需綜合考慮駐留時間和曝光時間來評價水交換能力．

圖7 高潮時刻釋放的粒子的無量綱駐留時間、無量綱曝光時間和返回系數(shù)(算例1和算例5)Fig.7 Dimensionless residence time，exposure time and return coefficient of particles released at high tide (case 1 and case 5)

圖8 2003年與2010年地形條件下粒子1年內(nèi)的運動軌跡Fig.8 Particle trajectories in a year under bathymetry of 2003 and 2010，respectively

根據(jù)圍填海工程對局部水交換能力的影響，本文選取分布在天津港、黃驊港、渤海灣中部及渤海灣灣口的示蹤粒子，進一步分析圍填海對渤海灣各處的水交換路徑的影響．圖 8給出了新舊地形下不同區(qū)域高潮時刻釋放的粒子1年內(nèi)的運動軌跡，圖中數(shù)字及連線代表粒子運動的先后順序及軌跡．可以發(fā)現(xiàn)岸線的變化使北塘口附近的 1號粒子在西北部順時針徘徊．2號粒子布置在新地形下的天津港出?？诟浇?，雖然粒子向東南方向徑直運動，但至少在灣內(nèi)停留了1年．從3號、5號和6號粒子在新地形下的運動軌跡也可以看出新地形下渤海灣南岸沿岸余流減弱，粒子離開渤海灣受阻，在渤海灣灣口輸移路徑變短．雖然渤海灣中部的 7號粒子的輸移路徑有所改變，但是長期停留在灣內(nèi)，且運動路徑較短，說明地形變化對這些地區(qū)的水交換影響不大．而在舊地形下，位于黃驊港東部的 4號粒子沿岸向北曲折前進，在新地形下可能受到黃驊港的阻隔作用，轉(zhuǎn)而向東北方向運動，朝灣口輸移．8號粒子向灣口運動的趨勢更加明顯，運動軌跡變長，原因可能是受新地形影響灣口北部原東向余流轉(zhuǎn)而向南流動以及南部西向余流受到黃驊港的阻隔而轉(zhuǎn)回東向，兩者共同作用加強了灣口南部的余流強度．從 10號粒子的運動軌跡也可以看出灣口南部的水交換能力在新地形下變強.而位于灣口北部的 9號粒子，雖然可能離開過渤海灣，但是整體向南運動的趨勢更加明顯，最后由南部灣口離開，所以在新地形下北部的水交換能力變?nèi)?，原因可能是北部曹妃甸淺灘處余流的消失，使北部水體向南輸移，需要更長的時間才能與外界進行交換.

4 結(jié) 論

根據(jù)拉格朗日粒子跟蹤法計算的結(jié)果，統(tǒng)計出了大規(guī)模圍填海前后渤海灣水體的半交換時間，從整體層面上評估了圍填海工程對渤海灣水交換能力的影響；進一步統(tǒng)計了渤海灣各處水質(zhì)點的駐留時間、曝光時間和返回系數(shù)，從局部層面上分析了渤海灣各處水體對其整體水交換的貢獻；結(jié)合新舊地形下余流場的變化分析了水交換能力變化的內(nèi)在機理；為進一步闡釋地形變化對渤海灣水交換特性的影響，選取典型位置的粒子作為標志點，分析了地形變化對粒子運動軌跡的影響．

(1) 渤海灣海岸帶圍填海工程的建設引起的岸線變化，使得天津港北部的逆時針環(huán)流消失，曹妃甸西南角的順時針環(huán)流依然存在，但是淺灘處的沿岸余流消失，由于受到黃驊港的阻隔，渤海灣西岸沿岸向東南方向余流有所減弱．

(2) 圍填海工程使渤海灣整體水交換能力變?nèi)酰虢粨Q時間由原先的300,d左右變?yōu)槌^600,d．

(3) 渤海灣西部、中部和南部沿岸地區(qū)對渤海灣整體水交換的貢獻??；黃驊港東北部水體對渤海灣整體水交換的貢獻較大，且以上地區(qū)的水交換能力受地形變化的影響較??；水交換能力受地形影響變化較大的地區(qū)有渤海灣西北部及灣口，西北部水體對整體的水交換貢獻更小，灣口北部的水交換能力變?nèi)?，南部的水交換能力增強．

(4) 對于渤海灣灣口這樣人為劃定的邊界，由于存在明顯的余環(huán)流，該處粒子的返回系數(shù)可能較高，用返回系數(shù)分析其局部水交換特性及其對整體水交換的貢獻時需結(jié)合駐留時間和曝光時間．

［1］ Sun Tao，Tao Jianhua. Experimental and numerical study of wave-induced long-shore currents on a mild slope beach[J]. China Ocean Engineering，2005，19(3)：469-484.

［2］ 魏 皓，田 恬，周 鋒，等. 海灣水交換的數(shù)值研究：水質(zhì)模型對半交換時間的模擬[J]. 青島海洋大學學報，2002，32(4)：519-525.Wei Hao，Tian Tian，Zhou Feng，et al. Numerical study on the water exchange of the Bohai Sea：Simulation of the half-life time by dispersion model[J].Journal of Ocean University of Qingdao，2002，32(4)：519-525(in Chinese).

［3］ Sun Jian，Tao Jianhua. Relation matrix of water ex-change for sea bays and its application[J]. China Ocean Eng，2006，20(4)：529-544.

［4］ 王金華，沈永明，石 峰，等. 基于拉格朗日粒子追蹤的渤海冬季與夏季環(huán)流及影響因素[J]. 水利學報，2011，42(5)：544-553.Wang Jinhua，Shen Yongming，Shi Feng，et al. Study of Bohai Sea circulation and its influencing factors during winter and summer based on Lagrangian particle tracing method[J]. Journal of Hydraulic Engineering，2011，42(5)：544-553(in Chinese).

［5］ Takeoka H. Fundamental concepts of exchange and transport time scales in a coastal sea[J]. Continental Shelf Research，1984，3(3)：311-326.

［6］ Zimmerman J T F. Mixing and flushing of tidal embayments in the Western Dutch Wadden Sea(Part I)：Distribution of salinity and calculation of mixing time scales[J]. Netherlands Journal of Sea Research，1976，10(2)：149-191.

［7］ de Brauwere A，de Brye B，Blaise S，et al. Residence time，exposure time and connectivity in the Scheldt Estuary[J]. Journal of Marine Systems，2011，84(3)：85-95.

［8］ Feleke A，Straud A，Badr A W. Modeling of residence time in the East Scott Creek Estuary[J]. Hydro-Environment Research，2008，2(2)：99-108.

［9］ Suh S W. A hybrid approach to particle tracking and Eulerian-Lagrangian models in the simulation of coastal dispersion[J]. Environmental Modelling & Software，2006，21(2)：234-242.

［10］ Yuan D，Lin B，F(xiàn)alconer R A. A modeling study of residence time in a macro-tidal estuary[J]. Estuarine，Coastal and Shelf Science，2007，71(3/4)：401-411.

［11］ Li X，Yuan D，Sun J. Simulation of water exchange in Bohai Bay[C]//Advances in Water Resources and Hydraulic Engineering. Nanjing，China，2009：1341-1346.

［12］ Yuan D，Lin B，Tao J，et al. Verification of a numerical model using field monitoring data for modelling Bohai Bay[C]// Proceedings of the 28th IAHR Congress.Graz，Austria，1999：346.

［13］ 秦延文，張 雷，鄭丙輝，等. 渤海灣岸線變化(2003—2011年)對近岸海域水質(zhì)的影響[J]. 環(huán)境科學學報，2012，32(9)：2149-2159.Qin Yanwen，Zhang Lei，Zheng Binghui，et al. Impact of shoreline changes on the costal water quality of Bohai Bay(2003—2011)[J]. Acta Scientiae Circumstantiae，2012，32(9)：2149-2159(in Chinese).