Radhakrishnan 等人的信号分割算法移植到 C 和 C#

• 下载二进制文件 - 43.6 KB
• 下载源代码 - 140.2 KB

引言

Radhakrishnan 等人提出了一种用于分割石油行业测井信号的算法（您可以在 此处 找到他们的文章）。该算法可用于将源信号划分为信号可视为恒定但带有噪声的区域。区域的边界和区域内的信号值是先验未知的。原始算法是用 FORTRAN-77 编写的。它已发表，但电子形式的源代码不可用。他们的算法已被移植到 C，此外还提供了一个 C# 包装器。这些移植是本文的重点。

需要注意的是，原始算法包含一个错误，该错误已在当前实现中得到纠正。该错误允许读取数组末尾之外的内存。纠正已在下面和源代码中注明。这是当前算法与原始算法唯一不同的地方。

除了算法和包装器之外，还提供了单元测试。本文中的示例来自单元测试。

本文将不深入探讨数学原理，因为算法的作者在这方面做得相当出色。但是，如果您没有阅读过原始论文，本文将提供足够的背景知识来理解后面的动机和示例。

背景

请看下面的图，假设我们想在给定由信号和一些噪声（不需要的信号随机干扰）组成的测量值（记录的值）的情况下，找到未知的信号（所需信息）。

这是“逆问题”的一种类型。在逆问题中，我们计算生成给定输出所需的输入。逆问题出现在许多科学和数学领域。Radhakrishnan 等人的算法是为一种特定类型的逆问题设计的；它将一个允许变化其值的测量信号分割成信号可视为恒定的区域，同时允许存在噪声。作者提出此算法的动机是为了分割石油和天然气行业钻取的井的测井数据。作者将他们的技术总结为：

该方法基于单变量状态变量模型，在该模型中，状态（信号）在且仅在段边界处以随机量变化。观测到的测井值是状态和随机噪声的总和。噪声被建模为高斯噪声，但不一定是平稳的。测井分割的逆问题是通过迭代应用卡尔曼滤波和估计理论的单一最可能替换（SMLR）技术来解决的。

库的设计

代码的设置是为了创建算法的三个可用版本：C 语言静态库、C 语言 DLL 和 C# 语言 DLL。这是为了提供一系列的最终用途选项。然而，为了尽量减少开发工作量并简化代码的维护，只实现了一个版本的算法；其他版本是通过包装该版本来创建的。该算法是用 C 语言编写的，因为 C 语言既高效又被广泛使用。这也使得代码可以在 C 和 C++ 中使用，或者包装后在 C# 甚至其他环境（如 Matlab（未提供））中使用。虽然包装会引入一些额外的计算开销，但这被认为是可以接受的权衡，以换取更少的工作量。原因如下：第一，包装只引入了少量的开销要求。第二，在 Matlab 或 C# 中工作时，绝对的最高计算效率通常不是必需的。这些环境通常选择在其易于开发性 outweighs 引入的额外开销时。

Visual Studio 解决方案包含五个项目。它们列在下面，并附有简要摘要。请注意，算法本身仅在 C Static Library 中。

• C Static Library - 包含论文中的信号分割算法，从原始 FORTRAN 源代码移植到 C
• C DLL - 创建代码的 DLL 版本。链接并调用 C Static Library。不包含除 DLL 导出函数以外的代码
• C Unit Tests - C Static Library 的单元测试
• C Sharp DLL Wrapper - 通过导入 C DLL 来创建 C# DLL。还提供接口类
• C Sharp Unit Tests - C Sharp DLL Wrapper 的单元测试

单元测试

没有创建示例项目来演示库的使用，但提供了包含示例的单元测试。单元测试利用 Visual Studio 的 Test Explorer/Unit Testing。此处不提供有关使用 Test Explorer/Unit Testing 的解释，因为其他地方有许多教程。一旦解决方案被打开和编译，应该会发现以下单元测试。

提供单元测试是为了在需要修改算法时可以进行测试。通常，用户除了确认下载和编译的库正常工作外，可能不会太多地需要执行单元测试。因此，用户将主要关注 LogSeparateTest 和 LogSeparateTestCS，它们演示了算法的使用。下面提供了有关它们使用方法的详细信息。这些函数重现了论文中的示例，以证明实现是正确的。在下图中，论文中提供的输出已与 LogSeparateTestCS 测试的输出并列显示。

Using the Code

引言

C 和 C# 版本的函数签名有所不同。C 语言版本复制了论文中的 FORTRAN 签名。C# 版本是“现代化”的。这主要包括将结果返回到单个类中，而不是通过输入参数传递。还使用了重载以允许默认值，并允许以列表或数组的形式输入。

为 C 版本提供了一个函数，其签名如下：

// Input
void SegmentSignal(double LOG[], int NSAMPS, double F, 
                   int ORDER, int ORDER1, int RMODE, int NITER,
// Output
double Q[], int &NQ, double FLTLOG[], double SEGLOG[], double R[], double &C,
double &D, int &NACT, int &IER);

总共有六个 C# 中的算法接口函数，三个用于将信号输入为 double 数组，三个用于将信号输入为 List<double>。下面提供了一个示例，以及返回类型 SegmentaitonResults 的签名。有关完整的调用选项列表，请参阅 C Sharp DLL Wrapper 项目中的 SegmentSignal.cs。

public static SegmentationResults Segment(double[] signal, double threshold,
	int jumpSequenceWindowSize, int noiseVarianceWindowSize,
	NoiseVarianceEstimateMethod noiseVarianceEstimateMethod, int maxSMLRIterations)

public class SegmentationResults
{
	...
	public double[]	BinaryEventSequence {get;}
	public int	NumberOfBinaryEvents {get;}
	public double[]	FilteredSignal {get;}
	public double[]	SegmentedLog {get;}
	public double[]	NoiseVariance {get;}
	public double	JumpSequenceVariance {get;}
	public double	SegmentDensity {get;}
	public int	Iterations {get;}
	public int	Error {get;}
	...
}

输入和输出

该算法有几个输入和输出参数。下面表格中显示了两种语言的参数以及简要说明。

C	C#	描述
LOG	signal	待分割的输入信号
NSAMPS	未使用	输入数组的长度（LOG 的长度）
F	threshold	分割阈值。用于调整算法分割信号的精细程度。值越小，段越多；值越大，段越少。
ORDER	jumpSequenceWindowSize	用于平滑输入测井数据以获得跳跃序列方差初始估计的移动平均窗口的长度。
ORDER1	noiseVarianceWindowSize	用于平滑噪声方差的移动平均窗口的长度。建议将其设置为 ORDER/jumpSequenceWindowSize 的大约一半。
RMODE	noiseVarianceEstimateMethod	噪声方差估计选项。您可以指定点估计（C 中为 0，C# 中为 NoiseVarianceEstimateMethod.Point）或平滑估计（C 中为 1，C# 中为 NoiseVarianceEstimateMethod.Smoothed）。点估计速度更快，但平滑估计可能提供更好的结果。
NITER	maxSMLRIterations	单一最可能替换（SMLR）迭代次数的上限。

下面表格中列出了输出参数。请注意，C# 参数在 SegmentationResults 的实例中返回。它们作为属性访问（例如 SegmentationResults.BinaryEventSequence）。

C	C#	描述
Q	BinaryEventSequence	二进制事件序列。一个与输入信号/测井数据大小相同的数组，在检测到二进制事件的位置为“1”，其他位置为“0”。二进制事件是算法确定信号值发生变化的地点。
NQ	NumberOfBinaryEvents	找到的二进制事件数（Q/BinaryEventSequence 中的“1”的数量）。
FLTLOG	FilteredSignal	信号的过滤估计。
SEGLOG	SegmentedLog	每个段的过滤测井数据（FLTLOG/FilteredSignal）的平均值。
R	NoiseVariance	估计的噪声方差。
C	JumpSequenceVariance	跳跃序列的估计方差。
D	SegmentDensity	估计的段密度。
NACT	迭代	实际 SingleMostLikelihoodReplacement 迭代计数。
IER	Error（错误）	错误标志。零：无错误。小于零：阈值后估计的事件密度无效。减少/增加阈值并重新运行。大于零：在 SMLR 计算似然比期间，对数参数在返回的样本号处变为零。可能存在更多此类样本，导致此问题。编辑/重新缩放数据值并重新运行。

从输出值来看，用户最可能感兴趣的是二进制事件序列（Q/BinaryEventSequence）和分割后的测井数据（SEGLOG/SegmentedLog）。它们包含了区域之间的边界位置和这些区域的“平均”值。

使用 C 语言版本

在 C Unit Tests 库的 LogSeparateUnitTest.cpp 文件中的 LogSeparateTest 单元测试下可以找到使用 C 语言版本的示例。该函数中与调用算法相关的部分如下所示。_data[Log] 变量包含输入数据，该数据是通过调用 ReadData 函数从文件中读取的。

// Input values specified in paper. These do not have to be variables. They can be hard coded.
double f	= 2.50;
int order	= 10;
int order1	= 6;
int rmode	= 0;
int niter	= 300;

// Output variables. These must be variables so they can be used to return data.
double Q[_numberOfDataPoints];
int numberOfBinaryEvents;
double FLTLOG[_numberOfDataPoints];
double SEGLOG[_numberOfDataPoints];
double R[_numberOfDataPoints];
double C;
double D;
int NACT;
int error;

// Call to separate the signal.
SegmentSignal(_data[Log], 100, f, order, order1, rmode, niter,
	            Q, numberOfBinaryEvents, FLTLOG, SEGLOG, R, C, D, NACT, error);

使用 C# 语言版本

在 C Sharp Unit Tests 库的 LogSeparateUnitTest.cs 文件中的 LogSeparateTestCS 单元测试下可以找到使用 C# 语言版本的示例。该函数中与调用算法相关的部分如下所示。_data[(int)InputType.Log] 变量包含输入数据，该数据是通过调用 ReadData 函数从文件中读取的。

请注意，在使用 C# 版本时，必须将 C 语言 DLL 放在同一个目录中。也就是说，它需要 SegmentSignalCS.dll 和 SegmentSignal.dll。此要求源于 C# 版本包装并调用 C 版本。

// Input values specified in paper.
private const double            _f          = 2.50;
private const int               _order      = 10;
private const int               _order1     = 6;
private const int               _niter      = 300;
NoiseVarianceEstimateMethod     _rmode      = NoiseVarianceEstimateMethod.Point;

...

// Call to separate the signal.
SegmentationResults results = SegmentSignal.Segment(_data[(int)InputType.Log], _f, _order,
                                                    _order1, _rmode, _niter);

算法的修改

在使用和测试算法的过程中，发现对于某些输入，可能会访问数组末尾之外的内存。在原始论文中，有问题的函数称为 GETSPK。当前实现中的新函数 GetSpikeyZone 如下所示：

// Original function name: GETSPK
void GetSpikeyZone(double Q[], int NSAMPS, int SEGSTR, int &SEGEND, bool &endofarray)
{
	// Function:
	// Identifies a zone of spikes.  A spiky zone is defined by 0,1,1,...,1,0.
	//
	// Input parameters:
	// Q: Binary event sequence in which the zone is to be identified.
	//
	// NSAMPS: Length of the array Q.
	//
	// SEGSTR: Sample number where the spiky zone starts.
	//
	// Output parameters:
	// SEGEND: Sample number where the zone ends.
	//
	// endofarray: End of data indicator.
	//		true: End of array reached.
	//		false: Not at end of the array.

	SEGEND = SEGSTR;
	if (Q[SEGEND] < 0.5)
	{
		// Modified from original paper. 
        // The "if" statement below was not part of the original
		// algorithm.
		//
		// It seems there can be a case where Q[i] < 0.5 for SEGEND = NSAMPS 
        // that causes us to miss
		// flagging the end of the array. 
        // Do a check here to ensure we don't overrun the end of
		// the array and over write memory we should not over write.
		if (SEGEND >= NSAMPS)
		{
			endofarray	= true;
			SEGEND		= NSAMPS - 1;
		}

		return;
	}

	do
	{
		SEGEND++;
	}
	while ((Q[SEGEND] > 0.5) && (SEGEND < NSAMPS));

	if (SEGEND >= NSAMPS)
	{
		endofarray = true;
	}

	SEGEND--;
}

此处添加了以下内容以防止返回时数组索引大于数组大小：

if (SEGEND >= NSAMPS)
{
	endofarray	= true;
	SEGEND		= NSAMPS - 1;
}

结论

本文仅提供了算法从 FORTRAN 到 C 和 C# 的移植，并提供了二进制文件和源代码供立即使用。对原始算法的唯一贡献是纠正了一个在某些条件下允许读取数组末尾之外的错误。

虽然算法最初的动机是用于分割测井数据，但该算法在其他地方也有应用。本文作者已将其用于分割测井数据和各种实验室实验产生的数据。它可以作为自动化信号分析过程的一部分；可以检测和处理分离的区域，而无需用户识别感兴趣的区域。

致谢

感谢 APS Technology 允许发表本文。

历史

• 2015 年 6 月 - 首次发布
• 2019 年 10 月 - 语法和拼写校正