裂缝测井识别

(一) 成像测井识别裂缝

裂缝在井壁电成像和声成像测井图上均表现为连续或间断的深色条带，其形状取决于裂缝的产状。垂直缝和水平缝分别为竖直的和水平的条带，斜交缝为正弦波条带状，网状缝为正弦波条带状、竖直的和水平的条带的组合特征。低角度裂缝、高角度缝和网状缝3种裂缝性储层的电成像测井响应特征如图4－21所示。

对裂缝性质的解释主要要注意天然裂缝与层理、各种诱导裂缝，如钻具振动形成的裂缝、重泥浆造成的压裂缝、应力释放裂缝和井眼崩落的区别。特别是应力释放裂缝，既可在岩心上出现，也可在井壁上出现。在成像图上的特征为1组接近平行的高角度裂缝，且裂缝面十分规则。在常规测井解释中，容易误解释为低孔高角度裂缝型储层。当出现在岩心上时，很容易给岩心描述带来错觉，必须注意识别。其方法是看裂缝中有无泥浆侵入的痕迹，无侵入者为应力释放裂缝。应力释放裂缝只有1组，且裂缝面较为完整而压裂缝或为3组，或为一组不完整的，且仅出现在两个对称方向上的高角度裂缝。

根据X1井EMI图像显示特征，在EMI测量井段内，裂缝较发育，类型有:斜交缝、高角度缝、网状缝(表4－6)。

表4－6 X1井石炭系裂缝统计表

续表

图4－21 X1井裂缝类型EMI图像

斜交缝:倾角小于90°的开口缝，包括高角度斜交缝(倾角≥70°)、低角度斜交缝(10°≤倾角<70°)，EMI图像显示为黑色正弦曲线。EMI测量井段内发育有多组(数条)斜交缝。

高角度缝:倾角等于90°的开口缝，EMI图像上表现为黑色竖线，缝宽不等，通常情况下两条线相互平行，延伸较长。

网状缝:由两组以上产状不同的裂缝相互切割的呈网状的开口缝组成。本井石炭系网状缝相对欠发育，网状缝基本是与斜交缝交替发育。

(二) 常规测井识别裂缝

理论研究结果表明，深、浅双侧向电阻率的大小及差异性质除受流体性质影响外，还严重地受到另外两个因素的控制。一是裂缝张开度、裂缝密度、裂缝产状、裂缝径向延伸等裂缝自身的特征二是岩石本身的电阻率。

1. 双侧向测井

(1) 裂缝性储层在深、浅双侧向上的响应特征

由于深、浅双侧向电阻率的大小及差异性质受流体性质、裂缝张开度、裂缝密度、裂缝产状、裂缝径向延伸以及岩石本身的电阻率影响。因此高角度裂缝(>75°)为主的储层来说，深、浅双侧向出现正差异，且比值随裂缝倾角、裂缝张开度、裂缝径向延伸度、裂缝纵向穿层长度的增大而增大。对于低角度裂缝(<75°)，深浅侧向出现负差异。此外还必须考虑到岩块电阻率Rb的影响，即对同样的裂缝，Rb越高，深浅双侧向的电阻率差异也越大。斜角缝、高角度缝和网状缝三种裂缝性储层的深、浅双侧向测井响应特征分别如图4－22～图4－24。

图4－22 X2井斜角缝测井响应特征

(2) 裂缝性储层在微球形聚焦测井曲线上的响应特征

微球形聚焦测井具有比双侧向的径向探测深度浅，垂直分辨率高的特点，因此它受井眼和泥饼的影响比双侧向测井大，但它分辨裂缝的能力却远比双侧向强。因此，当井眼较规则时，微球形聚焦测井在裂缝段将发生比双侧向较多的起伏，且在双侧向电阻率背景上来回变化，如图4－23。

2. 密度测井

密度测井测量的是岩石的体积密度，主要反映的是岩石的总孔隙度，而与孔隙的几何形态无关。由于密度测井仪为极板推靠式仪器，当极板接触到天然裂缝时会对密度测井产生较大影响。

图4－23 X3井高角度缝测井响应特征

图4－24 X4井网状缝测井响应特征

3. 补偿中子

与密度测井类似，补偿中子测量的也是岩石的总孔隙度，不受孔隙几何形态和分布的影响。补偿中子由于其测井探测深度较大，而成为确定非均质的裂缝性火山岩油藏总孔隙度的有效方法。在裂缝性火山岩剖面层段上，补偿中子显示为相对高的孔隙度值，而且裂缝越发育，中子孔隙度就越大。与其他常规测井类似，补偿中子也同样只能指示裂缝带的位置，不能确定裂缝的发育方向。

4. 声波时差

裂缝在声波时差曲线上的反应与井筒周围裂缝的产状及发育程度有关。声波时差对高角度缝没有反应，对低角度缝或网状裂缝，声波时差将相应增大。当遇到大的水平裂缝或网状裂缝时，声波能量急剧衰减而产生周波跳跃现象。因此利用声波时差可以识别水平裂缝或网状裂缝，但不能用于识别垂直裂缝。

(三) 地层倾角测井识别裂缝

地层倾角测井是探测天然裂缝的各种方法中较为有效的方法之一。用地层倾角测井资料识别裂缝的方法有:裂缝识别测井、电导率异常检测、定向微电阻率、双井径曲线等。

1. 裂缝识别测井(FIL)

地层倾角的微电阻率曲线常在高阻背景上以低的电阻率异常显示出裂缝。FIL是利用地层倾角的4条微电阻率曲线，按顺序排列组合相邻两极板的4组重叠曲线(1－2，2－3，3－4，4－1)，裂缝则以明显的高电导率异常显示出来。当任一极板通过充满高电导率泥浆的裂缝时，其电导率升高，重叠曲线出现幅度差。一般高倾角裂缝常以一组或两组明显的幅度差出现，垂直裂缝在两条曲线上有较长井段的异常而水平裂缝在4条重叠曲线上均有较短的异常。这种方法的缺点是不能准确地识别沉积构造和裂缝。

2. 利用电导率异常检测识别火山岩裂缝

该方法是利用地层倾角测量的原始记录在曲线对比垂向移动范围所确定的井段上，求出各极板与相邻两个极板电导率的最小正差异值，并把此值叠加在该极板的方位曲线上。作为判别裂缝的标志，这种方法排除了由于层理所引起的电导率异常，突出了与裂缝有关的电导率异常。在电导率异常检测DCA成果图上，不仅可以直接显示出裂缝的存在，而且直接给出了裂缝存在的方位。用该方法必须满足下列3个条件:①电导率超过一定的水准，②电导率数值之差足够大，③异常可以在极少数连续层位上探测到。

3. 双井径重叠法

双井径重叠是识别裂缝的一种重要方法，通常具有较好的使用效果。根据地层倾角测井曲线显示的定向扩径、椭圆形井眼及相对方位角曲线平直无明显变化等，可以划分出高角度裂缝层段。而且，根据扩径方位或椭圆形井眼的长轴方向，可以确定高角度裂缝的方向。一般双井径曲线值与钻头直径均相等为硬地层，双井径曲线值均小于钻头直径为渗透层，双井径曲线值均大于钻头直径为泥岩或疏松易塌层，双井径曲线值之一大于钻头直径，另一等于或小于钻头直径，呈椭圆形井眼，为高角度裂缝。

(四) 利用双侧向测井资料定性识别裂缝的实现方法

为了能有效识别出裂缝、优化单井射孔层段，从而更好地指导现场生产工作，在基于对众多的成像测井资料与常规测井资料进行对比分析后，建立了天然裂缝的常规测井解释模型。这种方法不同于裂缝孔隙度计算，是一种定性的判断方法，其主要方法是首先提取成像测井资料中典型的裂缝，然后对常规测井资料进行标定，从而提取裂缝在常规测井资料中的响应特征，然后针对这些特征进行编写识别程序，从而使用计算机对裂缝进行自动识别评价。

通过分析笔者发现当地层出现网状缝或其他类型的斜交缝的时候，微球测井曲线一般会比较迅速地下穿双侧向，在非裂缝处微球一般会悬浮于双侧向之上。其原理为:井壁的张开裂缝会导致微球电阻率值(RXO)的急剧下降(张开裂缝中充满泥浆所致)，依此可以识别裂缝发育井段。通过与成像测井对比发现，该方法可以识别多数的张开裂缝，但无法区别钻井诱导裂缝。

其识别图版为:

BRXO≤0.8并且BXOT≤0.8

则该井段为裂缝发育井段当RT≤70和RXO≤70的情况不能使用本方法判断。

其中:BRXO=RXO/RXO1

BXOT=RXO/RT。

式中:RXO———冲洗带电阻率

RT———地层真电阻率

BRXO———冲洗带电阻率变化幅度

BXOT———冲洗带电阻率与地层真电阻率幅度比

RXO1———RXO曲线上当前深度点的上一个采样点。

通过BRXO和BRXT的斜率大小来判断裂缝的存在。

图4－25是利用该方法进行裂缝识别的一个实例，图中第三道裂缝指示曲线即是根据双侧向和微球形聚焦测井曲线计算得到的，它只是一条定性指示曲线。无量纲，代表该深度存在裂缝或裂缝发育的相对程度。

通过利用常规测井曲线计算判断的储层裂缝段与成像测井拾取的裂缝层段对比，认为利用常规测井资料判断裂缝的方法是可行的，也是较为有效的。

(五) 基于测井曲线元的裂缝定量识别

针对火山岩裂缝性油气藏裂缝测井识别这一难题，在充分分析其裂缝曲线元及其变化特征的基础上，刻画了裂缝曲线元的数学特征，建立了基于测井曲线元的裂缝概率模型，进而来计算裂缝发育的概率。

裂缝的存在对电性、放射性等各种物理性质均有不同程度的影响，其影响可在测井曲线元的变化形态上造成异常响应。由于各种测井方法对裂缝的敏感程度并非完全相同，加之某些非裂缝因素也可能引起与裂缝相同的异常响应。所以，用一、二种常规测井方法识别裂缝的准确性往往很低，在井眼条件较差的情况下尤其如此，而多种测井信息综合反映裂缝的可能性明显增大。因此，本节利用多种常规测井信息来建立基于测井曲线元的裂缝概率模型，进而来对研究工区的裂缝进行定量识别。

1. 裂缝曲线元及其特征

在裂缝发育段，三侧向电阻率曲线和微电阻率都比上下相邻曲线段读值降低，但不同的层段降低的程度有所不同。也就是说，同样是裂缝发育段，曲线的形态还与岩性、层厚、泥浆电阻率、侵入深度等因素有关。而这些因素的影响都反映在曲线元的形态变化上。

图4－25 X5井利用双侧向和微球形聚焦测井识别裂缝实例

为了便于准确地刻画测井曲线的变化形态，引入了测井曲线元的概念。在测井曲线上，如果对曲线所考查的某一性质与邻近的曲线段明显不同，则把这样的一段曲线称为测井曲线元，简称为曲线元。记为

准噶尔盆地火山岩储层测井评价技术

或者记为C∶C∈［a，b］。通常a，b为测井曲线的左右刻度。

假设Ci－1∈［a，b］、Ci∈［a，b］、Ci+1∈［a，b］分别为相邻的3段曲线元，假设Ai－1，Ai，Ai+1表示相应曲线元的某一性质，ε1为一给定值，依据定义则有

准噶尔盆地火山岩储层测井评价技术

式(4－13)中，F并不表示一种单纯的映射，而是一种刻度，就是一种度量相邻曲线元某一特征的差别的方式。ε1是针对某一项待考察的指标给定的限定值，也就是划分不同曲线元的截至值。

若已知Ci∈［a，b］为一曲线元，［a，b］为该曲线的左右刻度。有时这一区间也可限定在该段曲线的最大最小值之间。

曲线元的数理统计分析主要计算曲线元的均值μ、极差J、数学期望E、方差σ2或标准差σ等。在进行数理统计分析之前，先要有一个合理的假定条件，对于xi∈［a，b］，i=0，1，…，n取任意值的几率都是相等的。因为对于某一段地层来说，在已知的值域内(比如曲线的左右刻度)，没有任何理由让某一项测井量(例如电阻率)只取某一值而不能为另外的值，也就是说，在值域测井量取任一个值的机会是均等的。因此对于任意点的概率Pi有

准噶尔盆地火山岩储层测井评价技术

在这一假设条件下，测井曲线元的数学期望和方差就可以计算了。

根据关键井的岩心标定，裂缝在测井曲线上的变化特征主要表现为三侧向电阻率曲线出现高值背景上的降低，深浅电阻率的幅度差也有所减小，同时微电极曲线也表现为同样的特征。自然伽马曲线没有增大或增大很小，通过滤波可以消除。把这样的曲线段称为裂缝曲线元。

2. 火山岩储层裂缝指标的定义

在实际处理过程中，考虑到火山岩地层岩性的复杂性，定义的裂缝指标有如下4种方式。

(1) 双侧向或双感应幅度差

直接在综合测井曲线图(对数坐标)上找到致密段和裂缝段的双侧向(或双感应)幅度差绝对值，ΔRb和ΔRf。当前处理深度的电阻率幅差指示的裂缝概率为:

准噶尔盆地火山岩储层测井评价技术

式中:P———裂缝概率

ΔR———当前深度的电阻率对数的幅度差绝对值。

(2) 井径测井曲线

准噶尔盆地火山岩储层测井评价技术

式中:CAL———当前深度的井径

CALf和CALb———裂缝层段和致密层段的井径。

(3) 微球形聚焦测井

对微球形聚焦电阻率测井曲线的对数lg(x)进行滤波处理，得到滤波后的测井曲线lg(x)'，提取剩余变化(Dx)，则裂缝概率:

准噶尔盆地火山岩储层测井评价技术

式中:Dxf和Dxb分别为裂缝层段和致密层段的剩余变化值Dx=lg(x)－lg(x)'。

(4) 其他曲线

对其他非电阻率测井曲线x进行滤波处理，得到滤波后的测井曲线x'，提取剩余变化Δx=x－x'。

准噶尔盆地火山岩储层测井评价技术

式中:Δxf和Δxb———裂缝层段和致密层段的剩余变化值。

上述指标除了双侧向取绝对值外，均考虑到裂缝特征在曲线上的方向性。

3. 裂缝识别

一般测井采样间距为12.5cm，岩心裂缝观察表明裂缝的一般长度在1个采样间距到几十个采样间距之间。因此，从采样间距上考虑，如果要利用常规测井曲线识别裂缝，必须至少有3个采样点，即2个采样间距构成的裂缝。因为如果是2个采样点，2个采样点读值的变化只可能是由大变小、相反或不变，3种之一，而不能构成裂缝曲线元由大变小再变大的变化特征。这是根据常规测井曲线判断裂缝存在的必要条件。一条裂缝的延伸长度必须至少大于12.5cm×2cm才能够被常规测井曲线识别到。从测井解释的角度说，常规测井资料识别的裂缝长度至少是25cm。

根据裂缝曲线元的特征编制了基于测井曲线元的计算机裂缝识别软件系统，其识别过程如下。

1) 测井数据录入。

2) 判断原始数据文件是否有三侧向曲线和微电极曲线，或二者之一。如果都没有，则无法进行裂缝识别，退出系统。

3) 测井资料校正和数据标准化。资料校正是正确识别的前提，数据标准化便于进行曲线元拟合和计算曲线元的数字特征。该软件主要采用了极差标准化、极差正规化和标准差标准化3种方法。

4) 读入分层数据。研究中只对砂岩储层段进行了裂缝识别，对泥岩未处理，所以在识别之前先要读入分层数据。如果该井还没有分层，必须先分层。

5) 曲线元滤波。滤波主要是消除曲线上微小的扰动，因为它会影响到对裂缝曲线元识别。滤波方法多采用加权滑动平均法，如钟型函数或汉明函数等，也可采用卡尔曼滤波。

6) 从第1层开始，逐点判断电阻率与邻近上下两个采样点关系，设存在xi－1，xi，xi+1为3个相邻的采样点，并且给定ε为一门限值，如果式(4－19)成立，则从采样点xi开始记录采样点数S1，直到式(4－19)成立，则必然有式(4－20)成立。同样开始记录采样点数S2，直到式(4－20)不成立

准噶尔盆地火山岩储层测井评价技术

记录满足式(4－19)和式(4－20)总采样点数(S=S1+S2)。对三侧向和微电极曲线依据式(4－19)和式(4－20)作判断。

7) 对同一层段的GR曲线进行判断，消除由于泥岩夹层引起电阻率降低而误判为裂缝的层段。

8) 计算裂缝存在概率Pf。

计算裂缝曲线元的极差(J)、数学期望(E)、方差(D)。因为裂缝曲线元的形态特征是电阻率在高值背景上的骤然降低，表明其极差很大，并且降低越明显，极差就越大。同时，其数学期望与极差的差值也随之增大。方差越大，裂缝存在的可能性越大，因此裂缝的概率与方差成正比。此外裂缝曲线元的曲线突变不会延续很长，否则这种突变成了一种渐变，也就不是裂缝了。定义单条曲线判断裂缝存在的概率计算式为

准噶尔盆地火山岩储层测井评价技术

在识别过程中，对于JD－581测井系列同时采用了深、浅三侧向及其幅度差、微电极及其幅度差、GR曲线、感应曲线和声波曲线对于CLS3700测井系列则将三侧向换为双侧向电阻率曲线。根据式(4－21)，每1条曲线都将给出1个单曲线裂缝存在概率，按照贝叶斯准则，所有曲线指示裂缝存在概率由下式计算

准噶尔盆地火山岩储层测井评价技术

式中，Pf，i———单曲线裂缝存在概率

l———参加裂缝判别的曲线条数。同时给定了一个经验参数εP作为裂缝存在概率的截止值，规定只有Pf>εP时，才认为该裂缝存在，并且被记录下来。

经本书研究确定，RXO权值变化范围为0.4～0.5、Rt和Ri幅度差的权值变化范围为0.25～0.35、井径变化范围为0.1～0.2、声波或密度变化范围为0.15～0.25。

4. 裂缝识别实例

基于上述方法，对研究盆地内的白X1井进行裂缝识别，其识别成果图如图4－26所示。由该图可知，在1717～1728m井段所计算的裂缝概率值较高，说明该段发育裂缝的概率较大而在1710～1717、1728～1750m井段处，所计算的裂缝概率值较低，说明该段发育裂缝的概率较小。裂缝概率反映地层中存在裂缝的概率大小，是对裂缝发育程度的判别。概率值越大，裂缝越发育反之，概率值越小，裂缝发育越差。从FMI图像上看，1715.0～1720.0m发育一组高角度缝，在1720.0～1722.0m发育一组高角度缝以及斜交缝，在1722.0～1728.8m发育一组雁状缝以及斜交缝。根据岩心描述裂缝统计(表4－7)可知，在该井段内，裂缝较发育，类型有:斜交缝、直劈缝、网状缝、充填缝以及雁状缝。由此可知，对缺乏成像测井和岩心描述等资料的情况下，该方法能够利用常规测井资料来较准确地识别其裂缝发育的井段。

图4－26 白X1井裂缝识别实例

表4－7 白X1井裂缝统计表

由于本区石炭系地层火山岩岩性复杂，裂缝的常规测井响应特征(如声波时差、中子孔隙度、深浅双侧向等)受岩性影响较大，容易将岩性的变化混淆为裂缝。而成像测井图可直观地反映裂缝的形状(如弯曲程度)、填充状况。从本区岩心资料和成像测井资料综合来看，利用成像测井来识别裂缝较为有效。

卡西欧5800计算器坐标的正反算程序（曲线），程序在计算器中约占2900字节，具体程序写法如下：

主程序 坐标正反算

程序名称：ZBZFS

Lb1 A↙

30→Dim Z :"0=ZS,1=FS"?Z ↙

If Z=0:Then Goto B:IfEnd↙ (Z=0进入里程点坐标正算）

If Z=1:Then Prog"FSLCZ "↙ (Z=1进入反算里程边桩）

Lb1 B ↙

Prog”ZBZS” ↙

子程序反算里程桩 名称：FSLCZ

Lb1 1↙

"XK="?X:"YK="?Y↙（输入任意测点的X、Y坐标）

“K=”?K ↙（试算里程，任意输入标段内里程点即可，也可以输入大致的估算里程加快速度）

Lb1 2↙

Prog"SJK ":T-90 →W:W<0=>360+W→W:Abs((Y-Z[19])*Cos(W)-(X-Z[18])*Sin(W))→S↙

If S<0.0001:Then Goto4:Else Goto 3:Ifend↙

Lb1 3↙

K+S→K: Prog"SJK":T-90 →W:W<0=>360+W→W:Abs((Y-Z[19])*Cos W-(X-Z[18])*Sin (W)→Q↙

If Q<0.0001 :Then Q→S: Goto 4:Else ifQ<S: Then K+Q→K: Goto 2 :Else if Q>S :Then K-Q→K:Goto 2:Ifend:Ifend:Ifend↙

Lb1 4↙

Pol (X-Z[18],Y-Z[19]:

"DP(-Z+Y)=":I◢(偏距）

"K=":K+S→K◢（里程）

Goto 1↙

子程序坐标正算 名称：ZBZS

Lb1 0 ↙

“XHS="?G(后视点X):"YHS="?L(后视点Y):"XZJ="?M(置镜点X):"YZJ="?N(置镜点Y):Pol(G-M,L-N):"DH=":I(后视距)◢J<0=>J+360→J:"FH=":JDMS◢(后视方位角)

Lb1 1↙

“K=”?K ：（输入所需计算里程）Prog"SJK "↙

XI : Z[18]↙

YI: Z[19] ↙

Pol(Z[18]-M,Z[19]-N): J<0=>J+360→J↙

“PJ=”?P↙ (输入桩与线路夹角)

“PD=”?D↙ （输入桩距中线的距离）

Z[18]+D*Cos(T+P) →Z[20] ↙

Z[19]+D*Sin(T+P) →Z[21] ↙

“X=”: Z[20] ◢ （放样坐标X）

“Y=”: Z[21] ◢ （放样坐标Y）

Pol(Z[20]-M,Z[21]-N):"D=":I◢（放样距）J<0=>J+360→J:"F=":JDMS◢（放样方位角）

Goto 1↙

子程序数据库 名称：SJK

if K<本段曲线终点里程 And K≥上段曲线终点里程 ：Then 本段曲线终点里程→Z[1] : 上段曲线终点里程（第一段曲线输起点的里程）→Z[2] ：1→O （注：左偏曲线输入-1→O,右偏曲线输入1→O）: 半径→R : 曲线偏角→A:第一缓和曲线→Z[6] : 第二缓和曲线→Z[7] : 交点 X→B :交点 Y→C : 小里程向交点方位角→E : 交点向大里程方位角→F : Prog”JSPB”:Return: Ifend↙

if…………Prog”JSPB”:Return:Ifend（曲线段分段输入）↙

补充直线段输入如下 (只需输线路的最后一段直线数据）

if K<本段直线终点里程 AndK≥本段直线起点里程:Then 本段直线终点里程→Z[3]:终点X→Z[16]:终点Y→Z[17]:方位角→E: Z[16]+ (K- Z[3])*Cos(E)→Z[18]:Z[17]+ (K- Z[3])*Sin(E)→Z[19] : Return:Ifend↙

子程序计算判别 名称：JSPB

Lb1 2 ↙（曲线要素计算）

Z[6]/2- Z[6]^3/(240*R^2)+ Z[6]^5/(34560*R^4) →Z[8] ↙ （M1）

Z[7]/2- Z[7]^3/(240*R^2)+ Z[7]^5/(34560*R^4) →Z[9] ↙ （M2）

Z[6]^2/(24*R)- Z[6]^4/(2688*R^3) →Z[10] ↙（P1）

Z[7]^2/(24*R)- Z[7]^4/(2688*R^3) →Z[11] ↙（P2）

π*A*R/180+0.5*( Z[6]+ Z[7])→Z[25] ↙（曲线总长）

90* Z[6]/(R*π) →Z[14] ↙ （第一缓和曲线总偏角）

90* Z[7]/(R*π) →Z[15] ↙（第二缓和曲线总偏角,可以省略）

Z[8]＋（R+Z[10])Tan(A/2)-(Z[10]-Z[11] )/Sin( A)→Z[12]↙(切线T1)

Z[9]＋（R+Z[11])Tan(A/2)+(Z[10]-Z[11] )/Sin (A）→Z[13]↙(切线T2)

B+ Z[12]*Cos (E+180)→ Z[16] ↙ （ZH点X）

C+ Z[12]*Sin(E+180)→ Z[17] ↙（ZH点Y）

Z[1]- Z[25]→Z[3] ↙(ZH点里程)

Z[3]+ Z[6]→Z[4] ↙(HY点里程)

Z[1]- Z[7]→Z[5] ↙(YH点里程)

Goto 3 ↙

LB1 3 ↙(判断里程点与曲线关系)

if K≤Z[3] AndK>Z[2] : Then Goto 4 : Ifend ↙

if K≤Z[4] AndK>Z[3] : Then Goto 5 : Ifend↙

if K≤Z[5] AndK>Z[4] : Then Goto 6 : Ifend ↙

if K≤Z[1] And K>Z[5] : Then Goto7 : Ifend ↙

Lb1 4 ↙（里程小于直缓点直线独立坐标）

K- Z[3] →Z[23] : 0→Z[24] : E→T : Goto8↙

Lb1 5 ↙（第一缓和曲线独立坐标）

K- Z[3] →H ↙

H-H^5/(40*R^2* Z[6]^2)+H^9/(3456*R^4* Z[6]^4) →Z[23] ↙

H^3/(6*R* Z[6])-H^7/(336*R^3* Z[6]^3) →Z[24] ↙

90*H^2/( R*π* Z[6]) →T ↙

if O>0 ：Then T +E→T : Eles E-T →T : T<0=>360+T→T : Ifend ↙

Goto 8 ↙

Lb1 6 ↙（圆曲线独立坐标）

K- Z[4] →H ↙

H*180/( R*π)+ Z[14]→T ↙

R*Sin(T)+ Z[8]→Z[23] ↙

R*(1-Cos(T))+ Z[10]→Z[24] ↙

if O>0 ：Then T +E→T : Eles E-T →T : T<0=>360+T→T : Ifend ↙

Goto 8 ↙

Lb1 7 ↙（第二缓和曲线独立坐标）

Z[1] -K →H ↙

H-H^5/(40*R^2* Z[7]^2)+H^9/(3456*R^4* Z[7]^4) →U↙

H^3/(6*R* Z[7])-H^7/(336*R^3* Z[7]^3) →V ↙

90*H^2/( R*π* Z[7]) →T ↙

Z[13]Cos(A)+ Z[12]-U*Cos(A)-V*Sin(A)→Z[23] ↙

Z[13]*Sin(A)-U*Sin(A)+V*Cos(A)→Z[24]↙

if O>0 ：Then F-T→T : T<0=>360+T→T : Else F+T →T : Ifend ↙

Goto 8 ↙

Lb18

if O<0 : Then- Z[24]→Z[24] : Ifend ↙

Z[16]+Z[23]*Cos(E) -Z[24]*Sin（E）→Z[18] ↙

Z[17]+Z[23]*Sin(E)+Z[24]*Cos(E)→Z[19] ↙

Return↙

function

fupin(w,c,l,r)

d1=4*l*c^4*r

d2=w.^2

d3=w.^3

d4=(2*c*r)*(-2*l*c+1./d2)

d5=(-4*l*c)./w+1./d3

d6=sqrt(d4.^2+d5.^2)

p=d1./d6

plot(w,p)

试试吧，你写一个函数，把里面的参数提出来，这样你在计算其他曲线的时候，就改相应的值就好了！不过我画出来感觉不太像，你按照这个思路试试！

---