第一步：标准化与等距重采样

原始的 GR 测井曲线不仅有量纲（API），而且因为仪器拉升，深度采样是不均匀的。算法的第一步，必须把数据“洗”成纯粹的、等间距的数学序列。

1. 消除量纲：Z-Score 标准化

对应代码：_safe_standardize(data)

数学公式：

符号含义：

• $x(z)$：在深度$z$ 处的原始 GR 值。

• $\mu$：整段曲线 GR 的算术平均值。

• $\sigma$：整段曲线 GR 的标准差（反映波动幅度）。

• $x'(z)$：无量纲、均值为 0、方差为 1 的相对波动信号。

  底层逻辑：把绝对的 GR 读数变成“相对变化”。不管这个地层是整体偏泥（GR高）还是整体偏砂（GR低），我们只关心它围绕局部均值是怎么波动的。

2. 强行等距：Pchip 重采样

对应代码：resample_to_equal_depth(depth, values)

物理现实：真实的深度序列是$z_1, z_2, z_3 \dots$，它们之间的差值$\Delta z_i = z_{i} - z_{i-1}$ 是变化的。

数学动作：代码强行生成一个绝对等距的理想深度轴$z_{eq}$，其固定步长为$\Delta z_{eq} = \text{median}(\Delta z_i)$（取原始步长的中位数）。

然后用 Pchip（分段三次埃尔米特插值）求出理想深度的值：

底层逻辑：傅里叶变换、EMD、希尔伯特变换，全部要求$X$ 轴必须是等距的。如果物理深度不等距，算法就会把空间上的“拉伸”误认为是地质上的“低频”，导致最终提取的旋回周期完全错乱。

第二步：EMD (经验模态分解) —— 频率的自适应剥离

第一层大逻辑：我们的终极目标是什么？

看这个公式：

通俗翻译：

• $x_{eq}(z)$：就是那团原始的“乱麻”（经过等距处理的 GR 曲线）。

• $\text{IMF}_k(z)$：就是被我们抽出来的一根根“纯粹的线”。

  • $\text{IMF}_1$ 是最细、抖动最快的丝线（高频噪波）。

  • $\text{IMF}_2, \text{IMF}_3$ 是中等粗细、抖动平缓的毛线（米级、十米级旋回）。

• $r(z)$：残差（Residue）。当把所有能抖动的线都抽走后，最后剩下的那根死气沉沉、几乎不怎么弯曲的光秃秃的“大铁柱子”（盆地大背景趋势）。

等式左边 = 等式右边。这意味着：EMD 不是凭空捏造数据，它只是把原始曲线拆开了。把剥出来的这些线重新叠在一起，还是原来那根一模一样的 GR 曲线。

第二层核心：自动化剥线机（Sifting 筛分）是怎么运转的？

怎么把最细、抖动最快的那根线（$\text{IMF}_1$）先抽出来？]

代码里的 for _sift in range(max_sift_iter): 就是这台剥线机在反复运转。我们来看它运转一次的三个物理动作：

动作 1：蒙床单与拉吊床（求包络线）

面对一根上下乱窜的曲线$h(z)$，机器怎么知道它的轮廓？

• 找波峰：把曲线上所有往上凸的最高点（极大值）挑出来，用一根有弹性的软尺（三次样条插值），顺滑地把这些最高点连起来。这就好比在曲线上方蒙了一层床单。这叫上包络线$E_{up}(z)$。

• 找波谷：把所有往下凹的最低点（极小值）挑出来，用软尺顺滑地连起来。这就好比在曲线下方挂了一张吊床。这叫下包络线$E_{low}(z)$。

数学解剖：三次样条是如何做到“连续可导”的？

• 现在，我们把图板上的钉子换成数学坐标系里的波峰点 $(z_1, y_1), (z_2, y_2), \dots, (z_n, y_n)$。

  这里的$z$ 是深度， $y$是波峰的 GR 值。

  注意：样条插值并不是用一个巨大的高维多项式去硬套所有的点！（那种叫拉格朗日插值，在点多的时候会产生恐怖的龙格现象，边缘会疯狂震荡甩飞）。

  样条插值的核心思想是**“分段治理”**。

  第一步：分段构造三次多项式

  在每一对相邻的波峰$(z_i, y_i)$ 和$(z_{i+1}, y_{i+1})$ 之间，我们单独给它分配一个小方程$S_i(z)$。

  为了模拟竹片的弹性，数学家发现最高次项只需用到$3$ 次就足够了：

  $$S_i(z) = a_i + b_i(z - z_i) + c_i(z - z_i)^2 + d_i(z - z_i)^3$$

  （注意：这里只有$z$ 的最高次幂是$3$ 次，而不是把所有$N$ 个点放到一起搞一个$N-1$ 次的超级方程）。

  如果我们有$N$ 个波峰，我们就有$N-1$ 段这样的小方程。每一段都有四个未知的参数：$a_i, b_i, c_i, d_i$。

  第二步：缝合接口，强行要求连续与可导

  这是最关键的一步。我们怎么保证第一段曲线$S_1$ 和第二段曲线$S_2$ 在它们交接的那个波峰点$(z_2, y_2)$ 上，接得像竹片一样天衣无缝？

  数学上，我们用三把“锁”把它们死死锁住：

  1. 位置锁（连续性）：必须接头

     第一段曲线走到$z_2$ 时的值，必须等于第二段曲线在$z_2$ 开始时的值，并且它们都必须等于那个真实的波峰高度$y_2$。

     $$S_1(z_2) = S_2(z_2) = y_2$$

     物理意义：竹片在这个钉子处没有断开。

  2. 一阶导数锁（斜率连续）：必须顺滑

     第一段曲线走到$z_2$ 时的斜率，必须绝对等于第二段曲线在$z_2$ 开始时的斜率。

     $$S'_1(z_2) = S'_2(z_2)$$

     物理意义：接口处没有折角，不可能出现“尖点”。竹片是顺着原本的方向滑过去的。

  3. 二阶导数锁（曲率连续）：必须受力均匀

     第一段曲线走到$z_2$ 时的曲率（弯曲的急缓程度），必须等于第二段曲线开始时的曲率。

     $$S''_1(z_2) = S''_2(z_2)$$

     物理意义：竹片在钉子处的受力是均匀过渡的，没有瞬间的“生硬扭曲”。这也正是它被称为“最平滑”插值的原因。

  当把这三个条件列出来，结合边界条件（比如假设最两端的波峰处没有外力弯折，二阶导数为 0，这叫自然边界），数学家就能通过解一个三对角线性方程组，瞬间求出所有分段方程里那些未知的$a, b, c, d$ 参数。

  就这样，一条完美的、分段的三次多项式包络线，被严丝合缝地拼接了出来。

  ---

  3. Pchip：保形插值，给“竹片”加一层刹车

  在提供的代码里，虽然常规 EMD 常用普通的 Cubic Spline，但用到了一种更硬核的防御机制：Pchip（分段三次埃尔米特插值多项式）。

  PchipInterpolator(depth, values)

  为什么要换用 Pchip？普通竹片（Spline）有什么致命缺陷？

  普通的 Cubic Spline 为了追求极致的二阶导数连续（全局最平滑），它在面对步长极端突变的点时，会**“用力过猛”**。

  设想一个场景：的三个波峰高度分别是 100, 105, 105。

  前两个波峰$y$ 值在长，后两个波峰$y$ 值平了。

  如果拿竹片硬压过去，为了保持曲率的绝对平稳过渡，竹片在穿过第二个 105 的时候，往往会向上高高地弹飞一下（比如弹到 120），然后再落回到第三个 105 上。

  这在数学上叫过冲现象（Overshoot）。在 EMD 里，一旦包络线发生了过冲，算出来的均值线就会严重扭曲，导致剥离失败。

  Pchip 是怎么解决这个问题的？

  Pchip 牺牲了一点点极致的平滑（它放弃了强制要求二阶导数连续），而换取了一个更为关键的物理特性：保形性（Monotonicity-preserving）。

  Pchip 就像是一根有记忆、带刹车系统的特殊竹片。它在拼接接口处的斜率（一阶导数）时，不再只考虑怎样最顺滑，而是增加了一个强制规定：

  如果在原始数据中，这几个波峰并没有继续往上长，那么我画的曲线也绝对不允许向上弹飞哪怕一毫米！

  它严格保持了原始极值点之间的单调性。波峰平了，它就乖乖地平着过去；波峰往下降了，它就老老实实地往下走。绝不制造原数据中不存在的虚假极值。

动作 2：找重心（计算均值面）

床单和吊床之间是有厚度的。我们把床单的高度和吊床的高度相加，除以 2，找到它们正中间的那条线。

• $m(z)$：这根均值线，其实代表了这团线当前的“粗大骨架”（较慢的低频趋势）。而那些细小的快速抖动，正骑在这个骨架上。

动作 3：把骨架抽走（剥离局部趋势）

• $h_{old}(z)$：没抽骨架前的曲线。

• $m(z)$：刚刚算出来的骨架（均值线）。

• $h_{new}(z)$：减法在这里的物理意义是“拍扁”。我们把那根弯弯曲曲的骨架$m(z)$ 强行拉直，按死在$0$ 刻度线上。

  此时，原本骑在骨架上的那些快速抖动，就被迫老老实实、对称地围绕着$0$ 轴上下翻飞了。

第三层判断：怎么知道剥干净了没有？（停止准则）

剥线机做完上面三个动作后，出来的$h_{new}(z)$ 就是一根完美的纯粹细线（IMF）了吗？

绝对不是。

因为一次“蒙床单、找中间、拍扁”的操作，往往不能把骨架抽干净。所以代码要用 for 循环，对着这个新出来的$h_{new}(z)$，再蒙床单、再找均值、再减去均值。就像拧海绵里的水一样，拧一次不干，得反复拧。

那机器什么时候停手，承认“好，这是一根完美的纯粹的线（IMF）”了？必须同时通过两个极其严苛的质检：

质检 1：物理形态必须“纯洁” (_is_imf)

一根纯洁的线（单一频率），它的物理形态必须极其规律：

穿过一次 0 轴，只能对应一个波峰（或波谷）。

• 如果一根曲线穿过 0 轴往上走，在掉下来穿过 0 轴之前，居然长出了两个甚至三个小鼓包，说明什么？说明这根线里还藏着更细的线！这不纯洁，打回去继续拧（继续 Sifting）。

• 数学定义：极值点（波峰+波谷）的总数，和穿过 0 轴的交点总数，最多只能相差 1。

质检 2：实在拧不出水了（数学收敛条件 SD）

怎么衡量海绵里没水了？看这次拧出来的水，和上次拧出来的水，差了多少。

把这个公式掰开：

• 分子$(h_{old}(z) - h_{new}(z))^2$：上一次的曲线减去这一次的曲线。这就是**“这一次操作拧出来的水分（偏差）”**。

• 分母$h_{old}^2(z)$：上一次曲线原本的体量。

• 相除：这其实是在算一个**“相对变化率”**。

• 小于 0.20：如果这次大动干戈剥了一遍，结果曲线跟上一次相比，变化连 20% 都不到（代码设定的阈值）。说明骨架已经被抽干了，剩下的完全是纯粹的快速抖动。

第四层闭环：剥完一根，然后呢？

一旦$h_{new}(z)$ 通过了上面两项质检，系统就大喊一声：“停！$\text{IMF}_1$诞生了！”

接下来怎么做？

把这根最细的线，从最开始的那团乱麻里彻底剪掉：

现在，原始曲线没有了高频噪音，变得平滑了一点。

剥线机拿着这个剩下的$\text{Residue}$，从头开始，继续蒙床单、找均值、拍扁……

去寻找第二细的线$\text{IMF}_2$。

如此循环往复，直到最后的$\text{Residue}$ 连 3 个波峰都凑不齐（无法蒙床单了），整个 EMD 宣告结束。

这就是 EMD。没有高深的玄学，只有极其笨拙但极其有效的物理分解。

这一步对应核心函数 emd()。测井曲线里混杂了仪器高频噪声、薄夹层、米级旋回、数十米级层序。EMD 的任务是把它们一层层剥开。

总目标公式：

符号含义：

• $x_{eq}(z)$：第一步处理好的等距标准化信号。

• $\text{IMF}_k(z)$：第$k$ 个本征模态函数（代表某种单一频率的波动，$k$越小频率越高）。

• $r(z)$：残差（单调不波动的全局大趋势）。

第三步：自适应带通滤波 —— 拼装目标旋回

第一层：我们在干什么？（数学与代码的对应）

核心代码：

target = np.sum(selected_imfs, axis=0)

数学公式：

我们把公式里的每一个零件掰碎了看：

• $z$：深度。代表我们是一米一米、一个数据点一个数据点在往下算。

• $k$：编号。代表第几根线（第几个 IMF）。

• $S$：下达的“订单”（即代码里的 selected_imfs_1based 参数）。假设传入的是 [3, 4]，那么集合$S$ 就是$\{3, 4\}$。

• $\sum$：累加求和。

• $T(z)$：目标曲线（Target）。也就是代码里的 target 变量。

底层物理动作（axis=0 是什么意思？）：

在计算机内存里，挑出来的$\text{IMF}_3$ 和$\text{IMF}_4$ 是两条长度完全一样（比如都是 9000 个深度点）的数据序列。

np.sum(..., axis=0) 并不是把它们首尾相接，而是把它们像两块长条积木一样上下对齐，然后在每一个相同的深度点上，把数值加起来。

• 如果在 2000 米处，$\text{IMF}_3$的值是$5$，$\text{IMF}_4$的值是$-2$。

• 那么组装后的$T(2000)$ 就是$5 + (-2) = 3$。

就这么简单。把选中的几根线，按深度点逐一揉成一根稍微粗一点的线。

第二层：为什么要这么干？（地质视角的“去粗取精”）

既然 EMD 费了九牛二虎之力把它们全拆开了，为什么我们又要挑几根把它们加回去？

因为原始的 GR 曲线里，包含了太多根本不关心的东西。我们来看看地质学家是怎么给这些线分类的：

1. 为什么扔掉$\text{IMF}_1, \text{IMF}_2$？（扔掉“高频”）

这两根线抖动得极其频繁（频率最高）。

• 物理本质：它们往往代表测井仪器在井下摩擦产生的机械噪音、井壁不规则导致的垮塌干扰，或者是地层里只有十几厘米厚的极薄粉砂岩夹层。

• 工程决定：我们要研究的是几十米级别的大旋回，这些高频的毛刺只会干扰我们的视线，必须当成“工业废料”扔掉。

2. 为什么扔掉$\text{IMF}_5, \text{IMF}_6 \dots$ 和残差$r(z)$？（扔掉“低频”）

这些线长得非常平缓，几百米才弯曲一次。

• 物理本质：它们代表了极其宏大的地质事件——比如整个盆地经历了长达几百万年的持续沉降，或者全球海平面发生了极其缓慢的上升。这种大背景趋势，会导致 GR 曲线整体慢慢变大或变小。

• 工程决定：这叫“低频基线”。如果还记得，捕捉这种大尺度斜坡，恰恰是上一篇博客里 INPEFA 算法最擅长干的事！而我们现在的 INCOSIN 算法，是为了找特定级别的周期性“节律”。大斜坡会掩盖小周期，所以必须把它们像“剔骨头”一样剔除掉。

3. 为什么留下$\text{IMF}_3, \text{IMF}_4$？（保留“中频”）

• 物理本质：当去掉了极高频的噪波，抽走了极低频的骨架，剩下的$\text{IMF}_3$ 和$\text{IMF}_4$ 加在一起，恰好对应了地质学上的三级或四级层序（数十米到百米级别的沉积旋回）。

• 结果：得到的这根$T(z)$ 曲线，就是干干净净的、只包含目标地层周期波动的纯净旋回曲线。

• 

优化方向：

如果在算法里把参数写死成 selected_imfs = [4, 5]（或者 [3, 4]），然后告诉别人“这就代表地质旋回”，这不仅是典型的**“拍脑门”工程学**，而且在严谨的数学地质和学术答辩上是绝对站不住脚的，一戳就破。

经验模态分解（EMD）最大的痛点和争议，就在于它的模态混叠（Mode Mixing）和物理意义的不确定性。

但是，通过对提供的这张高分辨率 EMD 解剖图进行客观的信号与地质尺度分析，我们可以推导出现阶段选择 IMF4 + IMF5 的潜在合理性。要让这个选择从“拍脑门”变成“有理有据”，必须从以下三个硬核的物理维度进行论证：

第一维度：视觉频率与地质厚度的“尺度映射”

抛开数学，我们先用肉眼直接读取图版上的空间频率（主导波长/地层厚度）：

的 Y 轴深度跨度大概是$2000m$ 到$4500m$（总长约$2500m$）。

1. $\text{IMF}_1 \dots \text{IMF}_3$（被抛弃的高频）

   • 形态观察：密集如杂草，过零点（Zero-crossings）极多。

   • 尺度估算：在这$2500m$ 内震荡了数百次，单一旋回厚度约为$1m \sim 5m$。

   • 地质判别：这种尺度的波动，绝大多数来源于测井仪器的纵向分辨率极限（机械噪音）、井壁微小不规则（扩径）、或是厘米到亚米级的粉砂/泥岩极薄互层。对于宏观储层对比，这些是无效高频噪波。

2. $\text{IMF}_4, \text{IMF}_5$（被选中的目标 T(z)）

   • 形态观察：呈现出非常清晰的“纺锤状”包络，波形舒缓且有规律。

   • 尺度估算：在这$2500m$ 内震荡了几十次。单一波长（一个完整旋回的厚度）大约在$20m \sim 80m$ 之间。

   • 地质判别：这在沉积学上极其致命！$20m \sim 80m$的沉积厚度，恰好完美对应了浅海/三角洲相的“四级层序（准层序组 Parasequence Sets）”或是米兰科维奇旋回（偏心率长周期）的地层表现。 这是石油工业中最核心的储层追踪级别。

3. $\text{IMF}_6 \dots \text{IMF}_8$ 和 Residue（被抛弃的低频）

   • 形态观察：波形极其宽缓。

   • 尺度估算：单一波长高达$200m \sim 500m+$。

   • 地质判别：这是三级层序（体系域级别）、甚至是二级构造沉降大背景。留着它们会产生严重的“低频漂移”，掩盖内部细节。

结论 1： 选$\text{IMF}_4 + \text{IMF}_5$，本质上是在做一个**“物理尺度带通滤波”**，截取的是$20m \sim 80m$ 地层厚度的岩性波动。这在特定的研究区（如珠江口盆地恩平/文昌组）是有地质依据的，但如果是其他沉积速率极快的盆地，这个尺度可能就跑到$\text{IMF}_5 + \text{IMF}_6$ 去了。

第二维度：方差贡献率（能量判定法）

如果要在论文或汇报中证明“选 4 和 5 是对的”，不能光靠肉眼看，必须拿出定量的数据支撑。

在信号处理中，我们用**能量（方差）**来判定一个 IMF 是真正的物理信号，还是随机噪音。

对于一个$\text{IMF}_k$，其能量（方差）定义为：

然后计算它占所有 IMF 总能量的百分比：

客观规律是：

• 纯粹的白噪声分解后，其 IMF 能量随着模态阶数$k$ 的增加，呈指数级递减（$\text{IMF}_1$能量最大，$\text{IMF}_8$极小）。

• 真实的、具有地质意义的周期信号，会打破这个指数递减规律，在某个特定的 IMF（比如$\text{IMF}_4$ 或$\text{IMF}_5$）上出现“能量异常突起”。

结论 2： 必须在代码里加一个模块，计算这口井各个 IMF 的方差占比。如果画出柱状图，发现$\text{IMF}_4$ 和$\text{IMF}_5$ 的能量出现了显著的波峰，这就成了选择它们的最硬核、最无可反驳的数学证据。

第三维度：与已知地质界面的绝对标定（测井与岩心的闭环）

再精妙的数学，如果脱离了岩石，就是空中楼阁。要彻底洗脱“拍脑门”的嫌疑，唯一的方法是实证。

在图版中，提取出的$\text{INCOSIN}$ 曲线（倒数第二道）在特定深度会出现极其干脆的$+1$（最泥）和$-1$（最砂）。

• 做法： 拿的岩心记录、录井描述、或者地震解释的已知层序界面（Sequence Boundaries, SB）和最大海侵面（Maximum Flooding Surfaces, MFS），直接怼到这个深度坐标上！

• 如果： 已知的最大泥岩段/测井高伽马段，刚好完美对应$\text{IMF}_4+\text{IMF}_5$ 组合算出来的$\text{INCOSIN} = 1$ 的波峰。

• 如果： 已知的砂岩冲刷面，刚好完美对应$\text{INCOSIN} = -1$ 的波谷。

结论 3： 这种与客观地质事实的“空间吻合”，是证明的 IMF 组合具有实际物理意义的终极闭环。

第四步：希尔伯特变换 (Hilbert) —— 一维到二维的数学升维

前置知识：欧拉公式（破解秘密的钥匙）

在弄懂算法之前，我们必须复习一个高中/大一的数学公式，欧拉公式：

这公式反过来写，余弦波（实数信号）可以表示为：

注意看这个等式！

在复数的世界里，一个普普通通的实数波$\cos(\theta)$，其实是由两个东西拼成的：

1. $\frac{1}{2} e^{j\theta}$：这是一个逆时针旋转的向量（我们叫它正频率）。

2. $\frac{1}{2} e^{-j\theta}$：这是一个顺时针旋转的向量（我们叫它负频率）。

核心天机：任何一个实数信号（包括的测井曲线），在频域里，天然包含完全对称的正频率和负频率。

解剖手术台：hilbert(target) 到底干了什么？

当在 Python 里执行 analytic_signal = hilbert(target) 时，CPU 其实在后台飞速执行了 3 个极其冷酷的动作。

动作 1：把数据送进频域 (FFT)

操作：对那一维的实数数组$T(z)$，执行一次快速傅里叶变换（Fast Fourier Transform, FFT）。

发生了什么：的深度域曲线变成了一张频谱图。在这个图上，因为$T(z)$ 是纯实数，所以它的频谱在正频率（$+f$）和负频率（$-f$）上，出现了完美对称的峰值。

就像上面那个欧拉公式写的，它被拆成了$\frac{1}{2} e^{j\theta}$ 和$\frac{1}{2} e^{-j\theta}$。

动作 2：极其野蛮的“外科切除” (造虚部的核心)

操作：在频谱图里，拿起一把刀，把所有小于 0 的负频率部分，直接砍掉（数值强行赋值为 0）。然后把大于 0 的正频率部分的振幅，直接乘以 2。

发生了什么惊天动地的数学质变？

我们回到刚才的欧拉公式：

原来是：$\cos(\theta) = \frac{1}{2} e^{j\theta} + \frac{1}{2} e^{-j\theta}$

现在，我们砍掉了负频率$e^{-j\theta}$，并把正频率$e^{j\theta}$ 乘以了 2。

得到的新信号直接变成了：

而根据欧拉公式，$e^{j\theta}$ 展开是什么？

看懂了吗！！！

我们仅仅是在频谱上执行了“砍掉一半，剩下一半翻倍”的野蛮操作。

我们在数学上，硬生生地把原本只有实部的$\cos(\theta)$，变成了一个自带虚部的$\cos(\theta) + j\sin(\theta)$！

而$\sin(\theta)$，不正是$\cos(\theta)$ 往后推迟了$90^\circ$ 的完美影子吗！

动作 3：退回现实坐标系 (IFFT)

操作：对被“切除并翻倍”后的新频谱$F_{new}(\omega)$，执行快速傅里叶逆变换（IFFT），把它变回到的测井深度域$z$。

此时，从 IFFT 输出的数组 Z(z)（也就是代码里的 analytic_signal），再也不是原来的 [1.2, 3.4, ...] 这种纯实数数组了。

它变成了一个复数数组：

[1.2 + 0.5j, 3.4 - 1.2j, ...]

• 它的实部 np.real(Z)，依然完美地等于输进去的原始目标曲线$T(z)$。

• 它的虚部 np.imag(Z)，就是那个在数学上严格向后相移了$90^\circ$ 的影子信号$\mathcal{H}[T(z)]$。

第一层：希尔伯特变换 (Hilbert) —— 制造一个完美的“影子”

要构建二维平面的指针，我们的原信号$T(z)$ 可以作为横坐标（$x$轴）。我们现在急需一个纵坐标（$y$轴）的信号。

这个$Y$ 轴信号不能随便给，它必须和原信号保持绝对的同步，但又要错开一定的节奏。

数学公式：

别被这个带无穷积分的公式吓到，我们把符号一个个扒光：

• $\mathcal{H}[T(z)]$：这是希尔伯特变换后的新曲线（我们将把它作为二维平面的虚部 /$Y$ 轴）。

• $\tau$ (Tau)：这是一个用来做积分的“替身变量”（哑变量）。可以把它想象成一个拿着放大镜的扫描仪，沿着整条测井曲线从头扫到尾。

• $z$：我们当前正在计算的那个具体深度点。

• $\frac{1}{z - \tau}$：这是核心过滤器。当扫描仪$\tau$ 在目标深度$z$ 的左边（上面地层）时，它是正数；在右边（下面地层）时，它是负数。这相当于对信号进行了一次极度精密的“加权翻转”。

• $\text{P.V.}$ (柯西主值)：这是一个数学安全帽。当扫描仪刚好扫到当前深度$\tau = z$时，分母变成了 0，计算会爆炸。$\text{P.V.}$ 的意思是：“在遇到除以 0 的那个无底洞时，我们在黑洞边缘小心翼翼地跳过去，取极限两边的对称抵消值。”

底层物理动作（它到底干了啥？）：

别管积分了，用一句话总结希尔伯特变换的物理本质：

它把原始测井曲线$T(z)$ 里包含的所有波形，极其精确地、没有任何时间/深度延迟地，在空间上向后生生推移了$90^\circ$（即$\frac{\pi}{2}$）。

如果的原信号是一个$\cos(z)$（余弦波，波峰开始），希尔伯特变换后，它就变成了$\sin(z)$（正弦波，从 0 开始）。

原信号和它的希尔伯特“影子”，永远相差四分之一个周期。

第二层：解析信号 (Analytic Signal) —— 从一维线到二维复平面

影子造好了，现在我们要把它们合体，完成伟大的“数学升维”。

数学公式：

符号含义：

• $Z(z)$：这叫“解析信号”。也就是代码里真正生成的那个变量 analytic_signal。

• $T(z)$：我们的原始旋回信号。我们把它放在复平面的**实轴（横坐标 /$X$ 轴）**上。

• $j$：虚数单位（代表$\sqrt{-1}$）。在这里，不要把它当成虚无缥缈的数字，要把它当成**“垂直于$X$ 轴的另一个维度”**！

• $\mathcal{H}[T(z)]$：刚刚算出来的、相移了$90^\circ$ 的影子信号。我们把它放在复平面的**虚轴（纵坐标 /$Y$ 轴）**上。

发生了什么物理质变？

在没有执行这一步之前，的测井数据只是一根在图纸上上下波动的线。

执行完这一步，由于引入了$j$，这根线被直接“拉”出了图纸！

随着深度$z$ 一点点往下走，原信号$T(z)$ 带着它的影子$\mathcal{H}[T(z)]$，在由实轴和虚轴组成的二维复平面上，画出了一个不断绕圈的螺旋。

此时，在这个二维平面里，每一个深度点，都不再是一个单调的数字，而变成了一根指着特定方向的“指针”（向量）。

$\cos(\theta) + j\sin(\theta)$ 绝对不等于$\sin(\theta)$！

它们根本不是同一个维度的东西。$\cos(\theta) + j\sin(\theta)$是一个二维的复数，而$\sin(\theta)$ 只是一个一维的实数。

我说“催生出了$\sin(\theta)$”，真正的意思是：我们利用$\cos(\theta) + j\sin(\theta)$ 这个复数结构，成功地把$\sin(\theta)$ 当作一个“配件”给制造出来，并挂在了$j$ 的后面！

第三层：几何提取 —— 毫不费力地拿走振幅与相位

既然每一个深度点$z$ 现在都变成了一根指针，那么这根指针天然就拥有两个任何人都能看懂的几何属性：长度 和 角度。

根据复数的欧拉公式极坐标表示：

这就是代码里紧接着提取的那两个核心物理量：

1. 瞬时振幅 (Amplitude) —— “指针的长度”

代码：amplitude = np.abs(analytic_signal)

公式：

底层逻辑：用勾股定理算一下这根指针有多长。

地质意义：这就是测井曲线在当前深度的包络线宽度。它代表当前这一层旋回的“震荡能量”有多猛。泥岩和砂岩差异越大的地层，这根指针就越长。

瞬时包络线$A(z)$ 是怎么算出来的？（勾股定理）

在代码里，当这两根曲线结合成 analytic_signal（解析信号）后，CPU 其实在脑海里画了一个直角坐标系。

• 把真实的信号$T(z)$，放在横轴（X 轴）上。

• 把影子信号$H[T(z)]$，放在纵轴（Y 轴）上。

此时，对于任何一个具体的深度点$z$，我们都得到了一个坐标点$(X, Y)$。

如果从原点$(0,0)$ 画一个箭头指向这个坐标点，就得到了一根向量（指针）。

什么是瞬时包络线$A(z)$？

它就是这根指针的绝对长度！

怎么算长度？用初中数学的勾股定理：

• 物理意义：因为$T(z)$ 和$H[T(z)]$ 刚好错开了$90^\circ$（一个大一个小），所以它们的平方和开根号，刚好填补了波形的坑洼，求出的$A(z)$ 能够完美、紧密地贴合在原始信号$T(z)$ 震荡的最外延。这就是包络线的来历！它代表了当前这一个瞬间，地层波动的最大潜在能量。

2. 瞬时相位 (Phase) —— “指针的角度”

代码：phase = np.unwrap(np.angle(analytic_signal))

公式：

底层逻辑：用反三角函数算一下，这根指针现在与横坐标（实轴）的夹角是多少度。

地质意义：这就是我们梦寐以求的“生命周期刻度”！

• 如果指针指在$0^\circ$，意味着原始信号$T(z)$ 刚好处在最大正值（波峰），旋回达到了最泥质的状态。

• 如果指针转到了$180^\circ$（即$\pi$），意味着原始信号$T(z)$ 处于最大负值（波谷），旋回达到了最砂质的状态。

• 如果指针在$90^\circ$ 或$270^\circ$，意味着刚好处于过零点的岩性突变期。

总结这一步

为什么要用 Hilbert 变换？

因为一维的数据无法算角度。

我们利用 Hilbert 变换极其巧妙地构造了一个相移$90^\circ$ 的$Y$ 轴数据，从而在数学空间里撑起了一个“二维复平面”。把原本的上下跳动，完美地映射成了复平面上的圆周运动。

把上下跳动变成了圆周运动后，我们就彻底摆脱了“绝对数值”的束缚，成功提取出了代表相对节律的绝对角度（相位$\theta$）。

第五步：INCOSIN 计算 —— 极其优雅的几何提取

第一层：中学几何的逆袭 —— 为什么除法就等于余弦？

我们在脑海里画一个直角三角形。

把复平面上的那根“指针”当成三角形的斜边。

数学公式还原：

符号物理映射：

• $|Z(z)|$（或者叫$A(z)$）：这是指针的长度，也就是这段代码里的 amplitude（包络线振幅）。它对应直角三角形的斜边。

• $\text{Real}(Z(z))$（或者叫$T(z)$）：这是指针在横轴（X轴）上的投影，也就是第三步拼装出来的原始目标旋回信号。它对应直角三角形的邻边。

初中数学告诉我们：余弦（$\cos\theta$） = 邻边$\div$ 斜边。

所以，要求相位角的余弦，根本不需要去算那个该死的角度$\theta$！直接拿原始信号（邻边），除以它的包络线振幅（斜边），天然就等于$\cos\theta$！

工程意义：用一次极速的除法，代替了昂贵的反三角函数（$\arctan$）和三角函数（$\cos$）计算。这在处理几万个深度的测井点时，速度是碾压级的。

第二层：除法背后的暴力美学 —— 振幅归一化 (Amplitude Normalization)

这不仅仅是个算得快的问题，这个“除号（$\div$）”在地质学上完成了一次暴力洗牌。

想象一下真实的地层：

• 在水动力剧烈变化的滨浅海，一段砂泥岩互层的旋回，GR 值可能在 40 到 150 之间剧烈震荡（振幅$A=55$）。

• 在水体平静的深湖区，一段纯泥岩内部的微弱气候旋回，GR 值可能只在 110 到 115 之间微微蠕动（振幅$A=2.5$）。

如果看原始曲线$T(z)$，那个 5 API 的微弱波动就像一条小虫子，肉眼根本看不见，完全被大波动碾压了。

除以振幅$A(z)$ 究竟干了啥？

就像是上帝伸出了手：

• 对着那个振幅 55 的大波浪说：“给我除以 55，最高只能到 1，最低只能到 -1。”

• 对着那个振幅 2.5 的小虫子说：“也给我除以 2.5，强行拉伸，最高也要到 1，最低也要到 -1。”

结果：不管原本的能量有多大、波动有多小，只要经过这一除，所有的旋回波幅，被极其蛮横地、绝对地锁死在了$[-1, 1]$ 这个区间内！ 那个原本看不见的小虫子，现在和惊涛骇浪一样清晰可见。

第三层：物理意义的彻底重塑 —— 纯粹的节律刻度尺

经过这一步，INCOSIN 曲线的物理意义发生了质的飞跃。

它不再代表 GR 放射性有多强（API 彻底被干掉了），它现在变成了一把纯粹标定“相对节律生命周期”的绝对刻度尺。

看这把尺子上的刻度：

• 当$\text{INCOSIN} = 1$ 时：代表$\theta = 0^\circ$（指针平躺在$X$ 轴正半轴）。这意味着不管当前地层的绝对 GR 是多少，它都已经爬到了这个局部旋回的最顶峰。地质上，这绝对对应着局部海侵最大面、最大泥质含量段、或是洪泛期。

• 当$\text{INCOSIN} = -1$ 时：代表$\theta = 180^\circ$（指针平躺在$X$ 轴负半轴）。这意味着它跌到了局部旋回的最谷底。地质上，这对应着局部水体最浅、砂质最纯、冲刷面或是层序边界。

• 当$\text{INCOSIN} = 0$ 时：代表$\theta = 90^\circ$ 或$270^\circ$（指针垂直）。余弦为 0，意味着它刚好越过均值线，地质上对应着相变的快速转换期（正从泥变砂，或从砂变泥）。

这就是整个算法。

在一堆混沌的、参差不齐的、充满噪音的自然伽马曲线中，通过提取 INCOSIN，生生造出了一台地质节拍器。它极其精准地在每一个深度点上敲击：这里是波峰，这里是波谷，这里是转换点。

第六步：逆向插值

对应代码：interpolate_back_to_original_depth(depth_original, depth_equal, incosin_eq)

底层逻辑：

第一步中，我们为了满足 EMD 和 Hilbert 的数学要求，凭空捏造了一个等距的$z_{eq}$ 深度轴。

现在我们在这个理想轴上算出了完美的$\text{INCOSIN}_{eq}$。但在现实应用中，我们需要把它填回那个因为仪器拉扯而不等距的原始深度列$z_{raw}$ 旁边。

数学动作：

再次使用线性插值（或 Pchip）：

这一步保证了输出的 CSV 文件中，INCOSIN 值能和原始的 GR、计算出的 INPEFA 在同一个深度点上完美对齐，毫无错位。

终极总结

这段 INCOSIN 代码，本质上是在做一套精密的地质信号分离手术：

1. $Pchip$ 强行等深（稳住频率标尺）

2. $\text{EMD}$ 逐层减均值（筛掉干扰频率）

3. $\sum \text{IMF}$ 拼装目标（提取旋回波形）

4. $\text{Hilbert}$ 积分相移$90^\circ$（一维升二维构建复向量）

5. $T(z)/A(z)$ 几何降维（除以振幅，榨出绝对的$[-1,1]$ 纯余弦节律）