[0001] 本发明涉及油气工程技术领域，具体涉及一种从DAS数据提取速度谱的方法、可读存储介质和电子设备。

[0002] DAS技术即分布式声波传感技术，是一种利用光纤瑞利后向散射光信息实现光纤沿线声波探测的技术。DAS技术主要由解调仪和光纤两部分构成。在油气工程领域，可以利用DAS技术获取油井的声波数据。应用时，将光纤与油井的套管耦合，解调仪不断地向连接的光纤发射激光脉冲，当光纤的某个位置受到井筒内流体的扰动而产生应变时，该位置处
的光纤折射率会发生改变，导致后向散射光的相位也会发生改变，由于瑞利后向散射光相
位变化与作用于光纤上的声波呈线性关系，通过解调光相位变化即可线性获取声波数据。
由DAS技术获取的声波数据即分布式声波传感数据，也即DAS数据。

[0003] 从DAS数据可以提取声速数据，现有技术中，首先对油井全井段的DAS数据进行数据块划分，再对每一数据块进行频率‑波数变换，计算得到各频率‑波数图。在频率‑波数图中，会出现大致呈线性分布的声速条带，由工作人员在声速条带内沿其大致倾斜方向手动
绘制斜率线，该斜率线的斜率即代表该数据块对应的声速数据。对于划分的所有数据块，均需进行上述手动操作，根据每一数据块的声速数据，拼接得到全井段静态速度谱，并由静态速度谱进一步获取声波的上行波速度谱与下行波速度谱。现有方法，声速数据获取的过程
重复繁琐，获取效率低，且由工作人员目测频率‑波数图中的声速条带，并大致绘制斜率线，由斜率线的斜率值得到的声速数据也不够客观准确，不能真实精确地反映油井中的声速数
据。

[0049] 下面将结合附图进一步说明本发明实施例。在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明的简化描述，而不是指示或暗示所指的装置或组件必需具有特定的方位、以特定的方位构造和操作，因此不能理解为对
本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。其中，术语“第一位置”和“第二位置”为两个不同的位置。

[0050] 图1为本发明一个实施例所述的从DAS数据提取速度谱的方法流程图。如图1所示，从DAS数据提取速度谱的方法包括如下步骤：

[0051] S101:将全井段DAS数据按照设定深度间隔划分为多个数据块，获取每一所述数据块对应的频率‑波数图。

[0052] 利用分布式光纤传感器获取油井的全井段DAS数据，图2为分布式光纤传感器布设示意图，如图2所示，油管11外套设套管12，将一根单模光纤13与套管12耦合，沿套管12延伸方向布设在套管12外，单模光纤13在井口的一端连接解调仪14，解调仪14既向光纤13发射
激光脉冲，同时又接收光纤13传回的光相位变化数据，并进一步将光相位变化数据转换为
声波数据，即全井段DAS数据。全井段DAS数据过于庞大，根据经验设定深度间隔，将全井段DAS数据按照深度间隔划分为多个数据块，例如，可以将深度间隔设定为200m。

[0053] 原始的DAS数据为时间‑空间域(T‑X)数据，经过傅里叶变换，将每一数据块分别变换到频率‑波数域(F‑K)。图3为本发明一个实施例所述的频率‑波数图，如图3所示，横坐标为波数，纵坐标为频率，背景为功率谱，也即傅里叶变换振幅的平方值，频率‑波数图中的亮点表示在特定频率和波数组合下能量较强。根据波数、频率与声波速度之间的函数关系，频率‑波数图上过原点画任意直线，则该直线的斜率值即代表一个声速值。因傅里叶转换过程中，会产生共轭结果，因此，波数轴上存在负值，不影响后续数据处理。

[0054] S102:对每一所述频率‑波数图进行网格划分。

[0055] 在波数坐标轴和频率坐标轴上，分别按照一定的间隔将频率‑波数图划分为依次设置的网格。网格的划分间隔可以结合频率‑波数的规格和经验等人为设定，不对最终获取的速度谱数据产生决定性的影响。

[0056] S103:获取每一网格内的功率谱。

[0057] 在将每一数据块经傅里叶变换，转换到频率‑波数域时，每一网格范围内都对应有功率谱，功率谱即傅里叶变换振幅的平方值。

[0058] S104:过所述频率‑波数图的原点设定多条不同斜率值的斜率线，每一所述斜率值均在预设的斜率值范围内。

[0059] 斜率值即声波在井筒中静态流体内的传播速度，预设的斜率值范围，根据声波在井筒中的不同流体中的传播速度进行限定。油井中还存在水和气体，声波在气体中的传播
速度是340m/s，在纯油中的传播速度是1050m/s，在纯水中的传播速度是1525m/s，在气水油混合流体中的速度则介于340‑1525m/s之间。从而，在频率‑波数图中，过原点设定一系列斜率线，使得所有斜率线的斜率值处于上述范围内。

[0060] S105:针对每一所述频率‑波数图，计算设定的每条所述斜率线的积分值，所述积分值为所述斜率线在穿过的每一所述网格内的长度值与该所述网格内的功率谱的加权平
均值。

[0061] 不同斜率值的斜率线穿过频率‑波数图中的不同位置的多个网格，穿过的每一网格内均对应特定长度的斜率线，且均有对应的功率谱，将一条斜率线穿过的每一网格的长
度值与对应的功率谱求积值，将每一网格对应的积值相加，再与斜率线的总长度求商，即该斜率线对应的积分值。

[0062] S106:以最大积分值对应的所述斜率线的斜率值作为其所在的所述数据块对应的静态速度值。

[0063] 积分值越大，则表明其对应的斜率线的斜率值越接近所在深度段的真实的声波速度。将每一斜率线对应的积分值之间进行比较，从中选出积分值的最大值，则该最大积分值对应的斜率线，其斜率值即该数据块对应的静态速度值。

[0064] S107:将每一所述数据块对应的所述静态速度值进行拼接，得到所述全井段DAS数据对应的静态速度谱。

[0065] 原始DAS数据按照设定深度间隔被划分为了多个数据块，获得每一数据块对应的静态速度值后，将所有不同深度段的静态速度值进行拼接，才能够得到油井全井段深度范
围内的静态速度谱。

[0066] S108:根据所述静态速度谱和井筒内的流体速度，获取井筒内的上行波速度谱与下行波速度谱。

[0067] 静态速度谱是声波在流体介质处于静止状态下的声波速度，在油井内，油水混合流体向地面方向输送，则根据上行波和下行波与流体流动方向之间的关系，结合静态速度
谱，可以计算得到井筒内的上行波速度谱与下行波速度谱。

[0068] 采用本发明从DAS数据提取速度谱的方法，可以用软件算法从DAS(分布式声波传感)数据中自动提取声波速度谱，自动化过程大大减少了人工干预，提高了数据处理效率，且在一定声波速度范围内设定一系列不同斜率值的斜率线(也即设定一系列不同声波速度
值)，并从中选取积分值最大的斜率线对应的速度值作为该井段的最终静态速度值，更加接近真实速度，提高了声速数据的准确度。另外，使用不同斜率线的方式，增加了方向灵活性，方向灵活性是指计算积分值时按照斜率线的方向进行积分，而不是整张频率‑波数图整体
范围内积分，从而提高了计算效率，减少了内存消耗。

[0069] 在将每一数据块分别进行傅里叶变换，获取对应的频率‑波数图的步骤之前，还可以包括如下步骤：

[0070] S201：对每一数据块进行数据预处理。

[0071] 原始的DAS数据中，可能存在异常值，对各数据块进行数据预处理，包括但不限于滤波、降噪处理等，以消除噪音和干扰，保证数据的有效性，保证最终结果的可靠性。

[0072] 可选地，所述对每一所述频率‑波数图进行网格划分的步骤中，具体包括：

[0073] S301:获取所述频率‑波数图的像素点；

[0074] S302:以所述频率‑波数图的原点为起点，沿第一方向，每间隔m个像素点，做频率坐标轴的频率分隔线，沿与所述第一方向垂直的第二方向，每间隔n个像素点，做波数坐标轴的波数分隔线，每相邻两个所述频率分隔线和相邻两个所述波数分隔线围合出一个所述网格。

[0075] 频率‑波数图中间隔分布有很多像素点，在划分网格时，频率方向和波数方向可以选用相同数量像素点作为网格边长，也可以选用不同数量像素点作为网格边长，例如，频率和波数方向可以分别将每5个像素点作为一个网格的边长。第一方向即沿频率坐标轴的方向，第二方向即沿波数坐标轴的方向，频率分隔线平行于波数坐标轴，波数分隔线平行于频率坐标轴，从而使得每一网格均呈现为矩形。

[0076] 可选地，每一所述网格内的功率谱，选用所述网格内所有所述像素点上的功率谱的平均值。

[0077] 在一个网格内，具有多个像素点，每一像素点上均有对应的功率谱的值，将每一网格内的所有像素点上的功率谱进行平均，能够降低奇异值，从而利用该功率谱平均值求取穿过该网格的斜率线的积分值，能够提高积分值结果的准确性。

[0078] 可选地，在所述将每一所述数据块对应的所述静态速度值进行拼接，得到所述全井段DAS数据对应的静态速度谱的步骤之前，还包括如下步骤：

[0079] S401：将每一所述数据块对应的所述静态速度值进行归一化处理。

[0080] 数据归一化通过将数据映射到相同的尺度范围内，不同特征之间的关系和变化将更容易观察和解释，可以帮助更好的理解数据的分布和趋势。将所有数据库对应的速度值
进行归一化处理后，得到归一化速度值，将得到的所有归一化速度值进行拼接后得到的全
井段范围内的静态速度谱，更能准确反映声波在油井的不同深度段之间的速度值的相对大
小。

[0081] 可选地，所述针对每一所述频率‑波数图，计算设定的每条所述斜率线的积分值，包括以下步骤：

[0082] S501：确定所述斜率线进入任一所述网格的进入点的坐标(x0，y0)，以及穿出该所述网格的穿出点的坐标(x1，y1)。

[0083] 在频率‑波数图中，横坐标为波数，纵坐标为频率；

[0084] 若所述斜率线从所述网格的侧边界穿出，则所述斜率线穿出所述网格侧边界的穿出点的坐标(xside，yside)为：

[0085] xside＝iΔx；

[0086] yside＝y0+m(xside‑x0)；

[0087] Δxside＝xside‑x0；

[0088] 其中，i为所述斜率线穿过的所述网格在波数坐标轴方向上的索引，取值为从1开始的正整数；

[0089] Δx为所述网格在波数坐标轴方向上的横向宽度；

[0090] m为所述斜率线的斜率值；

[0091] Δxside为所述斜率线穿过所述网格侧边界时在波数坐标轴上的增量；

[0092] 若所述斜率线从所述网格的上边界穿出，则所述斜率线穿出所述网格上边界的穿出点的坐标(xtop，ytop)为：

[0093] ytop＝jΔy；

[0094]

[0095] Δxtop＝xtop‑x0；

[0096] 其中，j为所述斜率线穿过的所述网格在频率坐标轴方向上的索引，取值为从1开始的正整数；

[0097] Δy为所述网格在频率坐标轴方向上的纵向宽度；

[0098] Δxtop为所述斜率线穿过所述网格上边界时在波数坐标轴上的增量；

[0099] 比较Δxside和Δxtop，

[0100] 若Δxside≤Δxtop，则所述斜率线穿出该所述网格的侧边界，所述斜率线穿出该所述网格的穿出点的坐标为：

[0101] x1＝iΔx；

[0102] y1＝y0+m(x1‑x0)；

[0103] 同时，所述斜率线穿过的所述网格在波数坐标轴方向上的索引i加1；

[0104] 若Δxside＞Δxtop，则所述斜率线穿出该所述网格的上边界，所述斜率线穿出该所述网格的穿出点的坐标为：

[0105] y1＝jΔy；

[0106]

[0107] 同时，所述斜率线穿过的所述网格在频率坐标轴方向上的索引j加1。

[0108] 图4为本发明一个实施例所述的频率‑波数图网格划分的示意图，图中括号内的数字表示网格的索引(i，j)。如图4所示，以图中的斜率线a为例进行说明，其斜率值为ma，所有的斜率线均经过频率‑波数图的原点O，因此，所有斜率线穿过的第一个网格均对应波数坐标轴方向上第一个网格，同时也是频率坐标轴方向上第一个网格，该网格索引(i，j)＝(1,
1)，斜率线进入索引为(1,1)的网格的进入点的坐标(x0，y0)即原点O坐标(0，0)。

[0109] 若斜率线从该网格的侧边界穿出，则穿出点的坐标(xside，yside)为：

[0110] xside＝iΔx＝Δx，此时i＝1；

[0111] yside＝y0+ma(xside‑x0)＝maΔx；

[0112] Δxside＝xside‑x0＝Δx；

[0113] 若斜率线从网格的上边界穿出，则穿出点的坐标(xtop，ytop)为：

[0114] ytop＝jΔy＝Δy，此时j＝1；

[0115]

[0116] 对于确定的斜率值ma， 对应图4中的长度值Δx1，确定地，Δx1>Δx，因此，Δxside＜Δxtop，斜率线a穿出索引为(1,1)的网格的侧边界，穿出点的坐标(x1，y1)＝(Δx，maΔx)。同时，因斜率线a穿出网格侧边界，因此，斜率线a要穿过的下一个网格对应的位置，波数坐标轴方向上索引i必然加1，而j不变，也即，斜率线a穿过的第二个网格的索引为(i，j)＝(2，1)。

[0117] 斜率线a穿出第一个网格的穿出点的坐标，即其进入第二个网格的进入点的坐标，因此，对于斜率线a穿过的第二个网格(i，j)＝(2，1)，其进入点坐标(x0′，y0′)＝(x1，y1)＝(Δx，maΔx)。

[0118] 若斜率线从该网格的侧边界穿出，则穿出点的坐标(xside，yside)为：

[0119] xside＝iΔx＝2Δx，此时i＝2；

[0120] yside＝y0+ma(xside‑x0)＝2maΔx；

[0121] Δxside＝xside‑x0＝Δx；

[0122] 若斜率线从该网格的上边界穿出，则穿出点的坐标(xtop，ytop)为：

[0123] ytop＝jΔy＝Δy，此时j＝1；

[0124]

[0125] 由图4可知， 因此，Δxside＞Δxtop，斜率线a穿出索引为(2，1)的网格的上边界，穿出点的坐标 同时，因
斜率线a穿出网格上边界，因此，斜率线a要穿过的下一个网格对应的位置，频率坐标轴方向上索引j必然加1，而i不变，也即，斜率线a穿过的第三个网格的索引为(i，j)＝(2，2)。

[0126] 依据以上步骤，依次类推，得到斜率线a穿过的每一网格的进入点坐标和穿出点坐标。

[0127] 另外，也可以采用其他方式判断斜率线从网格的侧边界还是上边界穿出。例如，网格划分确定后，网格的长度和宽度确定，则每一网格的四个顶点的坐标即可确定，在每一网格中，计算斜率线进入该网格的进入点与该网格右上角的顶点之间的连线的斜率值n，若连线的斜率值n大于斜率线的斜率值m，则该斜率线从该网格的侧边界穿出；若连线的斜率值n小于或等于斜率线的斜率值m，则该斜率线从该网格的上边界穿出。

[0128] S 502：计算所述斜率线在穿过的每一所述网格内的长度l，

[0129]

[0130] 在前述步骤中，已经获得斜率线穿过的每一网格的进入点坐标和穿出点坐标，根据勾股定理，即可计算得到斜率线穿过的每一网格内的长度l。

[0131] S503：计算所述斜率线的积分值I，

[0132]

[0133] 在每一数据块对应的频率‑波数图的网格范围内，对于每条斜率线，将斜率线穿过的每一网格内的长度l相加，即得到该斜率线的总长度L，将斜率线穿过的每一网格内的长度与该网格对应的功率谱相乘得到乘积值，将斜率线穿过的所有网格对应的乘积值相加，
再与该斜率线的总长度L相除，即得到该斜率线的积分值。

[0134] 按照上述过程，可以计算获取同一数据块中的每一斜率线的积分值，然后将所有斜率线的积分值进行比较，最大的积分值对应的斜率线，其斜率值即代表该数据块的声波
的静态速度值。

[0135] 同样的过程应用于不同井段深度的所有数据块，即可得到所有井段深度的声波的静态速度值，将所有井段深度的声波的静态速度值进行拼接，即得到全井段DAS数据对应的静态速度谱。

[0136] 可选地，根据所述静态速度谱v和井筒内的流体速度c，获取井筒内的上行波速度+ ‑
谱v与下行波速度谱v的步骤中，根据以下公式进行计算：

[0137] v+＝v+c；

[0138] v‑＝v‑c；

[0139] 从DAS数据获取的静态速度谱为声波在井筒中流体速度为零时的速度，而实际情况是，井筒中的流体向上朝向地面方向流动，向上传播的声波与流体流动方向相同，向下传+
播的声波与流体流动方向相反，因此，根据上述公式可以计算出上行波速度谱v 与下行波
‑
速度谱v。图5为本发明一个实施例所述的上行波速度谱示意图；图6为本发明一个实施例
所述的下行波速度谱示意图，图中空心圆圈代表该深度范围内积分最大值的斜率线对应的
斜率值计算出的上行/下行波速度值，实心圆点表示置信区间，声速范围表示根据各斜率线计算出的上行/下行波速度范围。

[0140] 进一步地，根据上行波速度谱v+和下行波速度谱v‑，还可用于双相流比例的计算，这是算法结果的直接应用，可以帮助更准确地监测和分析管道或油井中的流体组成。

[0141] 本发明还提供一种计算机可读存储介质，所述计算机可读存储介质中存储有可供计算机读取的指令信息，计算机读取所述指令信息后可执行上述任一实施例所述的从DAS
数据提取速度谱的方法。

[0142] 图7为本发明一个实施例所述的执行从DAS数据提取速度谱的方法的电子设备的硬件连接关系示意图。如图7所示，本发明还提供一种电子设备，包括至少一个处理器和至少一个存储器，至少一个所述存储器中存储有指令信息，至少一个所述处理器读取所述指
令信息后可执行上述任一实施例所述的从DAS数据提取速度谱的方法。

[0143] 如图7所示，电子设备还可以包括：输入装置23和输出装置24。处理器21、存储器22、输入装置23和输出装置24可以通信连接。存储器22作为一种非易失性计算机可读存储
介质，可用于存储非易失性软件程序、非易失性计算机可执行程序以及模块。处理器21通过运行存储在存储器22中的非易失性软件程序、指令以及模块，从而执行各种功能应用以及
数据处理，即实现上述任一实施例提供的从DAS数据提取速度谱的方法。存储器22可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储根据从DAS数据提取速度谱的方法的使用所创建的数据等。此外，存储器22可以包括高速随机存取存储器，还可以包括非易失性存储器，例如至少一个磁盘
存储器件、闪存器件、或其他非易失性固态存储器件。在一些实施例中，存储器22可选包括相对于处理器21远程设置的存储器，这些远程存储器可以通过网络连接至执行从DAS数据
提取速度谱的方法的装置。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。输入装置23可接收输入的用户点击，以及产生与从DAS数据提取速度谱的方法的用户设置以及功能控制有关的信号输入。输出装置24可包括显示屏等显示设备。
在所述一个或者多个模块存储在所述存储器22中，当被所述一个或者多个处理器21运行
时，执行上述任意方法实施例中的从DAS数据提取速度谱的方法。

[0144] 本发明还提供一种计算机程序产品，包括计算机程序/指令，该计算机程序/指令被处理器执行时实现上述任一实施例所述的从DAS数据提取速度谱的方法。

[0145] 最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；
而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。