潘 昀，竇希萍，焦增祥，王玉丹，丁 磊

(南京水利科學(xué)研究院 港口航道泥沙工程交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029)

自20 世紀(jì)50年代以來，圍繞長江口水流泥沙問題開始了大量研究工作，早期主要通過水文泥沙與河床演變分析、水槽和物理模型試驗(yàn)研究長江口泥沙運(yùn)動(dòng)機(jī)理，80年代末由于數(shù)值模擬技術(shù)的發(fā)展，建立了長江口潮流、波浪、泥沙數(shù)學(xué)模型，在河口形態(tài)［1］、攔門沙［2］、最大渾濁帶［3］、泥沙絮凝與沉降［4］、懸沙垂向分布［5］、泥沙擴(kuò)散［6］、浮泥［7］、航道回淤預(yù)測［8］等數(shù)學(xué)模型方面取得了突破性進(jìn)展，為長江口深水航道治理工程建設(shè)提供了技術(shù)支撐。長江口深水航道治理工程采用整治和疏浚相結(jié)合，分三期實(shí)施，1998年至2005年完成了一、二期工程，包括建設(shè)南北槽分流口工程;長分別為48 km 和49 km 的南、北雙導(dǎo)堤;總長約30 km的19 座丁壩。主要整治建筑物平面布置見圖2。一期工程完成后北槽航道水深達(dá)8.5 m，二期工程完成后水深達(dá)到10 m，三期工程主要為疏浚，自2006年9月開始疏浚，通過加長丁壩、增加疏浚船只，終于在2010年3月使北槽航道水深達(dá)到12.5 m。針對三期工程航道回淤嚴(yán)重的狀況，劉猛等［9］研究導(dǎo)堤越堤水沙對航道回淤影響;陳維等［10］研究南沙頭通道及橫沙通道的影響;沈琪等［11］分洪季特殊月對長江口深水航道的影響以及泥沙來源進(jìn)行分析;張宏等［12］基于神經(jīng)網(wǎng)絡(luò)方法研究深水航道海域懸浮泥沙濃度分布規(guī)律。由于北槽深水航道回淤涉及因素眾多，現(xiàn)有理論和模擬方法還難于解釋和描述這一現(xiàn)象，因此，利用不斷積累的北槽航道實(shí)測回淤資料，采用合適的方法進(jìn)行預(yù)測分析不失為一條途徑。

分形插值是分形幾何理論及其應(yīng)用研究中的一個(gè)重要內(nèi)容，它能很好地刻畫現(xiàn)實(shí)世界中普遍存在的粗糙的、復(fù)雜的、自相似、或自仿射的對象［13］。因此，它在地理地質(zhì)學(xué)、圖形與圖象處理、計(jì)算機(jī)動(dòng)畫仿真等許多領(lǐng)域都具有廣泛的應(yīng)用。由于分形插值函數(shù)具有很強(qiáng)的靈活性和穩(wěn)定性［14］，只要適當(dāng)調(diào)整其參數(shù)，所生產(chǎn)的插值曲線的維數(shù)就能取到1 和2 之間的任意值，因此它不僅可用來擬合光滑曲線和平穩(wěn)數(shù)據(jù)，更是在不光滑曲線和非平穩(wěn)數(shù)據(jù)的擬合中顯示出獨(dú)特的優(yōu)越性。河流動(dòng)力學(xué)領(lǐng)域，如劉德平［15］運(yùn)用分形理論和分形維數(shù)研究了長江漢口站4 個(gè)大洪水年份逐日平均水位、流量過程線的形態(tài)特征;馬宗偉等［16］分析了長江中下游河流的分析特征，說明了河流形態(tài)的分維大小與洪水發(fā)生的可能性高低之間的關(guān)系;Nikora［17］和鐘亮等［18］把分形理論應(yīng)用到床面粗糙形態(tài)的研究。

長江口水域的主要?jiǎng)恿σ蛩厥菑搅?、潮流和波浪。徑流受到三峽樞紐的調(diào)控，年內(nèi)雖有洪枯之分，但是年際間變化不大;外海潮汐呈周期性變化;波浪受季風(fēng)和臺(tái)風(fēng)影響，均具有周期性。由2011年至2013年長江口北槽12.5 m 深水航道月平均實(shí)測回淤量(圖1)可知，北槽航道沿程的回淤量具有自相似性(相似于圖1 中粗黑實(shí)線)，每月均有數(shù)日回淤量達(dá)到峰值和每年均有數(shù)月回淤量達(dá)到峰值。因此，可以采用分形理論進(jìn)行分析。鑒于現(xiàn)階段長江口深水航道回淤量預(yù)測涉及泥沙來源、動(dòng)力環(huán)境對泥沙回淤的影響機(jī)理等暫時(shí)無法解決的問題，本文嘗試?yán)梅中尾逯捣椒▽﹂L江口12.5 m 深水航道年回淤變化規(guī)律進(jìn)行探討，試圖揭示綜合因素對深水航道年回淤量的影響。

圖1 2011年至2013年長江口深水航道月平均回淤量Fig.1 Monthly backsilting quantity in Yangtze Estuary deepwater channel from 2011 to 2013

1 分形插值理論

傳統(tǒng)的數(shù)學(xué)插值擬合函數(shù)大都是多項(xiàng)式、有理函數(shù)或三角函數(shù)等初等函數(shù)性組合來表達(dá)，而分形插值函數(shù)則是用迭代函數(shù)系統(tǒng)(IFS)來實(shí)現(xiàn)的。經(jīng)典插值方法是針對歐氏幾何描述的對象，而山脈的輪廓、云彩的邊緣、粗糙表面的輪廓等圖像不能用歐氏函數(shù)很好地描述，用分形插值函數(shù)就能充分接近這些圖像，并且在適當(dāng)?shù)某叨确秶鷥?nèi)，它的分形維數(shù)趨近于這些數(shù)據(jù)的分形維數(shù)。

分形插值是構(gòu)造出一個(gè)迭代函數(shù)系，使目標(biāo)點(diǎn)集A 為插值函數(shù)的一部分。構(gòu)造目標(biāo)點(diǎn)集A 等于插值函數(shù)f(x)的迭代函數(shù)系{R2:Wi，i = 1，2，…，N}，方法如下:設(shè)迭代函數(shù)系中的每個(gè)函數(shù)Wi是仿射變換，其構(gòu)造表示為:

并且滿足條件:

根據(jù)式(1)和式(2)，對每一個(gè)i，仿射變換函數(shù)Wi的5 個(gè)常數(shù)ai、ci、di、ei與fi應(yīng)該滿足4 個(gè)線性方程式:

選用早熟、優(yōu)質(zhì)、高產(chǎn)，抗病蟲、抗逆性強(qiáng)、適應(yīng)性廣、商品性好的品種。當(dāng)前攀西地區(qū)種植較為普遍的品種為豫選一號品種，畝用種量2.5—3.5公斤。菜豆品種有紅花青殼、雙青玉豆、天馬架豆王、泰國架豆王、碧豐8號等。

四個(gè)方程有五個(gè)變量，所以存在一個(gè)自由變量，一般選di為自由變量，并稱之為垂直尺度因子。限定這個(gè)自由變量滿足－1 ＜di＜1(i = 1，2，…，N)。特別，如果取di= 0 ，就可以得到分段線性插值函數(shù)。當(dāng)di是任意限定的實(shí)數(shù)時(shí)，迭代函數(shù)系中的每個(gè)函數(shù)Wi的其他系數(shù)利用式(3)可以推導(dǎo):

2 時(shí)間尺度上回淤量的分形插值曲線

圖2 為長江口12.5 m 深水航道疏浚單元分布，2011年至2013年逐月航道回淤量見表1(其中負(fù)值表示回淤量，正值表示沖刷量)。采用分形插值方法，取垂直尺度因子－0.5≤di≤0.5(i=1，2，…，N)構(gòu)造長江口深水航道逐月的回淤量迭代函數(shù)系，對其每月回淤量進(jìn)行分形插值計(jì)算。對于航道回淤量而言，垂直尺度因子是一個(gè)涉及動(dòng)力(徑流、潮汐和風(fēng)浪)及泥沙特性(粒徑、沉速、粘性)的綜合影響因子，取不同的垂直尺度因子，可以得到不同的航道回淤量插值曲線，再對時(shí)間進(jìn)行積分，得到不同的年航道回淤量。所取di絕對值越大，則偏離實(shí)測回淤量越大。

根據(jù)式(1)～(7)，取di為－0.5、－0.4、－0.3、－0.2、－0.1、0、0.1、0.2、0.3、0.4 和0.5 對表1 中12 個(gè)月的航道回淤量進(jìn)行分形曲線插值。由式(1)～(7)可知，第n 次插值后數(shù)據(jù)點(diǎn)個(gè)數(shù)為N(N－1)n，并且與實(shí)測數(shù)據(jù)相比，插值后的曲線形成(N－1)n個(gè)與實(shí)測數(shù)據(jù)曲線完全相似的分形微元。圖3(a)為2011年長江口深水航道月回淤量的折線圖和月平均線圖。圖3(b)、(c)是根據(jù)2011年長江口深水航道每月實(shí)測回淤量，取不同的垂直尺度因子作三次插值后的曲線。圖3(b)曲線有156 個(gè)數(shù)據(jù)點(diǎn)，12 個(gè)分形微元;圖3(c)有1 872 個(gè)數(shù)據(jù)點(diǎn)，144 個(gè)分形微元;圖3(d)有22 464 個(gè)數(shù)據(jù)點(diǎn)，1 728 個(gè)分形微元。根據(jù)分形插值理論，在－1 ＜di＜1 內(nèi)經(jīng)過若干次插值后趨近于穩(wěn)定。由于分形插值方法收斂速度快，所以文中假設(shè)第三次插值后曲線基本趨于穩(wěn)定。同樣對2012年和2013年長江口深水航道每月回淤量進(jìn)行三次分形插值。從三年的實(shí)測回淤量來看，7月、8月和9月泥沙回淤量出現(xiàn)峰值，均比其他月份大數(shù)倍。因此，當(dāng)di＞0 時(shí)隨著di的增大，分形微元的峰沿Y 軸負(fù)方向增大;當(dāng)di＜0 時(shí)隨著di的減小分形微元的峰沿Y 軸正方向增大。第三次插值后所得回淤量插值數(shù)據(jù)點(diǎn)已非常龐大，插值曲線亦非常密集且基本趨于穩(wěn)定。

圖2 長江口12.5 m 深水航道示意Fig.2 The sketch of Yangtze Estuary 12.5 m deepwater channel

表1 2011年、2012年和2013年長江口深水航道每月回淤量Tab.1 Monthly backsilting quantity in Yangtze Estuary deepwater channel from 2011 to 2013

分形算法在做插值時(shí)，采用實(shí)測月份的淤積量作為控制節(jié)點(diǎn)，相鄰月份回淤的插值曲線呈現(xiàn)出控制點(diǎn)的“粗糙的、復(fù)雜的自相似”等特征，和實(shí)際的月航道淤積的時(shí)間序列分布特征并不一致卻又相似，即增加了月回淤量之間的分形維數(shù)。

圖3 2011年長江口深水航道月回淤量與分形插值曲線Fig.3 Monthly backsilting quantity in Yangtze Estuary deepwater channel and fractal interpolation curve in 2011

3 航道年回淤量計(jì)算和分析

3.1年回淤量計(jì)算公式構(gòu)建

插值出的回淤量在時(shí)間t 內(nèi)進(jìn)行積分即得到年回淤量V:

式中:Vi為第i 個(gè)月插值出的回淤量，ti為第i 個(gè)插值點(diǎn)X 坐標(biāo)，N 為初始實(shí)測數(shù)據(jù)點(diǎn)個(gè)數(shù)即為13，n 為分形插值的次數(shù)n=1，2，3。

式(8)計(jì)算的年回淤量即圖3(a)折線圖中折線與X 軸圍成的面積，而所有月回淤量的和即圖3(a)月平均線段與X 圍成的面積，兩種方法計(jì)算年回淤量略有差異。圖4 為2011年、2012年和2013年長江口深水航道年回淤量與di的變化關(guān)系圖，并對曲線進(jìn)行公式擬合。

2012年和2011年回淤量曲線相似，作一次插值所得到的直線的斜率相差較小，最小年回淤量為5 670 ×104m3，最大年回淤量為14 057 ×104m3。2013年作一次插值所得到的直線的斜率相差較大，長江口深水航道年回淤量最小值為3 860 ×104m3。

圖4 中實(shí)心點(diǎn)及曲線為作三次分形插值計(jì)算得出的年回淤量，空心點(diǎn)是所擬合公式的計(jì)算值，R 為相關(guān)系數(shù)(R1= R2= R3= 1 )。從圖4 中可知，不同垂直尺度因子di作一次插值得到的年回淤量為直線，作二次插值得到的年回淤量為二次曲線，作三次插值得到的年回淤量為三次曲線，且高次項(xiàng)的系數(shù)均等于作一次插值所得到的直線的斜率。

分形插值得到的年回淤量三次曲線和二次曲線均與作一次插值所得到的直線相切，且切點(diǎn)為垂直尺度因子為0 時(shí)實(shí)測回淤量根據(jù)式(8)計(jì)算得到的年回淤量，即Y 軸的截距。綜上所述，可以得出長江口12.5 m 深水航道年回淤量的計(jì)算公式:

式中:A 為作一次插值所得到的直線的斜率，B 為12 個(gè)月實(shí)測根據(jù)式(8)計(jì)算得到的年回淤量，n 為分形插值的次數(shù)n=1，2，3。式(9)中的系數(shù)A 與每年長江口深水航道月回淤量的曲線形狀有關(guān)系，實(shí)數(shù)B 是某年長江口深水航道的實(shí)測年回淤量。

圖4 長江口年回淤量及擬合公式Fig.4 Annual backsilting quantity in Yangtze Estuary channel and fitting formula

3.2 垂直尺度因子對航道年回淤量的影響

由于式(9)中系數(shù)A 和B 不確定因素較多，且直接受某年長江口深水航道實(shí)測回淤量的影響，所以對式(9)進(jìn)行如下處理:

式(10)中可以定義V年/B 為某年長江口深水航道相對回淤量。

圖5 相對年回淤量公式擬合Fig.5 The fitting formula of the relative annual backsilting quantity

圖5 為2011年、2012年和2013年相對回淤量各次插值數(shù)據(jù)擬合曲線。擬合曲線的相關(guān)系數(shù)均在0.95左右，于是可以得出長江口12.5 m 深水航道年相對回淤量計(jì)算公式:

式(11)中如果系數(shù)B 取多年平均值，則依據(jù)垂直尺度因子di可對長江口深水航道年回淤量進(jìn)行預(yù)測和評估。若n 取無窮大，且－1 ＜di＜1 式(11)可簡化為:

取B 為2011年、2012年和2013年深水航道年回淤量的平均值8 707.5 ×104m3。根據(jù)2011年、2012年和2013年深水航道年回淤量(V年)和式(12)，反求di依次為－0.026、0.177 和－0.104?？紤]年際間的影響，按照三年平均回淤量與最大和最小回淤量之比20%調(diào)整di的范圍為－0.12≤di≤0.20，則近些年長江口深水航道年回淤量的變化范圍:8 023.7 ×104m3≤V年≤10 303.2 ×104m3。當(dāng)取di= －0.5 時(shí)，其年最小回淤量為6 580 ×104m3;當(dāng)di= －1 時(shí)，得到年極限最小回淤量為5 516 ×104m3。

4 結(jié) 語

依據(jù)2011 至2013年長江口12.5 m 深水航道實(shí)測年回淤量資料，采用分形插值方法，取垂直尺度因子－0.5≤di≤0.5(i=1，2，…，N)對長江口深水航道回淤量預(yù)測進(jìn)行了初步研究。主要結(jié)論如下:

1)分形插值可以較好地反映長江口12.5 m 深水航道的月回淤量的變化態(tài)勢并通過調(diào)整垂直尺度因子，可以得到包含動(dòng)力和泥沙特性在內(nèi)的綜合因素對深水航道年回淤量的影響。

2)根據(jù)2011 至2013年12.5 m 深水航道年回淤量資料，確定了垂直尺度因子的變化范圍為－0.12≤di≤0.20。由此計(jì)算得到近些年的年回淤量的變化范圍為8 023.7 ×104m3≤V年≤10 303.2 ×104m3，其最小回淤量為6 580 ×104m3以及年極限最小回淤量為5 516 ×104m3。

3)由于目前長江口12.5 m 深水航道只有3年的回淤資料，因此該方法的有效性尚待更長時(shí)間資料的檢驗(yàn)。

［1］CHEN J Y，YU C X，XU H G.The model of the Changjiang Estuary during the last 2000 years，Estuarine comparison［M］.Massachusetts:Academic Press，1982:655-666.

［2］潘定安，孫介民.長江口攔門沙地區(qū)的泥沙運(yùn)動(dòng)規(guī)律［J］.海洋與湖沼，1996，27(2):279-286.(PAN Dingan，SUN Jiemin.The sediment dynamics in the Changjiang River Estuary mouth bar area［J］.Oceanologia et Limnologia Sinica，1996，27(2):279-286.(in Chinese))

［3］李九發(fā)，時(shí)偉榮，沈煥庭.長江河口最大渾濁帶的泥沙特性和輸移規(guī)律［J］.地理研究，1994，13(1):51-60.(LI Jiufa，SHI Weirong，SHEN Huanting.Sediment properties and transportation in the turbidity maximum in Changjiang Estuary［J］.Geographical Research，1994，13(1):51-60.(in Chinese))

［4］張志忠，阮文杰，蔣國俊.長江口動(dòng)水絮凝沉降與攔門沙淤積的關(guān)系［J］.海洋與湖沼，1995，26(6):662-669.(ZHANG Zhizhong，RUAN Wenjie，JIANG Guojun.The relationship between flocculating settlement in flowing water and deposition in the Changjiang River［J］.Oceanology et Limnology Sinica，1995，26(6):662-669.(in Chinese))

［5］時(shí)鐘，凌鴻烈.長江口細(xì)顆粒懸沙濃度垂向分布［J］.泥沙研究，1999(2):59-64.(SHI Zhong，LING Honglie.Vertical profiles of fine suspension concentration in the Changjiang Esturary［J］.Journal of Sediment Research，1999(2):59-64.(in Chinese))

［6］盛升國，施厚慶，于福生，等.長江口北槽航道拋泥區(qū)泥沙擴(kuò)散的試驗(yàn)研究［J］.泥沙研究，1986(3):37-47.(SHENG Shengguo，SHI Houqing，YU Fusheng，et al.Sedimentation committee of chinese society of hydraulic engineering accepts the consulting of sedimentation problems［J］.Journal of Sediment Research，1986(3):37-47.(in Chinese))

［7］徐建益，袁建忠.長江口深水航道建設(shè)中的浮泥研究及述評［J］.泥沙研究，2001(3):74-82.(XU Jianyi，YUAN Jianzhong.Study on the fluid mud in the Yangtze Estuary［J］.Journal of Sediment Research，2001(3):74-82.(in Chinese))

［8］DOU X P，LI T L，DOU G R.Numerical model of total sediment transport in the Yangtze estuary ［J］.China Ocean Engineering，1999，13(3):277-286.

［9］劉猛，吳華林，李為華，等.長江口深水航道工程南導(dǎo)堤越堤水沙運(yùn)動(dòng)觀測研究［J］.海洋工程，2013，31(3):129-140.(LIU Meng，WU Hualin，LI Weihua，et al.Observational study on flow and sediment transport overtopping the south leading jetty in the Yangtze Estuary Deepwater Channel Regulation Project［J］.The Ocean Engineering，2013，31(3):129-140.(in Chinese))

［10］陳維，匡翠萍，顧杰，等.南沙頭通道及橫沙通道對長江口深水航道的影響分析［J］.海洋科學(xué)，2013，37(4):75-81.(CHEN Wei，KUANG Cuiping，GU Jie，et al.Influences of the Nanshatou Passage and Hengsha Passage on sediment deposition in Deepwater Navigation Channel of the Changjiang River Estuary［J］.Marine Science，2013，37(4):75-81.(in Chinese))

［11］沈淇，顧峰峰，萬遠(yuǎn)揚(yáng)，等.長江口洪季北槽深水航道區(qū)域懸沙沉降速度估算［J］.海洋工程，2013，31(2):88-94.(SHEN Qi，GU Fengfeng，WAN Yuanyang，et al.Settling velocity of suspended sediment at the North Passage Deep-water Channel of the Yangtze Estuary in the flood season［J］.The Ocean Engineering，2013，31(2):88-94.(in Chinese))

［12］張宏，韓震.基于神經(jīng)網(wǎng)絡(luò)的長江口深水航道海域懸浮泥沙濃度研究［J］.海洋通報(bào)，2010，29(2):213-219.(ZHANG Hong，HAN Zhen.Research on suspended sediment concentration of the Yangtze River Estuary Deepwater Channel sea ares using neural network［J］.Marine Science Bulletin，2010，29(2):213-219.(in Chinese))

［13］BARNSLEY M F.Fractal functions and interpolation［J］.Constr.Approx.，1986，2:303－329.

［14］WANG H Y，YU J S.Fractal interpolation functions with variable parameters and their analytical properties［J］.J.Approx.Theory，2013，175:1－18.

［15］劉德平.分形理論在水文過程形態(tài)特征分析中的應(yīng)用［J］.水利學(xué)報(bào)，1998(2):20－25.(LIU Deping.The application of fractal theory in hydrologic process shape characteristic analysis［J］.Journal of Hydraulic Engineering，1998(2):20－25.(in Chinese))

［16］馬宗偉，許有鵬，李嘉峻.河流形態(tài)的分維與洪水關(guān)系的探討［J］.水科學(xué)進(jìn)展，2005，16(4):530－534.(MA Zongwei，XU Youpeng，LI Jiajun.River fractal dimension and the relationship between river fractal dimension and river flood［J］.Advance in Water Science，2005，16(4):530－534.(in Chinese))

［17］NIKORA V I.Fractal structures of river plan forms［J］.Water Resources Research，1991，27(6):1327－1333.

［18］鐘亮，許光祥.床面粗糙形態(tài)的二元分形插值模型［J］.水科學(xué)進(jìn)展，2011，22(5):662－668.(ZHONG Liang，XU Guangxiang.Bivariate fractal interpolation for estimating rough channel bedform［J］.Advance in Water Science，2011，22(5):662－668.(in Chinese))