長(zhǎng)江口疏浚土利用的數(shù)理優(yōu)化

周 瑜，鄭偉安

（華東師范大學(xué)國際航運(yùn)物流研究院，上海 200062）

1 前言

自19世紀(jì)40年代《南京條約》簽訂，以上海為代表的第一批通商口岸被迫開放以來，上海作為我國對(duì)外貿(mào)易往來的重要交通樞紐為眾人所熟知。特別是19世紀(jì)70年代，伴隨著黃浦江和蘇州河兩岸近代工業(yè)集聚區(qū)的興起，上海逐步奠定了其東方航運(yùn)中心的地位。新中國成立后，借助于上海市及長(zhǎng)江三角洲的緊密聯(lián)系，通過臨港工業(yè)的互動(dòng)、國際航運(yùn)中心建設(shè)的持續(xù)推動(dòng)，上海港保持著上升的發(fā)展勢(shì)頭，維持著其強(qiáng)有力的競(jìng)爭(zhēng)優(yōu)勢(shì)。

然而，在發(fā)展過程中，上海港及臨港產(chǎn)業(yè)正遭遇發(fā)展空間不足的瓶頸。

首先，在有限的地理維度中，上海和長(zhǎng)江口深水岸線及土地已經(jīng)達(dá)到了近乎飽和的狀態(tài)。具體而言，目前上海港主力港區(qū)分布于外高橋、洋山、羅涇，而臨港產(chǎn)業(yè)主要集中在江口南岸、杭州灣金山、臨港新城及長(zhǎng)江口長(zhǎng)興島南岸。因而，重新尋找到完整的深水岸線及近岸土地已然成為制約發(fā)展的難題。

其次，上海港集裝箱碼頭和散雜貨碼頭的能力已趨超負(fù)荷。以外高橋港區(qū)為例，2011年外高橋港區(qū)集裝箱吞吐量已是設(shè)計(jì)能力的1.39倍，而洋山港區(qū)吞吐量也已達(dá)到了設(shè)計(jì)能力的1.41倍。有研究推算，上海港需要維持每年2×106～3×106TEU增量才能達(dá)到長(zhǎng)江流域?qū)ν赓Q(mào)易發(fā)展的需求，而臨港臨?？捎霉I(yè)岸線土地資源已無法滿足。

再者，上海港作為建設(shè)中的國際航運(yùn)中心，功能尚未完善，尤其水水中轉(zhuǎn)能力較弱。目前部分集裝箱的水水中轉(zhuǎn)功能主要由洋山港兼顧，大宗散貨的中轉(zhuǎn)均有羅涇港區(qū)和外高橋港區(qū)兼顧，但受各港區(qū)自身的條件限制，中轉(zhuǎn)能力較弱（特別是與長(zhǎng)江流域各港口），無法使上海港成為完善的綜合運(yùn)輸樞紐。

最后，上海周邊大型港口不斷崛起，上海國際航運(yùn)中心的龍頭地位受到嚴(yán)峻挑戰(zhàn)。例如與鄰近的舟山港相比較，2011年，上海港包括內(nèi)河在內(nèi)的全港貨物吞吐量達(dá)到7.27×108t，而寧波—舟山港為6.94×108t，當(dāng)年兩港增速分別為11.4%和12%，從增速看，勢(shì)均力敵。但2012年，寧波—舟山港增速快速領(lǐng)先上海港，達(dá)到7.2%，遠(yuǎn)超上海港的1.1%，同時(shí)寧波—舟山港貨物吞吐量達(dá)到7.44×108t，首次突破7×108t，超過上海港800余萬噸，從而將上海港從雄踞多年的貨物吞吐量全球第一的位置上挑落下馬。從集裝箱增量看，2012年，上海港比2011年增加7.9×105TEU，而寧波—舟山增量為9.8×105TEU，此項(xiàng)上也是首次超過上海港。綜上所述，上海港面臨的問題集中在上海土地、岸線、航道等資源后備嚴(yán)重不足，制約了上海港及臨港產(chǎn)業(yè)發(fā)展空間及布局。在此港口資源發(fā)展壁壘的情況下，橫沙東灘的開發(fā)可謂是應(yīng)運(yùn)而生。

第一，橫沙東灘近期可提供的土地為1.7×105畝（1畝≈666.67m2大一區(qū)范圍），全部成陸后總共成陸的面積約5.8×105畝，可形成的岸線總長(zhǎng)在90 km以上，并且可利用直接面向外海的有利條件，采用大型挖入式港池等手段，增加深水岸線，解決了上海港所無法滿足的土地和深水岸線資源的問題。

第二，橫沙東灘自身擁有非常好的優(yōu)勢(shì)，作為長(zhǎng)江出海口的橋頭堡，可滿足江海聯(lián)運(yùn)的緊迫要求。一方面，橫沙東灘通江達(dá)海，具有良好的地理位置。它位于長(zhǎng)江出?？?，扼守我國海岸線與長(zhǎng)江黃金水道的T字形交點(diǎn)，若與長(zhǎng)興島（海洋裝備島）用短距離隧道或橋梁連通后，即可經(jīng)滬崇蘇陸上通道直抵上海浦東和蘇北。此外，橫沙東灘與洋山深水港水域距離約100多千米，與外高橋港區(qū)（南港南岸）水域距離約30千米。借助其地理優(yōu)勢(shì)，可逐漸形成上海國際航運(yùn)中心的港口群，確立“橫沙港—洋山深水港—南港”三足鼎立的格局，達(dá)到功能互補(bǔ)合理分配資源的目的。

另一方面，橫沙東灘可挖掘潛力大。第一，橫沙東灘已批準(zhǔn)促淤圈圍的面積為112 km2（1.7×105畝，大一區(qū)）。在其東側(cè)還約有270 km2（4.1×105畝）的灘涂存在，可作為中遠(yuǎn)期促淤圈圍規(guī)劃。第二，橫沙東灘北側(cè)有50 km以上岸線資源、緊貼北港航道，其中約14 km為目前10m以深的深水岸線，其余為目前7m左右深水岸線；南側(cè)約48 km岸線緊鄰長(zhǎng)江口北槽12.5m航道。隨著北港航道整治規(guī)劃實(shí)施，岸線水深將不斷增加，存在巨大的開發(fā)潛力，東部及東側(cè)亦可形成連接深水15～20m的大型港池或人工島。第三，航道資源豐富。橫沙東灘地處橫沙島，南貼長(zhǎng)江口北槽深水航道，北靠北港航道，西接長(zhǎng)江黃金水道，東臨東海。長(zhǎng)江口深水航道10m水深已向上貫通到南京長(zhǎng)江大橋，2010年3月長(zhǎng)江口的航道已達(dá)到12.5m水深，并將逐步向上延伸到南京長(zhǎng)江大橋。北港航道規(guī)劃為10m航道。橫沙通道水深維持10m。第四，有效可利用泥沙資源充足。長(zhǎng)江口豐水豐沙，大通站多年的年平均輸沙量為4.08×108t，近年來下泄泥量雖有所減少，但也有上億噸的泥沙下泄，長(zhǎng)江口口門地區(qū)的灘涂淤漲能力依然較強(qiáng)。第五，成陸速度快。在橫沙東灘南側(cè)，長(zhǎng)江口深水航道的年疏浚維護(hù)量可達(dá)7×107～8×107m3；在橫沙東灘北側(cè)，北港航道已在籌備建設(shè)，屆時(shí)也將產(chǎn)生大量的疏浚土。這些疏浚土就近吹填上灘，即可加快橫沙東灘的成陸過程，也可減少航道疏浚土的二次回淤現(xiàn)象，實(shí)現(xiàn)資源綜合利用。

因此，開發(fā)橫沙東灘，建設(shè)海港物流中心，完善國際航運(yùn)中心功能成為迫在眉睫的發(fā)展需要。從近期看，開發(fā)橫沙東灘，形成新港區(qū)及臨港產(chǎn)業(yè)區(qū)（長(zhǎng)江橋頭堡海港物流中心及臨海產(chǎn)業(yè)），與外高橋、洋山形成新的航運(yùn)中心格局；從長(zhǎng)遠(yuǎn)看，橫沙東灘、九段沙、南匯邊灘這些長(zhǎng)江口門區(qū)域的大片灘涂可綜合起來，形成陳吉余院士提出的長(zhǎng)江口亞三角洲，建設(shè)上海真正的港口群及海洋經(jīng)濟(jì)帶。

綜上，為了配合上海國際航運(yùn)中心的建設(shè)，利用長(zhǎng)江流域的疏浚土在橫沙島吹填土地，建立新的港口或?qū)⒊蔀橐粭l非常行之有效的道路。那么將面臨的數(shù)學(xué)問題就是，如何將長(zhǎng)江口岸的淤泥（圖1中兩條深色部分）吹填到它們所夾的粗線所圍區(qū)域內(nèi)，使得總工程耗費(fèi)最少。

圖1 橫沙東灘及長(zhǎng)江口淤泥示意圖Fig.1 East Hengsha shoaland the schematic diagram of siltat Yangtze River estuary

2 數(shù)學(xué)模型

將這個(gè)問題轉(zhuǎn)化為一個(gè)數(shù)學(xué)模型，就是困擾數(shù)學(xué)家232年的世界難題——“Monge-Kantorovich問題”。

1781年，Monge G提出了最優(yōu)輸運(yùn)問題（即Monge問題）[1]：假如要挖一個(gè)給定形狀的土坑，搬到另一處堆出一個(gè)給定形狀的建筑物，怎樣運(yùn)土可以使所作的功達(dá)到最小[2]，如圖2所示。

圖2 Monge問題背景Fig.2 M onge problem background

這個(gè)問題在一維的情況下早已經(jīng)解決。如果要把x處的土搬到y(tǒng)=f（x）處，記q（x）為x處的深度，p（y）為y處的設(shè)計(jì)高度，那么最優(yōu)的方案是f（x）是滿足p（y）dy=q（x）dx的單調(diào)遞增函數(shù)。但是在平面上的問題就要遠(yuǎn)遠(yuǎn)復(fù)雜了。主要原因是在平面上沒有像直線上的線性關(guān)系。

下面用數(shù)學(xué)語言把Monge問題重述如下：一個(gè)可分度量空間U，V上的測(cè)度μ、v和一個(gè)單位費(fèi)用函數(shù)C（x，y），同時(shí)記

問題轉(zhuǎn)化為：尋找一個(gè)變換T0∈Γ1使得

其中，T0為最優(yōu)映射，Ⅰ（T0）為最優(yōu)運(yùn)輸費(fèi)用。

通常單位費(fèi)用函數(shù)C（x，y）是非負(fù)的下半連續(xù)函數(shù)，常見的有：

線性費(fèi)用：C(x,y)=|x-y|；

平方費(fèi)用：C(x,y)=|x-y|2；

Lp費(fèi)用：C(x,y)=|x-y|p,p∈(1,+∞)；

C(x,y)=?(|x-y|)，?為某一凸函數(shù)。

其中|x-y|表示x、y之間的直線距離，Monge最初研究的是線性費(fèi)用情況。

自從Monge問題問世以來，它就成為了數(shù)學(xué)中幾個(gè)領(lǐng)域的數(shù)學(xué)家們所共同關(guān)心的經(jīng)典主題。微分幾何學(xué)家Appell P在1887年發(fā)表了一篇非常精彩的文章建立了一些最優(yōu)映射在平面和R3中的幾何性質(zhì)[3]。但由于Monge問題的不適定性，即使在一般歐式空間中，其最優(yōu)輸運(yùn)映射的存在性問題也不是那么顯而易見的。

為了解決這個(gè)困難，在1924年Kantorovich提出了輸運(yùn)問題的弱解形式[4]，Kantorovich當(dāng)時(shí)并不知道Monge的工作。Kantorovich問題簡(jiǎn)述為：找到U×V上的概率測(cè)度γ0∈Γ={在U×V上邊際分布為μ。v的概率測(cè)度全體}，使得：

他的本質(zhì)思想就是去尋找方案而不是輸運(yùn)弱解，也就是說，去尋找空間U×V中的概率測(cè)定，而不是兩個(gè)空間之間的映射。在兩者之間的聯(lián)系中，可以注意到，任意的輸運(yùn)映射T都可以推出一個(gè)方案γ，這個(gè)方案重點(diǎn)關(guān)注T在x×y中的圖譜，而且很容易證明反之亦成立。由于上述原因，任何一個(gè)輸運(yùn)都可以推出一個(gè)具有相同成本的方案這就導(dǎo)致了

因此，所有的Γ內(nèi)的極值點(diǎn)都是由輸運(yùn)推導(dǎo)而來的，那么可以直接從Kantorovich格式中得出輸運(yùn)映射的存在性。Kantorovich格式可以看成是原來Monge問題的一個(gè)弱解形式。

至此，人們將上述最優(yōu)輸運(yùn)問題統(tǒng)稱為“Monge-Kantorovich輸運(yùn)問題”，該問題又引發(fā)了學(xué)術(shù)界更多的關(guān)注，Kantorovich與他的學(xué)生Rubinstein 得到了著名的對(duì)偶公式[5，6]，將“Monge-Kantorovich問題”很好地引入了經(jīng)濟(jì)學(xué)、自動(dòng)化控制、運(yùn)輸學(xué)、流體力學(xué)、幾何學(xué)、形狀優(yōu)化、氣象學(xué)和金融數(shù)學(xué)等眾多領(lǐng)域。1991年Brenier Y用凸函數(shù)的梯度刻畫了最優(yōu)映射[7]，他的這篇文章建立了最優(yōu)輸運(yùn)問題與偏微分方程、概率論與泛函分析間的美妙聯(lián)系。從此最優(yōu)輸運(yùn)問題變得極其流行，應(yīng)用非常廣泛。

目前，這個(gè)問題的一維情況已經(jīng)被廣泛研究[2，8]，但是二維情況直到三年前才被鄭偉安和他的博士生沈銀芳所解決[9]，該文的想法是把原來的問題化為下列的二階擬線性橢圓偏微分方程的邊值問題。

這里未知函數(shù)u是一個(gè)二維的概率分布函數(shù)，ux與uy分別是u對(duì)x與對(duì)y的偏導(dǎo)數(shù)。藉此可以簡(jiǎn)單地算出最佳的平方距離E|X-Y|2，它的平方根給出Monge原問題的一個(gè)均方意義下的解。用此可以估算出原問題中欲求的“功”。

對(duì)這個(gè)數(shù)學(xué)問題感興趣的讀者可以參考陳木法院士的文章[10].

3 數(shù)學(xué)證明的簡(jiǎn)述

根據(jù)Kantorovich的表述，只需要找到兩個(gè)隨機(jī)的向量X=（X1，X2）與Y=（Y1，Y2），滿足下列條件：

1）X具有概率密度函數(shù)q（x1，x2），Y具有概率密度函數(shù)p（y1，y2）；

2）E|X-Y|2是所有滿足上述條件的隨機(jī)向量對(duì)中最小的。

文獻(xiàn)[9]中的想法是，假定（X，Y）就是欲找的最佳的一對(duì)?？梢詷?gòu)造新的隨機(jī)變量Z=（X1，Y2）。于是

有趣的是，此時(shí)X與Z只有一個(gè)分量不同，Z與Y也只有一個(gè)分量不同。如果知道Z的密度函數(shù)，則根據(jù)前述的一維解法，X與Z的關(guān)系就固定了，同理Z與Y的關(guān)系也固定了。于是問題簡(jiǎn)化為求Z的密度函數(shù)。更進(jìn)一步，如果Z的概率分布函數(shù)知道了，那它對(duì)體積的導(dǎo)數(shù)就是Z的密度函數(shù)，所以問題又進(jìn)一步簡(jiǎn)化為求Z的概率分布函數(shù)u如下。

實(shí)際上，E|X-Z|2與E|Z-Y|2都可以用一個(gè)含有u的積分式表示出來。所以用數(shù)學(xué)中常用的變分法就可以從（5）得到（4）。

圖3 區(qū)域網(wǎng)格劃分Fig.3 M esh p lotsof the region

4 數(shù)值模擬算法及結(jié)果

本文用文獻(xiàn)[9]中的結(jié)論作為求解的依據(jù)再加上線性規(guī)劃的思想來模擬計(jì)算工程量。在設(shè)計(jì)算法時(shí)，首先將要求解的區(qū)域分成若干個(gè)小正方形單元，如圖3所示，圖中共劃分了20×30個(gè)單元，每個(gè)小正方形單元的邊長(zhǎng)對(duì)應(yīng)實(shí)際情況中的2 km。劃分好網(wǎng)格后，將該實(shí)際問題轉(zhuǎn)化為可求解的數(shù)學(xué)模型，如圖4所示。

圖4 數(shù)學(xué)模型示意圖Fig.4 Schematic diagram of themathematicalmodel

設(shè)計(jì)算法如下：

1）用常用的線性規(guī)劃尋找每一大塊（4 km2的方格）搬運(yùn)的最優(yōu)方案；

2）在每個(gè)對(duì)應(yīng)的小單元中，用前面所述的“Monge-Kantorovich輸運(yùn)問題”解法分析每塊淤泥填到既定目標(biāo)區(qū)的均方意義下的功。

線性規(guī)劃問題：初始給定兩個(gè)矩陣(aij)N×N,(bij)N×N的最優(yōu)化問題可以化為

式中，表示aij的部分移動(dòng)到了bst,的花費(fèi)為。

其拉格朗日函數(shù)為：

由于問題的數(shù)據(jù)量較大，通過KKT條件直接求解矩陣的逆，計(jì)算量會(huì)很大。因此使用交替方向法[11]和并行運(yùn)算，求解該極小化問題：

初始化T0,λ0,μ0，選取步長(zhǎng)h

取k=1,2,3,…

計(jì)算{

取j=1,2,3,…

分塊計(jì)算{

更新：

而在每個(gè)對(duì)應(yīng)的小單元中，計(jì)算文獻(xiàn)[9]中的給出的變分問題：

其中

等價(jià)于求解極小化問題：

其中

r1，r2是常數(shù)。

使用交替下降法等計(jì)算方法，求解該變分問題的最優(yōu)值。算法如下：

初始化p0,q0,r0;

取k=1,2,3,…

計(jì)算{

更新：k=k+1．

在假設(shè)每平方千米單位深度的水里有a公斤泥沙，每頓泥沙搬運(yùn)一千米需要b元的情況下，用上述算法初步估計(jì)出工程總耗為22 000ab元。當(dāng)然這是在沒有任何水文資料與實(shí)際工程制約情況下的簡(jiǎn)單估算，實(shí)際問題難得多。

5 模型的擴(kuò)展

由于長(zhǎng)江口的淤泥是隨著時(shí)間不斷增長(zhǎng)的，而且還受到臺(tái)風(fēng)的影響，所以在實(shí)際應(yīng)用中會(huì)有更多的難度。因?yàn)檫@個(gè)問題還牽涉到了隨時(shí)間的變換。由于挖泥的形狀和次序影響以后新淤泥的堆積，目前還無法把它化為一個(gè)精確的數(shù)學(xué)題。但是可以從下列偏微分方程根據(jù)歷史數(shù)據(jù)確定懸浮泥沙的濃度。

式中，S為懸浮沙的濃度；u、v、w為x、y、z方向的流體速度；Kh、Kv為水平與豎直方向的擴(kuò)散系數(shù)；Ws為泥沙沉降速度。

把這個(gè)方程與前面提到的邊界問題結(jié)合起來，就可以把問題分時(shí)間階段化為前面的“Monge-Kantorovich輸運(yùn)問題”，例如把每三個(gè)月的淤泥形狀用前面的方法處理一下，再把它們累計(jì)起來，用計(jì)算機(jī)得到近似優(yōu)化值。

6 結(jié)語

隨著上海國際航運(yùn)中心的龍頭地位日益面臨嚴(yán)峻挑戰(zhàn)，為了配合上海國際航運(yùn)中心的建設(shè)，利用長(zhǎng)江口流域的疏浚土在橫沙島吹填土地，建立新的港區(qū)或?qū)⒊蔀橐粭l非常行之有效的道路。本文以長(zhǎng)江口疏浚土吹填橫沙東灘為背景，將其轉(zhuǎn)化為數(shù)學(xué)模型：“Monge-Kantorovich輸運(yùn)問題”，以此設(shè)計(jì)算法，進(jìn)而得到最有運(yùn)輸條件下的總工程耗費(fèi)。由于工程本身的規(guī)模巨大，通過尋找最優(yōu)的輸運(yùn)方案，能夠最大限度地減少施工費(fèi)用支出，降低工程造價(jià)成本。為了更好地模擬實(shí)際操作過程中所可能遇到的問題，全文最后還討論了模型進(jìn)一步的擴(kuò)展提升，提出了可將其轉(zhuǎn)化為一個(gè)與時(shí)間變化有關(guān)的問題，以加強(qiáng)模型的適用性。

[1] Monge G.Memoire sur la theorie des déblais et des remblais[J].Pairs ：Histoire de I′Acad.des Sciencesde Pairs，1781.666-704．

[2] Villani C.Optimal Transport[M].New York：Springer-Verlag，2009.

[3] Appell P.Mémoire sur les déblais et les remblais des systémes continus ou discontinues[J].Mémoires présentés par divers Savantsà l′Académie des Sciences de l′Institut de France，1887，29：1-208．

[4] Kantorovich L V.On a problem of Monge[J].UspehiMat Nauk，1948，3：225-226．

[5] Kantorovich L V.On the transfer of mass[J].Dokl Akad Nauk，SSSR.1942，37：227-229．

[6] Rubinstein G S，Kantorovich LV.On the space of completely additive functions[J].Vestnic Leningrad Univ，Ser MatMekh i Astron.1958，13（7）：52-59．

[7] Brenier Y.Polar factorization and monotone rearrangement of vector-valued functions[J].Comm Pure Appl.Math，1991，44：375-417．

[8] Rachev ST，Ruschendorf L.Mass Transportation Problems Volume I：Theory ，Volume II：Applications[M].New York ：Springer-Verlag，1998.

[9]Shen Y F，ZhengW A.On Momge-Kantorovich problem in the plane[J].CR Acad SciPairs，2010，1：267-271．

[10] 陳木法.談?wù)劯怕收撆c其他學(xué)科的若干交叉[J].數(shù)學(xué)進(jìn)展，2005，34（6）：661-672．

[11] Boyd S，Parikh N，Chu E，et al.Distributed optim ization and statistical learning via the alternating directionmethod ofmultipliers[J].Machine Learning，2010，3：1-122．