成毅 何晨 張?jiān)?/p>

1中國(guó)石油國(guó)際勘探開(kāi)發(fā)有限公司

2國(guó)家管網(wǎng)油氣調(diào)控中心

我國(guó)的地理環(huán)境是西高東低，山地、高原和丘陵約占陸地面積的67%，盆地和平原約占陸地面積的33%，這就決定了在我國(guó)的油氣管道項(xiàng)目中山區(qū)段的管道占有較大比重。受多種因素的制約，管道在山區(qū)管段的建設(shè)投資和運(yùn)營(yíng)代價(jià)要遠(yuǎn)高于平原地區(qū)的管道。因此，從2013 年開(kāi)始，為了滿(mǎn)足生產(chǎn)實(shí)際和降低成本以應(yīng)對(duì)低迷油價(jià)的需要，部分學(xué)者開(kāi)始研究三維地形條件下的管道路徑優(yōu)選問(wèn)題。

周軍針對(duì)氣田管網(wǎng)，基于數(shù)字化地形數(shù)據(jù)，在考慮了各類(lèi)環(huán)境因素以后，采用A*智能算法確定了氣體管道最優(yōu)路徑，并將曲面路徑的長(zhǎng)度與路徑投影為直線的長(zhǎng)度進(jìn)行了比較[1]。張思琦在三維地形條件下尋找最優(yōu)路徑時(shí)，將障礙、坡度等因素考慮在內(nèi)，障礙區(qū)域的高程值轉(zhuǎn)化成相對(duì)于整個(gè)區(qū)域較大的固定數(shù)值，從而使遺傳算法的適應(yīng)度函數(shù)值降低而在迭代的工程中被淘汰，進(jìn)而達(dá)到了避開(kāi)障礙的目的[2]。雷揚(yáng)首先將障礙區(qū)域提取成凸多邊形，然后采用模糊聚類(lèi)分析把注水區(qū)域劃分成多個(gè)子區(qū)域，最后利用智能算法確定了最優(yōu)路徑以及避障的距離[3]。

針對(duì)三維地形條件下的管道路徑優(yōu)選問(wèn)題，國(guó)內(nèi)外學(xué)者已經(jīng)開(kāi)展了大量的研究，并取得了豐碩的成果，但也存在一些難點(diǎn)和不足。三維地形的準(zhǔn)確表征是解決三維地形條件下管道路徑優(yōu)選問(wèn)題的前提[4]。現(xiàn)有研究主要采用三維曲面拼接的方式來(lái)表征地形，但是這種方式往往非常依賴(lài)曲面函數(shù)，魯棒性不理想，而且模擬出的地形與實(shí)際地形誤差較大（圖1，圖2）。因此，基于OpenGL 技術(shù)，使用空間插值方法建立DEM 模型從而進(jìn)行三維地形表征，并采用D-P 算法進(jìn)行路徑優(yōu)選。模擬實(shí)驗(yàn)結(jié)果表明，采用OpenGL 建立起的DEM 模型進(jìn)行路徑優(yōu)選與用三維曲面建立的地形模型相比，前者誤差更小，更加接近真實(shí)地形。D-P 算法可以有效地減少三維地形條件下管道長(zhǎng)度，降低管道建設(shè)投資費(fèi)用。

圖1 三維地形條件下管道路徑Fig.1 Pipeline path under 3D terrain condition

圖2 真實(shí)三維地形圖Fig.2 Real 3D topographic map

1 DEM 模型

DEM 是是當(dāng)今地理學(xué)、氣象學(xué)、計(jì)算機(jī)計(jì)算科學(xué)，尤其是地理信息科學(xué)研究的熱點(diǎn)[5]。DEM 是表示某一區(qū)域地形D 上的三維向量有限序列，用函數(shù)的形式可以描述為

DEM 模型有兩種表達(dá)形式：GRID 模型和TIN模型。TIN 模型在存儲(chǔ)數(shù)據(jù)時(shí)，不僅要存儲(chǔ)各個(gè)三角形頂點(diǎn)的高程信息，而且要存儲(chǔ)三角形各個(gè)頂點(diǎn)的拓?fù)潢P(guān)系，因此數(shù)據(jù)結(jié)構(gòu)復(fù)雜，效率較低[6]。另外，在構(gòu)建三角網(wǎng)時(shí)需要用到較為復(fù)雜的算法，這更增加了建模的難度。GRID 模型在存儲(chǔ)數(shù)據(jù)時(shí)，只需要存儲(chǔ)網(wǎng)格點(diǎn)的起始坐標(biāo)，網(wǎng)格行、列數(shù)以及各個(gè)網(wǎng)格節(jié)點(diǎn)的高程信息，因此，數(shù)據(jù)存儲(chǔ)和讀取方便，易于使用[7]。所以采用GRID 形式進(jìn)行DEM模型建模（圖3）。

圖3 GRID 模型Fig.3 GRID model

1.1 GRID 模型建模步驟

在現(xiàn)實(shí)生活中，采樣得到的數(shù)據(jù)往往是離散的隨機(jī)分布的，因此需要將離散數(shù)據(jù)網(wǎng)格化。離散數(shù)據(jù)網(wǎng)格化的一般做法是：根據(jù)精度要求，對(duì)地形進(jìn)行網(wǎng)格劃分，然后再利用空間插值的方法計(jì)算出每一個(gè)網(wǎng)格節(jié)點(diǎn)的值，進(jìn)而形成三維地形。GRID 模型建模流程如圖4 所示。

圖4 GRID 模型建模流程Fig.4 GRID model modeling process

1.2 GRID 模型數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)

GRID 模型中需要存儲(chǔ)的高程數(shù)據(jù)可以表示如公式（2）、（3）、（4）所示：

式中：X0、Y0為DEM 起始點(diǎn)橫縱坐標(biāo)；Dx、Dy為x方向和y方向上的間隔；Nx-1，Ny-1為DEM 網(wǎng)格的行列數(shù)。

可以看出，GRID 在存儲(chǔ)數(shù)據(jù)時(shí)，只需要存儲(chǔ)初始點(diǎn)坐標(biāo)，網(wǎng)格的間距，網(wǎng)格點(diǎn)的高程以及行、列數(shù)就可以得到當(dāng)前任何地形區(qū)域的高程信息。

2 模型建立

2.1 目標(biāo)函數(shù)

建立三維地形條件模型的目的是使管道投資建設(shè)費(fèi)用最少，這里的費(fèi)用主要與管道的長(zhǎng)度和管徑有關(guān)。為了簡(jiǎn)化計(jì)算和更好的關(guān)注問(wèn)題本身，假設(shè)管道的管徑一定，則管道的投資建設(shè)費(fèi)用只與管道的長(zhǎng)度有關(guān)，因此，以管道總長(zhǎng)度為目標(biāo)函數(shù)，建立如下數(shù)學(xué)模型:

式中：L為管道長(zhǎng)度，m；Xi，Yi，Zi為搜尋第i個(gè)網(wǎng)格時(shí)的高程數(shù)據(jù)，m。

2.2 約束條件

彈性敷設(shè)就是不使用彎頭，而是利用管道自身的彈性逐漸改變管道的走向。由于其具有應(yīng)力分布均勻，不存在過(guò)高的峰值應(yīng)力以及不需設(shè)置固定墩，施工方便等優(yōu)點(diǎn)，所以成為天然氣管道在起伏地區(qū)鋪設(shè)的重要方式。彈性敷設(shè)時(shí)，管道的曲率半徑應(yīng)該滿(mǎn)足管道的強(qiáng)度要求，垂直面上彈性敷設(shè)管道的曲率半徑還應(yīng)該大于管道在自重作用下產(chǎn)生的撓度曲線的曲率半徑，曲率半徑可以按照式（6）計(jì)算。

式中：R為管道彈性彎曲曲率半徑，m；α為管道的轉(zhuǎn)角，(°)；D為鋼管外徑，cm。

3 求解算法

3.1 Dijkstra 算法

Dijkstra 算法是目前公認(rèn)的效果較好的最短路徑尋優(yōu)算法[8]，它是典型的單源最短路徑算法，其在現(xiàn)實(shí)世界中的應(yīng)用非常廣泛，例如在道路管網(wǎng)鋪設(shè)中的運(yùn)用，智能交通系統(tǒng)中的運(yùn)用，管道鋪設(shè)中的運(yùn)用。它的基本思想是：對(duì)路徑中可能經(jīng)過(guò)的各個(gè)節(jié)點(diǎn)進(jìn)行編號(hào)，然后遍歷各個(gè)節(jié)點(diǎn)，計(jì)算出起點(diǎn)到終點(diǎn)可能經(jīng)過(guò)的路徑的距離，最后找出最短路徑[9]。由于遍歷了所有的節(jié)點(diǎn)，因此一定可以找到一條最短路徑。圖5 為一個(gè)無(wú)向帶權(quán)聯(lián)通圖，利用Dijkstra 算法，可求出v0與v3的最短路徑。

圖5 無(wú)向帶權(quán)聯(lián)通圖Fig.5 Undirected weighted connection graph

要從v0到達(dá)v3有5條路徑：v0v3，v0v1v3，v0v1v2v3，v0v4v3，v0v5v4v3，這五條路徑耗費(fèi)的距離分別為5、10、15、3、6，由此可以看出從v0到達(dá)v3的最短路徑為v0v4v3，距離為3。

從上述利用Dijkstra 算法求最短路徑的過(guò)程可以看出，Dijkstra 算法不限制搜索方向，可見(jiàn)其搜索量巨大，時(shí)間復(fù)雜度較高；而且Dijkstra 算法采用的是鄰接矩陣數(shù)據(jù)構(gòu)，占用的內(nèi)存較大，空間復(fù)雜度較高，這些都無(wú)疑會(huì)降低Dijkstra 算法的路徑尋優(yōu)效率。

3.2 人工勢(shì)場(chǎng)法

1986 年，Khatib 提出了人工勢(shì)場(chǎng)法，其基本思想是將機(jī)器人的規(guī)劃空間看成一個(gè)受障礙物的斥力場(chǎng)和目標(biāo)位置的引力場(chǎng)共同作用下的人工勢(shì)場(chǎng)，在斥力場(chǎng)和引力場(chǎng)的共同作用下，移動(dòng)機(jī)器人從高勢(shì)點(diǎn)向低勢(shì)點(diǎn)移動(dòng)，最后完成路徑尋優(yōu)[10]。

人工勢(shì)場(chǎng)分為兩部分：①由目標(biāo)位置產(chǎn)生的引力場(chǎng)，場(chǎng)強(qiáng)與機(jī)器人到目標(biāo)點(diǎn)的距離成正比，即機(jī)器人到目標(biāo)點(diǎn)的距離越大，場(chǎng)強(qiáng)越大；機(jī)器人到目標(biāo)點(diǎn)的距離越小，場(chǎng)強(qiáng)越小。機(jī)器人的受力方向是由機(jī)器人指向目標(biāo)點(diǎn)（Goal）[11]；②由障礙物產(chǎn)生的斥力場(chǎng)，場(chǎng)強(qiáng)與機(jī)器人到障礙物的距離成反比，即機(jī)器人到障礙物的距離越大，場(chǎng)強(qiáng)越小，機(jī)器人到障礙物的距離越小，場(chǎng)強(qiáng)越大，機(jī)器人的受力方向是由障礙物（Obstacle）指向機(jī)器人[12]。

機(jī)器人在人工勢(shì)場(chǎng)中的受力如圖6 所示。Obstacle 為障礙物，Start 為起點(diǎn)位置，Goal 為目標(biāo)點(diǎn)位置，F(xiàn)rep 為障礙Obstacle 對(duì)機(jī)器人的斥力，F(xiàn)att為目標(biāo)點(diǎn)對(duì)機(jī)器人的吸引力，F(xiàn)total 為障礙和目標(biāo)點(diǎn)對(duì)機(jī)器人的合力。由于在機(jī)器人尋找最短路徑的過(guò)程中，其受到斥力和引力兩方面因素影響，因此會(huì)被限制搜索方向，提高了搜索效率。

圖6 人工勢(shì)場(chǎng)法受力圖Fig.6 Artificial potential field stress diagram

綜上可得，Dijkstra 由于遍歷了整個(gè)區(qū)域的節(jié)點(diǎn)，因此，其必然可以找到一條最短路徑，但也正因此，在數(shù)據(jù)量巨大的情況下，其搜索效率較低[13]。人工勢(shì)場(chǎng)法在搜索路徑時(shí)，由于其受到斥力場(chǎng)和引力場(chǎng)兩方面因素的影響，因此可以限制搜索方向，提高了搜索的效率[14]?；贒ijkstra 算法和人工勢(shì)場(chǎng)法，提出了D-P 路徑搜索算法。

3.3 D-P 搜索算法步驟

（1）三維地形的建立。根據(jù)精度要求，對(duì)地形進(jìn)行網(wǎng)格劃分，然后再利用反距離權(quán)重插值的方法計(jì)算出每一個(gè)網(wǎng)格節(jié)點(diǎn)的值，進(jìn)而形成三維地形。樹(shù)形結(jié)構(gòu)在存儲(chǔ)數(shù)據(jù)時(shí)具有獨(dú)特的優(yōu)勢(shì)，它可以減少數(shù)據(jù)量，提高查找效率。因此，采用四叉樹(shù)數(shù)據(jù)結(jié)構(gòu)來(lái)存儲(chǔ)網(wǎng)格的起始點(diǎn)坐標(biāo)、網(wǎng)格間距以及各個(gè)網(wǎng)格交點(diǎn)的高程值信息。

（2）基準(zhǔn)路徑的建立。連接起點(diǎn)和目標(biāo)點(diǎn)，形成一條基準(zhǔn)路徑。以合適的寬度將基準(zhǔn)路徑拓寬，形成一條基準(zhǔn)路徑帶。

（3）地形復(fù)雜度的獲取。地形復(fù)雜度是描述地形陡峭崎嶇的重要指標(biāo)，不同的地形因子刻畫(huà)了不同的地形幾何特征。常見(jiàn)的地形復(fù)雜度因子有局部高差、局部標(biāo)準(zhǔn)差、局部褶皺度、局部全曲率四種。局部高差即目標(biāo)區(qū)域內(nèi)最高點(diǎn)與最低點(diǎn)的高程值之差，它反映了目標(biāo)區(qū)域內(nèi)地形的起伏程度。局部標(biāo)準(zhǔn)差就是目標(biāo)區(qū)域內(nèi)高程的標(biāo)準(zhǔn)差，它反映的是目標(biāo)區(qū)域內(nèi)地形平均起伏程度。局部褶皺度是地形的三維面積與投影平面面積之比，它描述了地形的平均褶皺程度。局部全曲率反映的是目標(biāo)地形的平均突變程度。利用人工勢(shì)場(chǎng)法搜索最短路徑，必須知道障礙區(qū)域的位置，因此采用局部全曲率來(lái)計(jì)算地形復(fù)雜度。

局部全曲率的表達(dá)式如式（7）所示。

式中：r為x方向的二階偏導(dǎo)數(shù)；t為y方向的二階偏導(dǎo)數(shù)；s為x，y方向上的偏導(dǎo)數(shù)。

利用局部全曲率計(jì)算出地形復(fù)雜度后，與預(yù)先設(shè)定好的閾值進(jìn)行比較，得出障礙的位置，障礙記為Obstacle。

（3）利用D-P 算法進(jìn)行路徑優(yōu)選。首先利用人工勢(shì)場(chǎng)法確定路徑搜索的方向，設(shè)搜索點(diǎn)位置為X=(x，y，z)，引力勢(shì)函數(shù)定義為

式中：Uatt為目標(biāo)點(diǎn)的引力場(chǎng)；k大于0 的引力勢(shì)場(chǎng)常量；X為搜索點(diǎn)的位置向量；Xg為搜索點(diǎn)在勢(shì)力場(chǎng)中的目標(biāo)點(diǎn)位置。

斥力場(chǎng)函數(shù)定義為

式中：Urep為障礙物的斥力場(chǎng)；μ為斥力場(chǎng)的常量；ρ為搜索點(diǎn)在空間位置與障礙物的最短距離；ρ0為單個(gè)障礙物影響的最大距離范圍，當(dāng)搜索點(diǎn)與障礙物的位置大于ρ0時(shí)，障礙物對(duì)于搜索點(diǎn)的運(yùn)動(dòng)不再產(chǎn)生影響。

在確定了搜索方向之后，利用Dijkstra 算法進(jìn)行路徑優(yōu)選。

4 實(shí)例

4.1 地形模擬

在Google Earth 中以1 200 m×1 800 m 為范圍隨機(jī)提取2×104個(gè)點(diǎn)的高程數(shù)據(jù)，然后以10 m 為分辨率對(duì)該區(qū)域進(jìn)行網(wǎng)格劃分，最后利用反距離權(quán)重插值法計(jì)算出每一個(gè)網(wǎng)格節(jié)點(diǎn)的高程值，生成三維地形（圖7）。管道起、終點(diǎn)坐標(biāo)以及邊界坐標(biāo)見(jiàn)表1。

圖7 OpenGL 建立的三維地形網(wǎng)格圖Fig.7 3D terrain grid map established by OpenGL

表1 管道起、終點(diǎn)坐標(biāo)以及邊界坐標(biāo)Tab.1 Starting,ending,and boundary coordinates of the pipeline

從建立的DEM 模型中隨機(jī)提取15 個(gè)點(diǎn)的高程數(shù)據(jù)，與從三維曲面生成的地形中提取的高程數(shù)據(jù)以及真實(shí)高程數(shù)據(jù)進(jìn)行對(duì)比結(jié)果如表2。

表2 高程數(shù)據(jù)及誤差對(duì)比Tab.2 Elevation data and error comparison

高程數(shù)據(jù)對(duì)比結(jié)果如圖8 所示，誤差對(duì)比結(jié)果如圖9 所示。從圖8 和圖9 可以看出，基于OpenGL建立的DEM 模型的高程數(shù)據(jù)與真實(shí)值的誤差比采用三維曲面拼接建立起的地形模型的高程數(shù)據(jù)誤差小，精度也更高。

圖8 高程數(shù)據(jù)對(duì)比Fig.8 Elevation data comparison

圖9 誤差對(duì)比Fig.9 Error comparison

4.2 路徑優(yōu)化

采用D-P 算法從管道的起點(diǎn)開(kāi)始搜索，得到起伏地形條件下的天然氣管道優(yōu)化路徑（圖10）。

圖10 地形起伏地區(qū)天然氣管道優(yōu)化路徑與直線路徑Fig.10 Optimal route and straight-line route of natural gas pipeline in undulating terrain

圖10 中紅色線代表管道直線路徑，黑色線代表優(yōu)化路徑。計(jì)算得到起伏地形區(qū)域天然氣管道的直線路徑和優(yōu)化路徑的總長(zhǎng)度如表3 所示。

表3 起伏地形區(qū)域直線路徑長(zhǎng)度和優(yōu)化路徑長(zhǎng)度對(duì)比Tab.3 Comparison of straight path length and optimized path length in undulating terrain area

分析上述實(shí)驗(yàn)數(shù)據(jù)可以得到:

（1）優(yōu)化后的路徑在三維空間上的總長(zhǎng)度比直線路徑在三維空間上的總長(zhǎng)度少208.4 m，優(yōu)化率達(dá)到了10.6%。

（2）由圖10 可以看出，優(yōu)化后的路徑避開(kāi)了起伏較大的障礙區(qū)域，沿程的地形起伏更加平緩。

5 結(jié)論

（1）基于OpenGL 建立的DEM 模型和通過(guò)三維曲面拼接建立起的地形模型相比，前者模型的高程數(shù)據(jù)更加接近地形真實(shí)地形，誤差更小，精度更高。

（2）基于Dijkstra算法和人工勢(shì)場(chǎng)法提出的D-P算法不僅可以降低管道總長(zhǎng)度，減少管道投資費(fèi)用，而且還可以有效地降低管道沿線的地形起伏程度。