Jelajahi Sumber

记录→npz 映射建成(它就是产品 spectra_meta.parquet) + 新增 spectra_lookup.py; 更正上一条的错判

用户令"把记录→npz 映射建起来" —— 实做之后发现**映射根本不用建**:
`scripts/rudong_tcm_spectra.py` 落盘时已同时写出 `spectra_meta.parquet`(每行一条谱:
turbine/sensor/meas_name/trigger_time/**shard/shard_row**/x_offset/x_delta), 取数配方就写在
该脚本头部: v = npz["values"][shard_row]; x = x_offset + arange(len(v)) * x_delta。
⇒ 更正上一轮 §7.5 里"要靠枚举顺序去对"的说法(那是我没读落盘侧就下的判断)。

· 新增 scripts/spectra_lookup.py: 把查表封成一条命令, 另带 --stats 看整窗谱库概况。
  实测(w0316): 谱件记录 420,742 条 · 分片 1,701 个 · 8 个测点 · 38 台 · 13 类 meas;
  低频谱件 FFT_4_Tr 是 0.01 Hz/线(401 点 ⇒ 0~4 Hz), 正是 1P(0.223 Hz) 能分辨的那一档。
  单条实测: WTG01/Main_bearing_rear/FFT_4_Tr @2026-03-31 23:04:15 →
  spectra/p00/FFT_4_Tr_n401_0000.npz 第 1 行, 401 点, 0.000~4.000 Hz; @0.223 Hz(最近线 0.220 Hz) = 0.000570981。
· 过程中修掉一次自己的解析错: shard 是**相对 <窗>/spectra/** 的(meta 记 p00/xxx.npz, 文件在
  <窗>/spectra/p00/xxx.npz); 第一版按"meta 的父目录"解析 → 误报"分片不在位"。已改并写进注释。
· §7.5 同步更正: 三步路径的第①步标记为**已完成**, 并注明 1P 那批(2026-01-27~02-01, w0127)的
  谱库当前**不在盘上**(只有 w0316 有谱), 要复算 1P 须先补 1 月窗的 index+spectra(原始件在
  data/raw/如东/windcms/, 按既有实测 index≈214s + spectra≈305s)。

位置不变: 仍然只差"理论频率附近怎么取观测峰"那一句口径 —— 现在测试台更完备了(记录→谱 → 频率轴 → 幅值, 一条命令)。
zhouyang.xie 3 minggu lalu
induk
melakukan
7313e4e4fe

+ 2 - 2
docs/振动六层链_接口规格与缺口_v0.1.md

@@ -219,11 +219,11 @@ MAD=0.80209 → z=−0.3116 · IQR/1.349=1.46899 → z=−0.1701 ·
 |---|---|
 | 覆盖范围 | 只两个测点:`Main_bearing_front` / `Main_bearing_rear`;662 行;时间 2026-01-27 11:35 ~ 2026-02-01 23:16(**w0127 首窗**,其索引在 `m5/tcm_index.parquet`,不在 `windows/w0127/`) |
 | 转速来源 | ★ 样件行的 `rpm` **不是**同一记录里 `Kurtosis` 行的 rpm:样件 `rpm=1601.0371`,而索引里 `trigger_time` 同一条的 `Kurtosis` 行 `rpm=1604.2537`。⇒ 1P 分析用的是**该记录自己的转速**(`RpmProfile` / `rawRPM` 那类),不是标量指标行带的转速。这条搞错,`f1P` 就对不上(速比 119.75 本身不变) |
-| 谱文件命名 | `spectra/p<NN>/<meas_name>_n<lines>_<seq>.npz`(例 `FFT_10000_Tr_n401_0000.npz`)—— 谱件**不带源件 UUID**,所以"某条记录 ↔ 某个 npz"要靠 `scripts/rudong_tcm_spectra.py` 的**枚举顺序**去对(该脚本在包内,规则可读) |
+| 谱文件命名 | ★**更正**:`spectra/p<NN>/<meas_name>_n<lines>_<seq>.npz`(例 `FFT_4_Tr_n401_0000.npz`);"记录 ↔ npz"**不需要靠枚举顺序去猜** —— `rudong_tcm_spectra.py` 落盘时已同时写出 `spectra_meta.parquet`(每行一条谱:turbine/sensor/meas_name/trigger_time/**shard/shard_row**/x_offset/x_delta),取数配方就是该脚本头部那句 `v = npz["values"][shard_row]; x = x_offset + arange(len(v))*x_delta`。本包新增 `scripts/spectra_lookup.py` 把它封成一条命令 |
 | 谱本身 | npz 内是幅值数组;频率轴由索引的 `x_offset`/`x_delta`/`lines` 给出(例 `FFT_16000: lines=400, x_offset=0, x_delta=40` ⇒ 40 Hz/线)⇒ 频率轴可算;`f1P=0.223 Hz` 落在**第 0 条线**上 —— 这正是"1P 在低频谱里极难分辨"的原因,也解释了为什么这一族要专门做(样件 `energy_pct` 只有 0.38%) |
 
 **下一次接手的最短路径**(三步,都在包内可做):
-① 读 `scripts/rudong_tcm_spectra.py` 的枚举顺序,建立"记录 → npz"映射;
+① ~~读枚举顺序建映射~~ **已完成**:映射本来就是产品(`spectra_meta.parquet`),本包已封装为 `scripts/spectra_lookup.py`(实测 `WTG01/Main_bearing_rear/FFT_4_Tr` → `spectra/p00/FFT_4_Tr_n401_0000.npz` 第 1 行,0.01 Hz/线,401 点);
 ② 用 `x_offset + i*x_delta` 重建频率轴,在 `f1P` 附近按候选口径(最近线 / ±2% 带内取最大 / 抛物线插值)取幅值,与样件 `a1P` 对拍;
 ③ `a1P` 一旦对上,`struct`(疑似 `a1P / fleet` 归一分)、`energy_pct`、`g1/g2/g4` 闸门与 `verdict` 有望按同样办法逐个对拍。
 

+ 134 - 0
scripts/spectra_lookup.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+r"""谱库查表:**记录 → npz 分片**(2026-09-17 用户令"把映射建起来")。
+
+## 结论先说:这个映射**已经是产品的一部分**,不需要重建
+
+`scripts/rudong_tcm_spectra.py` 落盘时同时写了 `spectra_meta.parquet`,每行一条谱:
+
+    turbine · sensor · meas_name · trigger_time · shard · shard_row · x_offset · x_delta · …
+
+取数配方(该脚本头部原文):
+
+    v = np.load(<store>/<shard> + '.npz')['values'][shard_row]
+    x = x_offset + arange(len(v)) * x_delta
+
+所以"某条记录对应哪个 npz 的哪一行"**直接查表就有**;本器把它封装成一条命令,
+让"在特征频率附近取幅值"这类下一步工作只需几行代码。
+
+## 用法
+
+    python scripts/spectra_lookup.py --window w0316 --turbine WTG01 --sensor Gear_HS_generator_side \
+        --meas FFT_10000_Tr [--at 0.223] [--list]
+    python scripts/spectra_lookup.py --window w0316 --stats        # 各 meas 的谱件/记录数、频率轴分辨率
+
+退出码: 0 找到 · 2 窗/记录/谱件缺(会指明缺哪一件)
+"""
+from __future__ import annotations
+
+import argparse
+import pathlib
+import sys
+
+import numpy as np
+import pandas as pd
+
+ROOT = pathlib.Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+from src import paths as P                                              # noqa: E402
+
+
+def meta_path(win: str) -> pathlib.Path:
+    """window → spectra_meta.parquet(两份落点都认:窗根与 spectra/ 下)。"""
+    m5 = P.m5()
+    if win == 'w0127':                       # 首窗特例: 索引在 m5 根, 谱在 m5/spectra
+        cands = [m5 / 'spectra' / 'spectra_meta.parquet', m5 / 'spectra_meta.parquet']
+    else:
+        w = m5 / 'windows' / win
+        cands = [w / 'spectra_meta.parquet', w / 'spectra' / 'spectra_meta.parquet']
+    for c in cands:
+        if c.is_file():
+            return c
+    return cands[0]
+
+
+def store_of(mp: pathlib.Path) -> pathlib.Path:
+    """npz 分片所在根:meta 里 shard 是**相对该根**的(含 p<NN>/ 子目录)。
+
+    ★ 两份 meta 落点对应同一个根:`<窗>/spectra_meta.parquet` 与 `<窗>/spectra/spectra_meta.parquet`
+      的 shard 都是相对 **`<窗>/spectra/`** 的(实测:meta 记 `p00/FFT_4_Tr_n401_0000.npz`,
+      文件在 `<窗>/spectra/p00/…`)。第一版按 "meta 的父目录" 解析,于是报"分片不在位"。
+    """
+    return mp.parent if mp.parent.name == 'spectra' else (mp.parent / 'spectra')
+
+
+def load(win: str, turbine: str, sensor: str, meas: str, trigger: str | None = None):
+    mp = meta_path(win)
+    if not mp.is_file():
+        raise SystemExit(f'[X] 没有谱库元数据: {P.rel(mp)} —— 该窗的谱还没摄入 '
+                         f'(scripts/rudong_tcm_spectra.py);注意 w0127(1 月窗) 的谱库当前不在盘上。')
+    m = pd.read_parquet(mp)
+    sel = m[(m['turbine'] == turbine) & (m['sensor'] == sensor) & (m['meas_name'] == meas)]
+    if sel.empty:
+        raise SystemExit(f'[X] 表里没有 {turbine}/{sensor}/{meas} 的谱')
+    if trigger:
+        t = pd.to_datetime(sel['trigger_time'], errors='coerce')
+        sel = sel.loc[[(t - pd.to_datetime(trigger)).abs().idxmin()]]
+    r = sel.iloc[0]
+    npz = store_of(mp) / str(r['shard'])
+    npz = npz.with_suffix('.npz') if npz.suffix != '.npz' else npz
+    if not npz.is_file():
+        raise SystemExit(f'[X] 分片不在位: {P.rel(npz)}(meta 记的 shard={r["shard"]})')
+    with np.load(npz, allow_pickle=False) as z:
+        v = z['values'][int(r['shard_row'])]
+    x = float(r['x_offset']) + np.arange(len(v)) * float(r['x_delta'])
+    return dict(row=r, npz=npz, v=np.asarray(v, dtype=float), x=x,
+                x_offset=float(r['x_offset']), x_delta=float(r['x_delta']))
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description='谱库查表: 记录 → npz 分片 (取数配方见脚本头)')
+    ap.add_argument('--window', default='w0316')
+    ap.add_argument('--turbine', default=None)
+    ap.add_argument('--sensor', default=None)
+    ap.add_argument('--meas', default=None)
+    ap.add_argument('--trigger', default=None, help='指定 trigger_time(不指定取该组合第一条)')
+    ap.add_argument('--at', type=float, default=None, help='看这个频率(Hz)处的幅值(最近线)')
+    ap.add_argument('--stats', action='store_true', help='只打印该窗谱库的概况')
+    a = ap.parse_args()
+    mp = meta_path(a.window)
+    if not mp.is_file():
+        print(f'[X] 谱库元数据不在位: {P.rel(mp)}')
+        print(f'    该窗的谱还没摄入。重新摄入: python scripts/rudong_tcm_index.py --root <decode 目录> '
+              f'--out <窗>/index.parquet 然后 python scripts/rudong_tcm_spectra.py --root <同> --out <窗>/spectra')
+        return 2
+    m = pd.read_parquet(mp)
+    if a.stats or not a.turbine:
+        print(f'== 谱库概况 · 窗 {a.window} · {P.rel(mp)} ==')
+        print(f'   谱件记录 {len(m)} 条 · 分片 {m["shard"].nunique()} 个 · 测点 {m["sensor"].nunique()} 个'
+              f' · 机组 {m["turbine"].nunique()} 台')
+        g = m.groupby('meas_name').agg(条数=('shard_row', 'size'), 分片=('shard', 'nunique'),
+                                       dx=('x_delta', 'first'), x0=('x_offset', 'first'))
+        print(g.sort_values('条数', ascending=False).head(18).to_string())
+        print(f'   取数配方: v = npz["values"][shard_row]; x = x_offset + arange(len(v))*x_delta')
+        return 0
+    got = load(a.window, a.turbine, a.sensor, a.meas, a.trigger)
+    r = got['row']
+    print(f'== 记录 → npz ==')
+    print(f'   {r["turbine"]} / {r["sensor"]} / {r["meas_name"]} @ {r["trigger_time"]}')
+    print(f'   分片 {P.rel(got["npz"])} · 第 {int(r["shard_row"])} 行 · {len(got["v"])} 点'
+          f' · x: {got["x_offset"]} + i*{got["x_delta"]} Hz ⇒ {got["x"][0]:.3f} ~ {got["x"][-1]:.3f} Hz')
+    print(f'   幅值: min={got["v"].min():.6g} max={got["v"].max():.6g} 均值={got["v"].mean():.6g}')
+    if a.at is not None:
+        i = int(np.argmin(np.abs(got['x'] - a.at)))
+        print(f'   @{a.at} Hz(最近线 {got["x"][i]:.3f} Hz, 第 {i} 条): {got["v"][i]:.6g}')
+    return 0
+
+
+if __name__ == '__main__':
+    for _s in (sys.stdout, sys.stderr):
+        try:
+            _s.reconfigure(errors='replace')
+        except Exception:
+            pass
+    sys.exit(main())