中國東部沿海管道海岸—斷裂帶復合效應

方莉，夏鵬程，梁志珊*，熊樹海，畢武喜，藍衛

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

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

0 引言

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

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

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

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

1 三維埋地管道PFFEM 理論

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

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

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

(1)大地磁場FFT分解

將地磁水平分量和垂直分量分別作為大地模型的初始磁源，那么首先對地磁臺磁場數據東向和北向分量hx、hy進行分解為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，因此可看作磁準靜態場，與渦流電流密度相比，位移電流可以略去不計[20]。一次電流源Jp產生了一次場(一次電場Ep和一次磁場Hp)，一次場在大地導體中產生二次場(二次電場Es和二次磁場Hs)，內部無源。時諧磁場中以e-jωt為時諧因子，其中j為虛數單位。則二次電磁場滿足

式中，假設介質的磁導率為真空磁導率μ0=4π×10-7A/m，介質的介電常數為真空介電常數ε0=8.85×10-12F/m，γ表示導體電導率。

由赫姆霍茲定理可知導電區電場和磁場有唯一解[21]，在渦流場求解控制方程，根據區域結構和邊界條件，需要考慮計算機運行時間和計算精度選擇適當的電磁位對，常見方法有A法和T法。T法引入位函數T和Ω，T為矢量電位，Ω為標量磁位，求解過程未知數相比于A法較少，且能達到滿意精度。根據麥克斯韋方程，J、H與T、Ω之間的關系是

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

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

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

(3)管道電場IFFT合成

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

在管道j點處的x、y電場分量值由磁場x、y分量產生的電場分量對應相加得到

則t時刻，管道j點處的電場幅值為

(4)管道PSP計算

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

2 海岸—斷裂帶復合效應

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

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

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

圖3 不同斷裂帶電導率下的海岸—斷裂帶效應電場Fig.3 Electric field effect of coastal-fault zone under different conductivity

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

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

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

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

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

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

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

3.1 數據選擇

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

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

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

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

表2 3D電導率模型結構參數Table 2 Structural parameters of 3D conductivity model

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

3.3 電場分布特征

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

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

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

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

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

3.4 管道PSP時空特性

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

3.4.1 管道PSP空間分布特征

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

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

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

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

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

3.4.2 管道PSP時間分布特征

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

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

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

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

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

4 結論

本文從數學模型和實際數據兩方面驗證海岸—斷裂帶復合效應對沿海管道的影響，主要結論如下：

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

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

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

進一步的，應監測沿海管道的地磁感應電流，結合空間天氣的預報，提前找出災變點，排除隱患。