趙江艷， 劉 文， 張明陽， 徐強(qiáng)強(qiáng)， 謝 平,4

(1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072； 2. 武漢理工大學(xué) 內(nèi)河航運技術(shù)湖北省重點實驗室，武漢 430063； 3.國家水運安全工程技術(shù)研究中心，武漢 430063； 4. 國家領(lǐng)土主權(quán)與海洋權(quán)益協(xié)同創(chuàng)新中心，武漢 430072)

長江屬于季節(jié)性大型河流，受季節(jié)性因素和通航樞紐下泄流量等因素影響，航道水位變化明顯，尤其是長江中游航道，年最高水位與年最低水位的多年平均值之差高達(dá)13 m，對船舶通航安全的影響顯著。[1]通過對長江航道的事故進(jìn)行統(tǒng)計分析發(fā)現(xiàn)，2013年和2014年枯水期發(fā)生的水上事故均占當(dāng)年事故總數(shù)的50%以上，其中長江中游航道是長江干線的礙航瓶頸，有近20處礙航淺灘及10多個重點淺水道，通航條件十分不穩(wěn)定，極易發(fā)生擱淺事故，從而造成礙航局面。[2]因此，對水位變化下的長江中游船舶擱淺風(fēng)險進(jìn)行分析具有實際應(yīng)用價值。

船舶擱淺風(fēng)險分析已在水上交通工程領(lǐng)域引起廣泛關(guān)注，但不同通航環(huán)境下的船舶擱淺風(fēng)險的定量計算依然是亟待解決的問題。通常采用的研究方法有故障樹方法、貝葉斯方法及概率模型等。[2-8]AMROZOWIEZ[6]將船舶擱淺事故分為動力性擱淺和漂移性擱淺，分別建立故障樹，確定失效概率；KITE-POWELL等[7]根據(jù)統(tǒng)計學(xué)中的貝葉斯方法提出船舶擱淺的風(fēng)險計算模型，可有效預(yù)測船舶擱淺風(fēng)險，定量計算相應(yīng)的直接經(jīng)濟(jì)損失和間接經(jīng)濟(jì)損失；ARSHAM等[8]利用船舶交通和航道復(fù)雜度建立評估船舶擱淺概率的模型等。從以上研究成果中可看出，大多數(shù)研究是從整體的角度或從船舶的擱淺現(xiàn)象出發(fā)，闡述有關(guān)船舶擱淺分析的總體看法。但是，船舶擱淺分析中的影響因素與研究區(qū)域的交通條件、自然條件緊密相關(guān)，因此無法推廣到所有地區(qū)。此外，由于交通環(huán)境和自然環(huán)境會隨時間變化，因此船舶的擱淺分析應(yīng)是一個適應(yīng)新的風(fēng)險因子的動態(tài)過程。對此，本文將船舶擱淺問題視為首次穿越問題，根據(jù)隨機(jī)振動的首次穿越失效理論建立首次穿越失效(船舶擱淺)模型。首次穿越問題是指一個受隨機(jī)參激或隨機(jī)外激作用的動力學(xué)系統(tǒng)，其響應(yīng)在相空間上發(fā)生隨機(jī)擴(kuò)散。當(dāng)系統(tǒng)響應(yīng)第一次離開某個安全域或允許域時，就稱系統(tǒng)發(fā)生首次穿越。首次穿越在工程[9]、物理[10]、生物[11]、醫(yī)學(xué)[12]及經(jīng)濟(jì)[13]等學(xué)科中廣泛存在。QUY等[14]將首次穿越理論引入到船舶擱淺的概率模型中，初步驗證該理論在解決船舶擱淺問題上的適用性，但在具體的建模過程中僅考慮航速變化對船舶吃水的影響，沒有解決航道水位變化對船舶擱淺風(fēng)險影響的問題。

對此，本文重點研究長江中游航道水位變化對船舶通航安全的影響。為保證分析結(jié)果的可靠性，運用隨機(jī)振動的首次穿越失效理論構(gòu)建船舶擱淺概率計算模型，定量分析不同水位下船舶的擱淺概率；同時，給出船舶擱淺概率分別為0和1時對應(yīng)的水位臨界值，用以驗證實際的設(shè)計最低通航水位是否符合通航安全的要求。

1 問題描述與建模

船舶擱淺事故是一種受多種因素影響的水上交通事故，在分析船舶擱淺風(fēng)險的影響要素時，除了不同時期的水位變化之外，人為因素、船舶因素等不受水位變化影響的因子也會對船舶的擱淺風(fēng)險造成一定的影響。若進(jìn)行全面考慮，則會增加求解問題的難度，甚至無法求解。因此，本文在船舶擱淺風(fēng)險的分析中主要考慮水位變化對船舶擱淺的影響，參考波浪載荷下船舶運動響應(yīng)研究中的譜分析方法的基本思路分析船舶的垂向運動，基于首次穿越理論構(gòu)建船舶擱淺風(fēng)險計算模型。

1.1 水位變化對船舶擱淺的影響

擱淺事故表明，航道水深不足是導(dǎo)致船舶發(fā)生擱淺的直接原因。[15]由于航道深度在一定時期內(nèi)是不變的，因此航道水位的變化也可間接描述航道水深的變化，從而建立水位與船舶擱淺風(fēng)險之間的關(guān)系。富余水深、航道水深及船舶吃水之間的關(guān)系見圖1，其中：H為航道水深；W為當(dāng)時當(dāng)?shù)貙嶋H水位；d為圖示水深；D為船舶吃水；Δh為富余水深。富余水深主要包括船舶航行下沉量和觸底安全富余量2部分，其主要作用是防止船舶觸底，保證船舶安全航行。有關(guān)系式

H=W+d

(1)

允許船舶的最大吃水

Dmax=H-Δh=W+d-Δh

(2)

一般而言，增大富余水深Δh的取值會提高船舶的通航安全，同時會降低航道的通航能力。因此，為平衡通航安全與通航能力之間的關(guān)系，須對富余水深Δh進(jìn)行合理取值。《內(nèi)河通航標(biāo)準(zhǔn)》[16]根據(jù)航道等級規(guī)定了相應(yīng)的富余水深。

由式(2)可計算任意水位下允許船舶的最大吃水，當(dāng)實際吃水大于該最大吃水時，船舶發(fā)生擱淺。由生產(chǎn)實踐可知，即使船舶的載貨量一定，船舶的吃水也會隨著通航環(huán)境的變化而變化。這是因為船舶航行時不僅會受波浪的影響產(chǎn)生升沉運動的響應(yīng)位移，還會受周圍水流的影響出現(xiàn)下沉現(xiàn)象。船舶航行下沉量(即船舶動吃水)是指船舶航行中的吃水較靜止時的增加量，本文采用實際航運中推薦的公式，即

(3)

式(3)中：CB為方形系數(shù)；v為船速，kn。

因此，計算船舶擱淺概率的前提是確定當(dāng)時當(dāng)?shù)氐乃患按暗膭討B(tài)吃水。對于載貨量一定的單艘次船舶來說，水位越高，船舶的擱淺風(fēng)險越小。

1.2 水位變化下船舶擱淺風(fēng)險計算模型

1.2.1建模的基本假設(shè)

為避免模型自身的復(fù)雜性帶來不確定性風(fēng)險，參考波浪載荷下船舶運動響應(yīng)研究的譜分析方法[17-20]，確定建模的2個基本假設(shè)：

(1)波浪表面是一個各態(tài)歷經(jīng)的零均值高斯隨機(jī)過程；

(2)船舶為一個線性系統(tǒng)。

假設(shè)(1)與隨機(jī)波浪理論的思想一致，認(rèn)為在任意較短的時間內(nèi)，波浪符合零均值的平穩(wěn)正態(tài)隨機(jī)過程，其統(tǒng)計特性[21-23]可用功率譜密度來確定，即

(4)

式(4)中：S(ωi)為波浪功率譜密度函數(shù)，表示一個隨機(jī)波浪所包含的所有頻率及其對應(yīng)的能量。

在波浪載荷的作用下，船舶會產(chǎn)生位移、切力、彎矩及應(yīng)力等各種響應(yīng)。[24]若將波浪的作用記為Q(t)，則有

(5)

式(5)中：mQ為Q(t)的均值；An和Bn為互不相關(guān)的隨機(jī)變量。An和Bn有均值為零、方差相等的特征，即mAn=mBn=0，DAn=DBn=Dn。

波浪的方差為

(6)

(7)

因此，當(dāng)Δω→0時，可得

(8)

式(7)和式(8)中：SQ(ω)為波浪譜。

對于單自由度系統(tǒng)，假定受一簡諧干擾力作用，則振動方程為

(9)

式(9)中：q(t)為廣義位移；ω0為系統(tǒng)固有頻率；ξ為黏性阻尼系數(shù)；Q0為干擾力幅值；ω為干擾力頻率。

該系統(tǒng)的響應(yīng)為

q(t)=q0eiωt=H(ω)Q0eiωt=H(ω)Q(t)

(10)

根據(jù)假設(shè)(2)，將垂向運動的船舶視為單自由度的線性系統(tǒng)，將波浪作為系統(tǒng)的輸入量，則船舶垂向運動的響應(yīng)即為系統(tǒng)的一種輸出，其值由波浪載荷和響應(yīng)函數(shù)共同決定。響應(yīng)函數(shù)的物理意義是各種不同頻率的單位波幅余弦波產(chǎn)生的船舶響應(yīng)。若定義H(ω)的模|H(ω)|為響應(yīng)函數(shù)，則有

|H(ω)|=q/Q0

(11)

利用傅里葉變換可知

(12)

因此，作為船體—波浪線性系統(tǒng)的隨機(jī)輸出，波浪載荷在短時間內(nèi)也服從零均值正態(tài)分布，其幅值遵循Rayleigh分布特性。[25-26]

1.2.2船舶垂向運動建模

建立坐標(biāo)系O-xyz見圖2，當(dāng)將船體作為剛體看待時，船舶的搖蕩運動是沿3個坐標(biāo)軸方向的位移和轉(zhuǎn)動，具有6個自由度。在不考慮舵力作用的情況下，船舶的六自由度運動可分解為縱向和橫向2組耦合運動，其中：縱向運動由升沉、縱搖和縱蕩組成；橫向運動由橫蕩、橫搖和艏搖組成。

縱向耦合運動與橫向耦合運動之間的影響很小，且自由狀態(tài)下船舶的縱蕩運動對船舶升沉運動位移的分析影響不大。因此，在求解船舶在波浪中的升沉運動位移時，不用考慮船舶的縱蕩運動，只需研究船舶升沉和縱搖的耦合運動。

船舶升沉和縱搖的耦合運動的微分方程為

(13)

(14)

(15)

式(15)中：

水動力系數(shù)可通過計算近似得到。由此可得

(16)

式(16)中：Az=E-1M；Bz=E-1N；隨機(jī)波浪W為有色噪聲。

若選擇船舶升沉運動的位移和縱搖傾角為測量狀態(tài)，則測量方程為

Y=CzX+V

(17)

本文利用波浪擾動成形濾波器計算船舶受到的波浪干擾，其核心思想是將高斯分布的白噪聲通過一個波浪擾動成形濾波器G(s)，從而得到作為隨機(jī)波浪擾動仿真過程的有色噪聲。[17]若將成形濾波器視為船舶運動系統(tǒng)的一部分，則可將白噪聲信號視為新系統(tǒng)的擾動。新系統(tǒng)的組成見圖3。

圖3中：GZ(s)為波浪對船舶縱向擾動力的成形濾波器；GM(s)為波浪對船舶縱向擾動力矩的成形濾波器。通常選取成形濾波器

(18)

式(18)中：s為Laplace算子；G0為增益；ω0為中心頻率；ξ為阻尼系數(shù)。wZ(t)和wM(t)均為白噪聲信號，其均值為0，方差為1；Z(t)為船舶升沉運動干擾力；M(t)為縱搖干擾力矩；q(t)為最終求得的船舶垂向響應(yīng)位移。

因此，當(dāng)波浪的白噪聲信號wZ(t)和wM(t)已知時，可通過上述成形濾波器得到船舶升沉運動干擾力Z(t)和縱搖干擾力矩M(t)，隨后結(jié)合船舶運動的狀態(tài)方程求出船舶垂向運動響應(yīng)位移q(t)。

可運用Simulink工具箱對上述新系統(tǒng)進(jìn)行仿真，從而得到船舶垂向運動響應(yīng)位移q(t)。仿真流程圖見圖4。

1.2.3船舶擱淺風(fēng)險量化建模

由“1.2.1”節(jié)可知，船舶在波浪激勵下產(chǎn)生響應(yīng)，由于激勵是服從高斯分布的，運動響應(yīng)是線性的，最后得到的船舶運動也符合高斯分布。因此，若運用船舶擱淺的數(shù)學(xué)關(guān)系式D>H，即船舶吃水大于航道水深，則船舶擱淺問題可看成首次穿越失效問題：將航道水深表示成隨時間變化的閾值，當(dāng)船舶垂向的位移響應(yīng)首次穿越該閾值時，認(rèn)為船舶發(fā)生擱淺。首次穿越概率(即船舶的擱淺概率)可通過相關(guān)的隨機(jī)理論計算得到。

通過分析可知，船舶在時域內(nèi)穿越的失效閾值可表示為

φ(t)=W(t)+h-D

(19)

式(19)中：W為當(dāng)時當(dāng)?shù)貙嶋H水位；h為航道圖示水深；D為船舶吃水。

圖5為船舶擱淺(首次穿越失效)示意，若在船舶航行期間(即[0，T]內(nèi))，其升沉運動響應(yīng)位移q(t)穿過φ(t)，則船舶首次穿越失敗，即發(fā)生擱淺事故。

令β=φ(t)，則船舶響應(yīng)的期望穿閾率[27]可表示為

(20)

式(20)中：m0和m2分別為船舶運動響應(yīng)譜的零階矩及二階矩。m0和m2的計算式可表示為

(21)

(22)

式(21)和式(22)中：Sq(ω)為描述在頻率域ω上響應(yīng)x(t)的均值分布的功率譜密度函數(shù)，即船舶運動響應(yīng)譜。船舶運動的響應(yīng)譜實質(zhì)上是一個功率譜密度函數(shù)，描述信號功率隨頻率的分布狀況。

根據(jù)“1.2.2”節(jié)可求出船舶垂向運動的響應(yīng)位移q(t)。將q(t)視為一個隨時間變化的隨機(jī)信號。以Ts對隨機(jī)信號進(jìn)行等間隔抽樣，得到隨機(jī)序列q(nTs)，記為X(n)；隨后通過直接法、間接法、Bartlett法或Welch法計算得到其譜密度函數(shù)。

若應(yīng)用首次穿越理論，則船舶首次穿越失效的概率(即船舶的擱淺概率)可表示為

P(β|T)=1-exp(-vβT)

(23)

2 實例計算與結(jié)果分析

選取一艘滿載排水量為2 858 t的船舶作為研究對象，運用本文提出的方法計算其擱淺風(fēng)險。該船的具體參數(shù)為：船長75 m；垂線間長71.33 m；型寬13.4 m；型深5.4 m；設(shè)計吃水4.2 m；方形系數(shù)0.704；水線面系數(shù)0.980；中橫剖面系數(shù)0.980；滿載排水量2 858 t。

根據(jù)水文和氣象條件，將內(nèi)河船舶航行區(qū)域劃分為A，B，C等3級(見表1)。長江中游主要屬于B級航區(qū)，根據(jù)該航區(qū)的波浪要素及船舶尺度確定

表1 我國內(nèi)河各級航區(qū)波浪要素

波浪成形濾波器的參數(shù)見表2。

表2 波浪成形濾波器的參數(shù)

將上述船舶尺度、圖示水深(3.2 m)及船舶與波浪的遭遇頻率(π/5)等數(shù)據(jù)作為輸入，結(jié)合MATLAB/Simulink仿真平臺(如圖4所示)，便可求出波浪干擾力、波浪干擾力矩及船舶升沉運動的響應(yīng)位移。設(shè)定仿真時間為60 s，波浪對船舶升沉運動的干擾力和干擾力矩分別見圖6及圖7。從圖6和圖7中可看出：船舶受到的波浪載荷符合正態(tài)分布，其幅值符合Rayleigh分布。因此，當(dāng)將船舶視為一個線性系統(tǒng)時，船舶的運動響應(yīng)過程是各態(tài)歷經(jīng)的平穩(wěn)隨機(jī)過程，運動幅值也符合Rayleigh分布(見圖8)。

運用基于MATLAB的譜估計程序?qū)崿F(xiàn)船舶升沉運動響應(yīng)位移到船舶響應(yīng)譜的轉(zhuǎn)化，最終得到船舶升沉運動的響應(yīng)譜見圖9。隨后運用首次穿越理論計算得到船舶擱淺概率與水位的關(guān)系見圖10。

從圖10中可看出：隨著航道水位下降，船舶擱

淺風(fēng)險逐漸增大；隨著航道水位升高，船舶擱淺風(fēng)險逐漸減小。船速越大，船舶擱淺概率為0對應(yīng)的水位臨界值越大；船速越小，船舶擱淺概率為0對應(yīng)的水位臨界值越小。具體臨界值見表3。

表3 不同船速v下的水位臨界值

通過分析可知：船速越大，船舶的運動響應(yīng)越顯著，船舶升沉運動的響應(yīng)幅值越大；同時，船速越大，船舶航行下沉量越大，故其對通航水位的要求越高。因此，利用上述方法不僅可確定某一水位下船舶的擱淺概率，還能由船舶擱淺風(fēng)險的容忍度確定航道的最低通航水位，從而驗證實際的設(shè)計最低通航水位是否符合要求。

此外，由于當(dāng)時當(dāng)?shù)氐膶嶋H水位與圖示水深之和即為航道水深，故該方法也可為航道部門確定航道維護(hù)水深提供一定的參考。

3 結(jié)束語

本文提出一種基于首次穿越失效理論的船舶擱淺風(fēng)險的定量計算方法，建立水位變化下的首次穿越失效(船舶擱淺)模型。同時，以長江中游航道的典型船舶為例，定量計算其在淺水航段的擱淺風(fēng)險，最終得到船舶擱淺概率與航道水位的對應(yīng)關(guān)系。從結(jié)果上看：在同一船速下，航道水位越低，船舶的擱淺風(fēng)險越大；船速越大，船舶擱淺概率為0和1對應(yīng)的水位臨界值越大，與實際情況一致。研究結(jié)果表明，該方法不僅能確定某一水位下船舶的擱淺概率，還可確定船舶擱淺概率為0和1時對應(yīng)的水位臨界值。此外，利用該方法還能求得任一可接受擱淺概率下的最低水位值。因此，該方法可用于確定船舶擱淺風(fēng)險容忍度下的最低通航水位，為最低通航水位的設(shè)計提供新的思路，為航道部門確定航道維護(hù)水深和航運經(jīng)營人對船舶進(jìn)行配載提供參考。

[1] 劉婷.水位變化對長江中游航運物流成本的影響機(jī)理及模型研究[D]. 武漢: 武漢理工大學(xué), 2014.

[2] 張笛. 枯水期長江通航風(fēng)險評價和預(yù)測方法研究[D]. 武漢: 武漢理工大學(xué), 2011.

[3] 陳婷婷, 施朝健, 雷琴. 船舶擱淺模糊故障樹的灰色模糊關(guān)聯(lián)分析[J]. 中國航海, 2015, 38(3):83-88.

[4] 吳金龍, 劉大剛, 范中洲,等. 大風(fēng)浪條件下渤海海區(qū)重點船舶風(fēng)險評估模型[J]. 中國航海, 2012, 35(2):89-92.

[5] MAZAHERI A, MONTEWKA J, KUJALA P. Modeling the Risk of Ship Grounding—A Literature Review from a Risk Management Perspective[J]. WMU Journal of Maritime Affairs, 2014, 13(2): 269-297.

[6] AMROZOWIEZ M D. The Need for a Probabilistic Risk Assessment of the Oil Tanker Industry and a Qualitative Assessment of Oil Tanker Groundings[Z]. Massachusetts Institute of Technology, 1996.

[7] KITE-POWELL H L, JIN D L, JEBSEN J L, et al. Investigation of Potential Risk Factors For Groundings of Commercial Vessels In U S Ports[J]. International Journal of Offshore & Polar Engineering, 1999, 9(1)：16-21.

[8] ARSHAM M, JAKUB M, PENTTI K, et al. Assessing Grounding Frequency Using Ship Traffic and Waterway Complexity[J]. Journal of Navigation, 2015, 68(1):89-106.

[9] ZHANG L Y, SUN L M, SHANG Z Q. Real-Time Reliability Assessment Based on Acceleration Monitoring for Bridge[J]. Science China Technological Sciences, 2016,56(8): 1-11.

[10] SHARMA K, KUMAR N. First-Passage Time: Lattice Versus Continuum[J]. Physical Review E, 2012, 86(3): 032104.

[11] HU Z, CHENG L, BERNE B J. First-Passage Time Distribution in Stochastic Processes with Moving and Static Absorbing Boundaries with Application to Biological Rupture Experiments[J]. Journal of Chemical Physics, 2010, 133(3): 034105.

[12] SIROVICH L, KNIGHT B. Spiking Neurons and the First Passage Problem[J]. Neural Computation, 2011, 23(7): 1675-1703.

[13] BAKSHI G, PANAYOTOV G. First-Passage Probability, Jump Models, and Intra-Horizon Risk[J]. Journal of Financial Economics, 2010, 95(1): 20-40.

[14] QUY N M, VRIJLING J K, GELDER P, et al. On the Assessment of Ship Grounding Risk in Restricted Channels [C]. 8th International Conference on Marine Sciences and Technologies, 2006: 294-299.

[15] 任亞磊, 牟軍敏, 李亞軍, 等. 進(jìn)港航道中船舶擱淺概率的蒙特卡洛仿真[J]. 船舶力學(xué), 2014, 18(5): 532-539.

[16] 交通運輸部.內(nèi)河通航標(biāo)準(zhǔn)：GB 50139—2014[S]. 北京: 中國計劃出版社, 2014.

[17] 彭秀艷, 張文穎, 何志勇. 基于譜估計方法的船舶縱向運動建模預(yù)報研究[J]. 自動化技術(shù)與應(yīng)用, 2008, 27(11): 1-4.

[18] DINIS M S, PIERSON J W J. On the Motion of Ships in Confused Seas[Z]. Routledge, 1953,22(23):8144-8154.

[19] 李成成. 船舶在波浪中運動時域數(shù)值模擬研究[D]. 上海：上海交通大學(xué), 2010.

[20] 杭棟棟. 基于船舶運動響應(yīng)數(shù)據(jù)的海浪估計方法研究[D]. 哈爾濱：哈爾濱工程大學(xué), 2013.

[21] 徐德倫, 于定勇. 隨機(jī)海浪理論[M]. 北京：高等教育出版社, 2001:111-115.

[22] 俞聿修. 隨機(jī)波浪及其工程應(yīng)用[M]. 大連: 大連理工大學(xué)出版社, 2011:60-81.

[23] KIMURA A. Statistical Properties of Random Wave Groups[Z]Coastal Engineering,1980:2955-2973.

[24] 王學(xué)亮. 大型起重船在海浪中的運動響應(yīng)研究[D]. 天津: 天津大學(xué), 2003.

[25] 吳小平. 基于切片理論的波浪載荷直接計算[J]. 船舶與海洋工程, 2010(4):21-25.

[26] PEREZ T, BLANKE M. Ship Roll Damping Control[J]. Annual Reviews in Control, 2012, 36(1):129-147.

[27] QUY N M, VRIJLING J K, VAN GELDER P. Risk-and Simulation-Based Optimization of Channel Depths: Entrance Channel of Cam Pha Coal Port [J]. Simulation, 2008, 84(1): 41-45.