中国东部沿海管道海岸—断裂带复合效应

2023-11-11方莉夏鹏程梁志珊熊树海毕武喜蓝卫

石油科学通报 2023年5期

方莉，夏鹏程，梁志珊*，熊树海，毕武喜，蓝卫

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幅值更大，波动更加剧烈，容易超出控制范围，加速管道腐蚀。

因此，针对中国东部沿海管道的阴极保护工程应该重点关注近海附近的日照站等，考虑海岸和断裂带复合效应，建议采用双极性恒电位仪加去耦合器的保护措施。