中國東部沿海管道海岸—斷裂帶復(fù)合效應(yīng)

方莉，夏鵬程，梁志珊*，熊樹海，畢武喜，藍(lán)衛(wèi)

1 中國石油大學(xué)(北京)信息科學(xué)與工程學(xué)院，北京 102249

2 中石油管道科技研究中心，廊坊 065000

0 引言

太陽活動和星際空間共同作用生成地球的磁層—電離層電流，這個變化的電流產(chǎn)生變化的地磁場，在埋地管道設(shè)備中引起地磁感應(yīng)電流(Geomagnetically Induced Current，GIC)，直接影響管道與土壤之間的管地電位(pipe-soil potential，PSP)。管道GIC和PSP異常加速管道設(shè)備的腐蝕，威脅管道輸送系統(tǒng)的安全運(yùn)行，增大了維護(hù)成本[1]。

近年來，中國石油管道公司在管道運(yùn)行中發(fā)現(xiàn)，山東日照陰極保護(hù)站管道PSP干擾大，恒電位儀經(jīng)常出現(xiàn)異常甚至報警[2]，若管道長時間處于欠保護(hù)狀態(tài)，漏電處將加快腐蝕。因此，油氣管道地磁暴干擾監(jiān)測和地磁暴對油氣管網(wǎng)和電網(wǎng)的致災(zāi)機(jī)理研究被中石油管道公司列為“十二五”重點(diǎn)研發(fā)計劃項(xiàng)目。中國石油大學(xué)(北京)承擔(dān)了埋地管道研究，研制了國內(nèi)首臺埋地管道GIC和PSP監(jiān)測裝置并在中石油西氣東輸、陜京線、日東線等沿線陰保站投入使用[3]，為石油管道GIC和PSP預(yù)測提供了實(shí)測數(shù)據(jù)支持。同時，使用回路電流法計算管道GIC和PSP[4]。經(jīng)調(diào)查研究后發(fā)現(xiàn)，日照東臨黃海，西側(cè)有郯廬大斷裂[5]，猜測日照站PSP異常原因是受到海岸和斷裂帶兩者的相互作用。對于沿海地區(qū)存在斷裂帶的復(fù)雜地質(zhì)情況下，位于海岸—斷裂帶之間區(qū)域電場受到海岸和斷裂帶雙重電導(dǎo)率橫向突變的影響，本文將之稱為海岸—斷裂帶復(fù)合效應(yīng)。

國外，Boteler D H和Pirjola R等學(xué)者早前研究表明地磁感應(yīng)電流在電網(wǎng)和埋地管道中的存在和危害[6]。2002 年Boteler D H采用平面波法通過分層波阻抗法得到大地電場對均勻長輸管道模型(DSTL)的管道PSP在時域的變化[7]，此方法和復(fù)鏡像法[8]都忽略了地質(zhì)結(jié)構(gòu)電導(dǎo)率的橫向變化，不適用于沿海地區(qū)。Gilbert和Weaver[9-10]后來提出了一維、二維和薄殼模型計算半空間線電流在海岸線附近的地電場分布，尚未計算三維空間電場分布。Kuvshinov[11-12]應(yīng)用快速積分方程法計算海岸附近的地電場，但是解析法求解地面電磁場只適用于簡單形態(tài)地質(zhì)模型。2013 年Kuvshinov基于衛(wèi)星磁測數(shù)據(jù)建立全球模式下的三維磁暴感應(yīng)地電場模型[13]，研究只針對海岸附近電場異常。國內(nèi)已有研究[14-15]針對埋地管道GIC的特性和評估，劉連光等[16]基于大地分層模型和均勻傳輸線管道模型建立了管道PSP磁暴干擾評估大地模型，但這兩種模型均存在缺陷，大地分層模型未考慮到橫向電導(dǎo)率差異，均勻傳輸線管道模型中未涉及大地電阻。綜上，目前國內(nèi)外還未有對沿海復(fù)雜地質(zhì)下埋地管道PSP時空分布特性的研究報道。這就需要探索一種適應(yīng)性更強(qiáng)、具備可操作性和實(shí)際預(yù)測意義的計算方法。

埋地管道PSP研究需要分別建立地電場計算模型和管道模型，有限元法(FEM)是目前較為簡單實(shí)用且有效的電磁數(shù)值分析方法，適應(yīng)各種復(fù)雜地質(zhì)結(jié)構(gòu)且精度較高，已成功應(yīng)用橫向電導(dǎo)率不均勻大地電場研究[17]。傅里葉變換(FT)分析合成技術(shù)通過對數(shù)據(jù)的時頻轉(zhuǎn)換，對信號實(shí)時處理發(fā)揮重要作用，為大地電磁實(shí)時仿真提供了理論方法。LZS-DSPL模型[18]在Boteler的DSTL模型基礎(chǔ)上進(jìn)行改進(jìn)，更加符合實(shí)際情況。因此本文提出PFFEM法(Pipeline Finite element method based on Fourier transform)，針對2016 年10 月12-14 日山東地區(qū)日東線管道PSP地磁擾動響應(yīng)，以地磁臺提供的實(shí)際記錄磁場數(shù)據(jù)為基礎(chǔ)，計算日東線管道PSP分布，將沿海日照站與內(nèi)陸兗州站的計算值與監(jiān)測值進(jìn)行對比，發(fā)現(xiàn)位于海岸-斷裂帶區(qū)域的埋地管道PSP的幅值水平和波動幅度更大，結(jié)果揭示了海岸-斷裂帶復(fù)合效應(yīng)下沿海埋地管道PSP的分布特性。

1 三維埋地管道PFFEM 理論

如圖1 所示，建立三維大地埋地管道PSP預(yù)報的PFFEM理論模型。該計算模型基于傅里葉變換的管道PSP有限元算法，以電磁Maxwell T-Ω方程三維有限元法為基礎(chǔ)[19]，結(jié)合快速傅里葉變換和反變換的數(shù)學(xué)思想，以地面磁場數(shù)據(jù)正演埋地管道PSP時空分布，將此方法稱為管道傅里葉—有限元法(PFFEM，Pipeline Finite element method based on Fourier transform)方法。PFFEM方法將時域磁場H(t)經(jīng)快速傅里葉變換FFT分解成頻域磁場H(w)(t為時間，w為頻率)，作為磁場源輸入三維電導(dǎo)率模型進(jìn)行有限元計算得到頻域電場E(w)，電場經(jīng)傅里葉反變換IFFT后合成時域電場E(t)，經(jīng)過管道PSP計算方法得到PSP時空分布結(jié)果。

圖1 三維埋地管道PFFEM計算方法流程圖Fig.1 Method steps of PFFEM

計算管道PSP時空分布方法主要分為以下4 個步驟：

(1)大地磁場FFT分解

將地磁水平分量和垂直分量分別作為大地模型的初始磁源，那么首先對地磁臺磁場數(shù)據(jù)東向和北向分量hx、hy進(jìn)行分解為nx、ny個諧波之和，即

式中，Hx,i、Hy,i分別表示時域hx、hy第i個諧波的幅值，ωx,i、ωy,i表示hx、hy第i個諧波的頻率，φx,i、φy,i表示hx、hy的第i個諧波的相位。

(2)大地電場FEM計算

磁暴地磁場的頻率為0.0001～0.01 Hz，因此可看作磁準(zhǔn)靜態(tài)場，與渦流電流密度相比，位移電流可以略去不計[20]。一次電流源Jp產(chǎn)生了一次場(一次電場Ep和一次磁場Hp)，一次場在大地導(dǎo)體中產(chǎn)生二次場(二次電場Es和二次磁場Hs)，內(nèi)部無源。時諧磁場中以e-jωt為時諧因子，其中j為虛數(shù)單位。則二次電磁場滿足

式中，假設(shè)介質(zhì)的磁導(dǎo)率為真空磁導(dǎo)率μ0=4π×10-7A/m，介質(zhì)的介電常數(shù)為真空介電常數(shù)ε0=8.85×10-12F/m，γ表示導(dǎo)體電導(dǎo)率。

由赫姆霍茲定理可知導(dǎo)電區(qū)電場和磁場有唯一解[21]，在渦流場求解控制方程，根據(jù)區(qū)域結(jié)構(gòu)和邊界條件，需要考慮計算機(jī)運(yùn)行時間和計算精度選擇適當(dāng)?shù)碾姶盼粚ΓＲ姺椒ㄓ蠥法和T法。T法引入位函數(shù)T和Ω，T為矢量電位，Ω為標(biāo)量磁位，求解過程未知數(shù)相比于A法較少，且能達(dá)到滿意精度。根據(jù)麥克斯韋方程，J、H與T、Ω之間的關(guān)系是

其中，J表示渦流場電流源，H表示渦流場磁場，E表示渦流場電場。

式(4)中“－?Ω”項(xiàng)為常量，記為H0，表示電流源在導(dǎo)體表面產(chǎn)生T0的等效磁場源。

對于地磁H0,x分解到的每一個余弦分量記為Hi(ωi)，將分解的每個磁場分量作為有限元計算需要的磁場計算模型的地表面邊界條件，對大地求解區(qū)進(jìn)行四面體網(wǎng)格自適應(yīng)劃分，在滿足最大迭代步數(shù)和誤差精度的迭代要求下，求解方程得到頻域大地電場x方向和y方向分量。

(3)管道電場IFFT合成

式(6)中，Ex,i,j、Ey,i,j為第i個頻率下管道j點(diǎn)處的電場x、y分量幅值，ωx′,i、ω′y,i為第i個頻率下管道j點(diǎn)處的電場x、y分量的相位，φx′,i、φy′,i為第i個頻率下管道j點(diǎn)處的電場x、y分量的相角。

在管道j點(diǎn)處的x、y電場分量值由磁場x、y分量產(chǎn)生的電場分量對應(yīng)相加得到

則t時刻，管道j點(diǎn)處的電場幅值為

(4)管道PSP計算

取出管道每千米間隔的地電場值代入中石油日東線分布源傳輸線LZS-DSPL模型，將地磁作用等效至管道電路模型進(jìn)行計算。把單位長度管道作為一個基本單元，一條長直管道可看作n段基本單元串聯(lián)；采用基爾霍夫定律回路電路法或節(jié)點(diǎn)電壓法建立GIC和PSP矩陣方程，得到管道管地電位PSP在磁暴期間的時空分布。

2 海岸—斷裂帶復(fù)合效應(yīng)

為分析斷裂帶電導(dǎo)率對海岸—斷裂帶區(qū)域的電場幅值變化的影響，我們假定海岸—斷裂帶三維模型為長方體結(jié)構(gòu)，長1000 km，寬200 km，深500 km。海岸線位于模型上表面中央，斷裂帶在距離海岸250 km處，斷裂帶區(qū)域深度為25 km，寬度為250 km，海水區(qū)域深度為5 km，寬度為500 km，將海水電導(dǎo)率取為4 S/m，大地連同斷裂帶底層和海洋底層一體，取大地電導(dǎo)率為0.0025 S/m，如圖2 所示。

圖2 海岸—斷裂帶電導(dǎo)率模型示意圖Fig.2 Schematic diagram of the conductivity model of coastal-fault zone

對模型上表面施加平行于海岸線方向的磁場H，幅值為0.8 A/m。考慮到斷裂帶電導(dǎo)率對感應(yīng)電場的影響，將斷裂帶區(qū)域電導(dǎo)率參數(shù)分別設(shè)為0.1 S/m、0.01 S/m、0.001 S/m和0.0001 S/m以 及0.0025 S/m，分別進(jìn)行有限元計算后取出以上5 種參數(shù)模型中橫跨斷裂帶與海岸線的上表面中線處(圖2 中虛線)電場值。根據(jù)斷裂帶與陸地電導(dǎo)率之間的大小關(guān)系，繪制成2 種對比圖，結(jié)果如圖3 所示。

圖3 不同斷裂帶電導(dǎo)率下的海岸—斷裂帶效應(yīng)電場Fig.3 Electric field effect of coastal-fault zone under different conductivity

假定無斷裂帶時，即斷裂帶與大地電導(dǎo)率都為0.0025 S/m時，由圖3(a)和圖3(b)可以看出，電場在陸地側(cè)由海岸線向遠(yuǎn)處逐漸下降穩(wěn)定至0.32 V/km。

當(dāng)存在斷裂帶時，有2 種情況分別討論如下：

(1)當(dāng)斷裂帶區(qū)域電導(dǎo)率高于中間陸地電導(dǎo)率時，海岸—斷裂帶之間的陸地均在海岸和斷裂帶邊界處產(chǎn)生極大值，斷裂帶對中間區(qū)域電場有“提升”作用，且斷裂帶電導(dǎo)率越高，海岸—斷裂帶區(qū)域電場“提升”的程度越大，如圖3(a)所示；

(2)當(dāng)斷裂帶區(qū)域電導(dǎo)率低于中間陸地電導(dǎo)率時，海岸—斷裂帶區(qū)域電場在海岸處產(chǎn)生極大值，在斷裂帶邊界處產(chǎn)極小值，斷裂帶對中間區(qū)域電場有“下拉”作用，且斷裂帶電導(dǎo)率越低，海岸—斷裂帶區(qū)域電場“下拉”越明顯，如圖3(b)所示。

表1 分別列出海陸邊界處陸地側(cè)、陸地—斷裂帶邊界處陸地側(cè)和陸地—斷裂帶邊界處斷裂帶側(cè)的電場極值，從表中數(shù)據(jù)得出，在海岸線處陸地側(cè)電場極值幾乎未受到斷裂帶電導(dǎo)率影響，而陸地—斷裂帶邊界處陸地側(cè)電場極值隨著斷裂帶電導(dǎo)率的減小而降低，陸地—斷裂帶邊界處斷裂帶側(cè)電場隨著斷裂帶電導(dǎo)率的減小而升高，該結(jié)果表明隨著季節(jié)因素的變化，各邊界電場會有所不同。

表1 邊界處電場極值變化Table 1 Extreme value changes of the electric field at the boundary

3 日東線管道時空特性分析

3.1 數(shù)據(jù)選擇

根據(jù)國家地磁臺站提供的從2016-10-12 22:12:00開始至2016-10-15 08:50:00 結(jié)束的地磁場數(shù)據(jù)，檢測到其Kp指數(shù)達(dá)到6，選取此段時間近20 萬個磁場水平分量和垂直分量Hx和Hy每秒采樣數(shù)據(jù)作為我們研究的磁擾動源(兩條虛線時刻之間)，如圖4 所示。應(yīng)用PFFEM方法，為符合精度要求，將磁場數(shù)據(jù)分解為56個主要幅頻分量，分別做有限元計算，取出每次磁場分量的電場結(jié)果疊加得到電場時空分布，輸入管道模型計算此次磁暴期間的管道PSP的時空分布。

圖4 2016-10-11 至16 日磁場分量Fig.4 Component of the magnetic field from 2016-10-11 to 16th (a)Data of Hx(b) Data of Hy

3.2 電性結(jié)構(gòu)建模

山東位于中國東部沿海，油氣管網(wǎng)密集，日東線原油管道于2013 年1 月開始投入運(yùn)行，其東始于山東日照的嵐山港，向西止于山東菏澤的東明縣，全長426 km，目前有日照、兗州、平邑3 個陰極保護(hù)站。日東線所經(jīng)大地巖石圈范圍內(nèi)包含郯廬斷裂帶，斷裂帶中段呈現(xiàn)低阻特性，研究區(qū)域見圖5。

圖5 華東沿海地質(zhì)及管道俯瞰圖Fig.5 Geological and pipeline survey of east China coast

為構(gòu)造較為真實(shí)的地質(zhì)三維電導(dǎo)率模型，根據(jù)華北地區(qū)的巖石圈結(jié)構(gòu)及電磁測深數(shù)據(jù)資料[22-24]，將電性結(jié)構(gòu)劃分為華北板塊、蘇魯造山帶、蘇北盆地、濟(jì)陽垇陷、膠北隆起、斷裂帶、膠萊盆地和黃海8 個區(qū)域。模型表面分塊如圖6 所示，黑色粗折線為日東線管道。設(shè)定x方向?yàn)橄蚰?，y正方向向東，z方向?yàn)橄蛏?，選取東西范圍880 km，南北范圍450 km，地質(zhì)層厚度為500 km。模型各區(qū)域分層分塊電導(dǎo)率參數(shù)見表2。

表2 3D電導(dǎo)率模型結(jié)構(gòu)參數(shù)Table 2 Structural parameters of 3D conductivity model

圖6 三維模型地表面分塊示意圖Fig.6 Blocks diagram of 3D model surface

3.3 電場分布特征

為模擬山東沿海范圍內(nèi)的在地磁暴時的電磁感應(yīng)環(huán)境，對模型上表面施加頻率為0.003 Hz、幅值為1000 nT的磁場，三維電導(dǎo)率模型及網(wǎng)格劃分如圖7 所示，經(jīng)過地磁有限元計算后可得到電場分布。

圖7 三維模型分層分塊及網(wǎng)格劃分Fig.7 3D conductivity model hierarchical block and grid generation

大地表面的地電場如圖8 所示，電場幅值在海岸處產(chǎn)生突變，位于海岸和斷裂帶之間的蘇魯造山帶電場較強(qiáng)，在海岸拐角處幅值E最大可達(dá)6 V/km。圖9(a)為模型上表面橫跨海陸的中線處(圖7 中黑色虛線)電場分布曲線，由于華北板塊與蘇魯造山帶之間的郯廬斷裂帶區(qū)域電阻率較低，顯示該區(qū)域電場減弱；由于海岸與斷裂帶之間的蘇魯造山帶受到兩側(cè)海岸-斷裂帶復(fù)合效應(yīng)的提升作用，電場增強(qiáng)。圖9(b)為日東線管道沿線處x、y方向地電場分量Ex和Ey，在海岸-斷裂帶復(fù)合效應(yīng)影響下均有提升，尤其是東向分量Ey影響顯著。

圖8 三維模型地表面電場幅值分布Fig.8 Electric field distribution on 3D model surface

圖9 電場分布Fig.9 Electric field distribution

3.4 管道PSP時空特性

考慮管道沿線大地電導(dǎo)率和管道走向兩大重要因素建立擴(kuò)展傳輸線LZS-DSPL管道模型，將經(jīng)過傅里葉逆變換合成后的管道電場時空分布輸入管道網(wǎng)絡(luò)，對2016 年10 月12-15 日磁擾動期間日東線管道PSP進(jìn)行計算，分別對時間分布和空間分布特征研究如下。

3.4.1 管道PSP空間分布特征

選取從2016-10-12 22:12:00 開始至2016-10-15 08:50:00 結(jié)束的地磁場數(shù)據(jù)可得到這段時間內(nèi)的電場E分布和管道PSP分布，經(jīng)觀察不同時刻的PSP空間分布，發(fā)現(xiàn)存在兩種相反的變化趨勢。因此，選取計算時間區(qū)間中800 min(10-13 11:32:00)和1500 min(10-13 23:12:00)2 個時刻管道電場E和PSP的計算值，進(jìn)行對照研究，如圖10 所示。

圖10 管道沿線PSP和電場空間分布Fig.10 Pipeline PSP and electric field distribution along pipeline

圖10(a)和(b)表示在本文所選數(shù)據(jù)計算下的兩時刻的管地電位PSP的空間分布，圖10(c)和(d)表示兩時刻的管道電場E的空間分布。日照站距離海岸僅5 km左右，處于海岸—斷裂帶之間區(qū)域，而兗州站在管道165 km處，特將這兩處陰保站位置在圖10 中標(biāo)示，并在圖中用豎虛線表示斷裂帶所在位置，橫虛線表示PSP零點(diǎn)。

T=800 min和T=1500 min管 道PSP的 空 間 分 布代表了2 種計算結(jié)果顯示分布趨勢，圖10(a)中，T=800 min的管道PSP在管道始端數(shù)值大，管道尾端降低至零值以下；圖10(b)中，T=1500 min的管道PSP分布在管道始端數(shù)值小，管道尾端上升至零值以上。這兩種相反的PSP分布趨勢分別對應(yīng)著不斷變換的負(fù)值管道電場和正值管道電場。

另 外， 在 圖10(a)和 圖10(b)中， 管 道 在A點(diǎn)(33 km)，B點(diǎn)(157 km)，C點(diǎn)(205 km)，D點(diǎn)(239 km)和F點(diǎn)(355 km)處 轉(zhuǎn) 折， 而E點(diǎn)(328 km)和G點(diǎn)(359 km)分別是管道經(jīng)過斷裂帶的兩側(cè)臨界點(diǎn)。圖10(a)和圖10(b)顯示管道PSP在A,B,C,D點(diǎn)處出現(xiàn)“拐點(diǎn)”，而在圖10(c)和圖10(d)中由于A,B,C,D點(diǎn)處于電導(dǎo)率較為均勻的區(qū)域，管道電場分布較為一致，由此可見管道PSP轉(zhuǎn)折程度與管道鋪設(shè)角度有關(guān)。同理，分布在斷裂帶之間的F點(diǎn)附近電場均勻，而PSP在F點(diǎn)處發(fā)生突變說明受到管道角度因素的影響。在E點(diǎn)和G點(diǎn)處管道角度不變，但在這兩點(diǎn)處管道PSP的變化顯然與斷裂帶電場的突變有關(guān)。

3.4.2 管道PSP時間分布特征

日東線兗州站和日照站安裝有PSP監(jiān)測裝置，分別將這2 處陰保站的管地電位PSP計算值與監(jiān)測值對比如圖11 所示。

圖11 日照與兗州的PSP計算值和測量值對比Fig.11 Comparison of the monitoring PSP in Rizhao with Yanzhou

在圖11(a)中，將兗州站PSP計算值與監(jiān)測值對比，發(fā)現(xiàn)變化基本一致，都在-1.3 V上下浮動，管地電位的最小值出現(xiàn)第1500 min左右，PSP計算值最小-1.68 V，測量值最小-1.65 V，PSP計算值最大出現(xiàn)在790 min達(dá)到-1.05 V，測量值最大為-1.03 V；在圖11(b)中，將日照站PSP計算值與監(jiān)測值對比，PSP實(shí)際值波動較大，與計算值相比有偏差，但出現(xiàn)突變的時間和趨勢基本相符，且計算值和測量值整體都在-1 V上下浮動。這在一定意義上驗(yàn)證了此模型的正確性。圖11(c)為日照、兗州兩站管道PSP的計算值對比，圖11(d)為兩站管道PSP實(shí)測值對比。計算值和測量值對比均顯示：日照處的PSP幅值高于兗州站PSP，日照站PSP幅值在-1 V上下波動，兗州站在-1.4 V上下波動；兩站波動趨勢一致，但是日照站PSP振幅更大，在此次磁暴最劇烈時日照PSP的波動幅度是兗州站的10 倍左右，極值接近(-3,+1) V。山東地區(qū)對管地電位PSP的管理要求是穩(wěn)定在-1.4至-0.85 V之間，可以看出日照站PSP超出控制范圍的218%。

這些結(jié)果說明，海岸和斷裂帶復(fù)合效應(yīng)在中石油日東線管道PSP產(chǎn)生了負(fù)面影響，在海岸—斷裂帶區(qū)域的日照站管道PSP幅值更大，波動更加劇烈，容易超出控制范圍，加速管道腐蝕。

因此，針對中國東部沿海管道的陰極保護(hù)工程應(yīng)該重點(diǎn)關(guān)注近海附近的日照站等，考慮海岸和斷裂帶復(fù)合效應(yīng)，建議采用雙極性恒電位儀加去耦合器的保護(hù)措施。

4 結(jié)論

本文從數(shù)學(xué)模型和實(shí)際數(shù)據(jù)兩方面驗(yàn)證海岸—斷裂帶復(fù)合效應(yīng)對沿海管道的影響，主要結(jié)論如下：

(1)海岸—斷裂帶復(fù)合效應(yīng)是使海岸—斷裂帶區(qū)域電場得到“提升”或“下拉”。

(2)兗州站的管道時變PSP計算結(jié)果與實(shí)測值在趨勢和幅值上基本保持一致，驗(yàn)證了管道PSP有限元算法的準(zhǔn)確性。

(3)在一次磁擾動期間日東線PSP空間分布特點(diǎn)表現(xiàn)為日照站PSP波動幅度是兗州站的10 倍左右，極值接近(-3,+1) V，超出控制范圍(-1.4,-0.85) V的218%，說明海岸和斷裂帶的復(fù)合效應(yīng)對沿海管道產(chǎn)生的負(fù)面影響很嚴(yán)重。建議采用雙極性恒電位儀加去耦合器的保護(hù)措施。

進(jìn)一步的，應(yīng)監(jiān)測沿海管道的地磁感應(yīng)電流，結(jié)合空間天氣的預(yù)報，提前找出災(zāi)變點(diǎn)，排除隱患。