实验室多舱环境采集：TCP/IP以太网温湿度传感器时钟同步与数据校准

┌─────────────────────────────────────────────────────────────┐
│                 时钟源层                                     │
│  GPS/北斗 原子钟 / 本地铷钟（可选）                          │
│       │                                                     │
│  ┌────▼──────────┐                                         │
│  │ PTP Grandmaster│ (支持 PTP + NTP 同时输出)               │
│  │ (Meinberg/Orolia│                                        │
│  └────┬──────────┘                                         │
├───────┼─────────────────────────────────────────────────────┤
│       │                                                     │
│  ┌────▼──────────┐      ┌──────────────┐                   │
│  │ 核心交换机     │      │ NTP 服务器    │ (从 PTP 同步)     │
│  │ (边界时钟)     │      └──────┬───────┘                   │
│  └────┬──────────┘             │                           │
│       │                         │                           │
│  ┌────▼──────────┐      ┌──────▼───────┐                   │
│  │ 接入交换机     │      │ 普通设备      │ (PC/服务器/NTP)   │
│  │ (透明时钟)     │      └──────────────┘                   │
│  └────┬──────────┘                                          │
│       │                                                      │
│  ┌────▼──────────┐      ┌──────────────────┐               │
│  │ 传感器（PTP）  │      │ 传感器（无 PTP）  │ (NTP 客户端)  │
│  │ 高端型号支持    │      │ 通过网关打时间戳  │               │
│  └───────────────┘      └──────────────────┘               │
└─────────────────────────────────────────────────────────────┘

"""
sensor_time_config.py - 批量配置传感器 NTP
"""

import asyncio
from pymodbus.client import AsyncModbusTcpClient

# 传感器支持 NTP 配置寄存器（示例，需按厂家点表）
NTP_SERVER_IP_REG = 0x50    # 4个寄存器存 NTP 服务器 IP（4字节）
NTP_ENABLE_REG = 0x54       # 1=启用, 0=禁用
NTP_SYNC_INTERVAL_REG = 0x55 # 同步间隔（秒）

async def config_ntp(devices, ntp_server_ip: str):
    """批量配置 NTP"""
    # IP 转 4 个 16 位寄存器
    octets = [int(x) for x in ntp_server_ip.split(".")]
    ip_regs = [ (octets[0]<<8)|octets[1], (octets[2]<<8)|octets[3] ]
    
    for dev in devices:
        c = AsyncModbusTcpClient(dev["ip"], port=502)
        try:
            await c.connect()
            # 写 NTP 服务器 IP
            await c.write_registers(NTP_SERVER_IP_REG, ip_regs, slave=dev["unit"])
            # 写同步间隔 3600s
            await c.write_register(NTP_SYNC_INTERVAL_REG, 3600, slave=dev["unit"])
            # 启用 NTP
            await c.write_register(NTP_ENABLE_REG, 1, slave=dev["unit"])
            print(f"{dev['ip']}: NTP configured")
        except Exception as e:
            print(f"{dev['ip']}: failed - {e}")
        finally:
            c.close()

# 无 NTP 支持的传感器：记录设备时标与采集端时标的偏移，用于后期补偿
async def measure_clock_offset(dev_ip, samples=10):
    """测量设备时钟偏移"""
    offsets = []
    c = AsyncModbusTcpClient(dev_ip, port=502)
    await c.connect()
    for _ in range(samples):
        recv_ts = time.time()
        rr = await c.read_holding_registers(0x04, 2, slave=1)  # 设备时标寄存器
        if not rr.isError():
            dev_ts = (rr.registers[0] << 16) | rr.registers[1]
            offsets.append(recv_ts - dev_ts)
        await asyncio.sleep(1)
    c.close()
    return sum(offsets) / len(offsets) if offsets else None

"""
timestamping.py - 采集端精确时间戳
"""

import time
import asyncio
from dataclasses import dataclass

@dataclass
class TimestampedReading:
    dev_id: str
    temperature: float
    humidity: float
    device_ts: float       # 设备内部时标（如果有）
    recv_ts: float         # 采集端接收时间（PTP 同步后）
    ingest_ts: float       # 应用层处理时间
    clock_offset: float    # device_ts - recv_ts（用于检测漂移）
    quality: str           # good / clock_drift / device_unsynced

class PreciseCollector:
    def __init__(self, ptp_sync_threshold_ms=10):
        self.ptp_sync_threshold_ms = ptp_sync_threshold_ms
        self.clock_offsets = {}  # dev_id -> rolling average offset
    
    def get_hw_timestamp(self) -> float:
        """获取硬件时间戳（如果网卡支持 SO_TIMESTAMPING）"""
        # Linux 上可以通过 socket 选项获取硬件时间戳
        # 此处简化为系统单调时钟
        return time.time()
    
    def classify_quality(self, dev_id, device_ts, recv_ts) -> str:
        """根据时钟偏移分类数据质量"""
        if device_ts is None:
            return "device_unsynced"
        
        offset_ms = abs(device_ts - recv_ts) * 1000
        if offset_ms > self.ptp_sync_threshold_ms:
            return "clock_drift"
        return "good"
    
    async def on_packet(self, dev_id, raw_temp, raw_hum, device_ts, data):
        recv_ts = self.get_hw_timestamp()
        quality = self.classify_quality(dev_id, device_ts, recv_ts)
        
        reading = TimestampedReading(
            dev_id=dev_id,
            temperature=raw_temp / 10.0,
            humidity=raw_hum / 10.0,
            device_ts=device_ts,
            recv_ts=recv_ts,
            ingest_ts=time.time(),
            clock_offset=(device_ts - recv_ts) if device_ts else None,
            quality=quality
        )
        
        # 更新滚动偏移量
        if device_ts:
            offset = device_ts - recv_ts
            hist = self.clock_offsets.setdefault(dev_id, [])
            hist.append(offset)
            if len(hist) > 100:
                hist.pop(0)
        
        return reading

"""
calibration.py - 偏移校准
"""

import numpy as np
from dataclasses import dataclass
from typing import List, Dict

@dataclass
class CalPoint:
    sensor_val: float
    reference_val: float
    timestamp: float

@dataclass
class CalResult:
    offset: float          # 偏移量
    slope: float = 1.0     # 斜率（单点校准时为 1.0）
    rmse: float = 0.0      # 均方根误差
    sample_count: int = 0
    valid: bool = False

class OffsetCalibrator:
    def __init__(self, stability_threshold=0.1):
        self.stability_threshold = stability_threshold  # 标准器稳定阈值
    
    def calibrate(self, points: List[CalPoint], method="mean") -> CalResult:
        """计算偏移量"""
        if len(points) < 3:
            return CalResult(offset=0.0, valid=False)
        
        sensor_vals = np.array([p.sensor_val for p in points])
        ref_vals = np.array([p.reference_val for p in points])
        
        # 稳定性检查：标准器读数方差应小于阈值
        if np.std(ref_vals) > self.stability_threshold:
            return CalResult(offset=0.0, valid=False)
        
        if method == "mean":
            # 简单平均偏移
            offsets = sensor_vals - ref_vals
            offset = np.mean(offsets)
            rmse = np.sqrt(np.mean(offsets**2))
        elif method == "ols":
            # 最小二乘拟合 y = mx + b，取截距为偏移
            slope, intercept = np.polyfit(sensor_vals, ref_vals, 1)
            offset = intercept
            slope = slope
            rmse = np.sqrt(np.mean((ref_vals - (slope * sensor_vals + intercept))**2))
        else:
            raise ValueError(f"Unknown method: {method}")
        
        return CalResult(
            offset=offset,
            slope=slope if method == "ols" else 1.0,
            rmse=rmse,
            sample_count=len(points),
            valid=True
        )

# 校准结果存储（按设备、按日期、按环境舱）
class CalibrationStore:
    def __init__(self, db_path="calibration.db"):
        self.db_path = db_path
        self._init_db()
    
    def _init_db(self):
        import sqlite3
        conn = sqlite3.connect(self.db_path)
        conn.execute("""
            CREATE TABLE IF NOT EXISTS calibrations (
                dev_id TEXT,
                chamber_id TEXT,
                cal_date TEXT,
                cal_type TEXT,
                offset REAL,
                slope REAL,
                rmse REAL,
                sample_count INTEGER,
                operator TEXT,
                cert_ref TEXT,
                PRIMARY KEY (dev_id, cal_date)
            )
        """)
        conn.commit()
        conn.close()
    
    def save(self, dev_id, chamber_id, result: CalResult, operator, cert_ref):
        import sqlite3
        conn = sqlite3.connect(self.db_path)
        conn.execute(
            "INSERT OR REPLACE INTO calibrations VALUES (?,?,?,?,?,?,?,?,?,?)",
            (dev_id, chamber_id, time.strftime("%Y-%m-%d %H:%M:%S"),
             "offset", result.offset, result.slope, result.rmse,
             result.sample_count, operator, cert_ref)
        )
        conn.commit()
        conn.close()
    
    def get_latest(self, dev_id) -> CalResult | None:
        """获取最新有效校准"""
        import sqlite3
        conn = sqlite3.connect(self.db_path)
        row = conn.execute(
            "SELECT offset, slope, rmse, sample_count FROM calibrations "
            "WHERE dev_id=? ORDER BY cal_date DESC LIMIT 1", (dev_id,)
        ).fetchone()
        conn.close()
        if row:
            return CalResult(offset=row[0], slope=row[1], rmse=row[2],
                           sample_count=row[3], valid=True)
        return None

"""
cross_calibration.py - 多舱交叉校准
"""

import numpy as np
from typing import Dict, List

class CrossCalibrator:
    def __init__(self, topology: Dict[str, List[str]]):
        """
        topology: {chamber_id: [adjacent_chamber_ids]}
        描述舱室间的物理连接关系
        """
        self.topology = topology
    
    def detect_drift(self, readings: Dict[str, Dict[str, float]], 
                     window_minutes=60) -> Dict[str, float]:
        """
        readings: {chamber_id: {"temp": float, "hum": float}}
        返回每个舱室的漂移评分（0=正常，>1=可能漂移）
        """
        drift_scores = {}
        
        for chamber, vals in readings.items():
            neighbors = self.topology.get(chamber, [])
            if not neighbors:
                drift_scores[chamber] = 0.0
                continue
            
            # 与相邻舱室比较
            neighbor_temps = [readings[n]["temp"] for n in neighbors if n in readings]
            neighbor_hums = [readings[n]["hum"] for n in neighbors if n in readings]
            
            if not neighbor_temps:
                drift_scores[chamber] = 0.0
                continue
            
            # 温度偏差（与邻居平均值的差异）
            temp_deviation = abs(vals["temp"] - np.mean(neighbor_temps))
            # 湿度偏差
            hum_deviation = abs(vals["hum"] - np.mean(neighbor_hums))
            
            # 漂移评分（基于物理合理性）
            # 相邻舱室温差通常 < 2℃，湿度差 < 10%RH（除非有独立控制）
            temp_score = temp_deviation / 2.0
            hum_score = hum_deviation / 10.0
            
            drift_scores[chamber] = max(temp_score, hum_score)
        
        return drift_scores
    
    def consensus_calibration(self, readings: Dict[str, Dict[str, float]],
                               reference_chambers: List[str]) -> Dict[str, float]:
        """
        以参考舱室（已知校准的）为基准，计算其他舱室的修正量
        """
        corrections = {}
        
        # 参考舱室的平均值作为基准
        ref_temps = [readings[r]["temp"] for r in reference_chambers if r in readings]
        ref_hums = [readings[r]["hum"] for r in reference_chambers if r in readings]
        
        if not ref_temps:
            return corrections
        
        base_temp = np.mean(ref_temps)
        base_hum = np.mean(ref_hums)
        
        for chamber, vals in readings.items():
            if chamber in reference_chambers:
                corrections[chamber] = {"temp_offset": 0.0, "hum_offset": 0.0}
            else:
                corrections[chamber] = {
                    "temp_offset": base_temp - vals["temp"],
                    "hum_offset": base_hum - vals["hum"]
                }
        
        return corrections

"""
apply_calibration.py - 实时数据校准
"""

class CalibratedCollector:
    def __init__(self, cal_store, cross_calibrator=None):
        self.cal_store = cal_store
        self.cross_calibrator = cross_calibrator
        self.active_corrections = {}  # dev_id -> (offset, slope)
        self.cross_corrections = {}   # chamber_id -> (temp_offset, hum_offset)
    
    def load_calibrations(self, dev_ids):
        """加载所有设备的校准参数"""
        for dev_id in dev_ids:
            cal = self.cal_store.get_latest(dev_id)
            if cal and cal.valid:
                self.active_corrections[dev_id] = (cal.offset, cal.slope)
    
    def apply(self, reading: TimestampedReading) -> TimestampedReading:
        """应用校准"""
        corr = self.active_corrections.get(reading.dev_id)
        if corr:
            offset, slope = corr
            reading.temperature = reading.temperature * slope + offset
            # 湿度通常只做偏移，不做斜率
            # reading.humidity = reading.humidity * hum_slope + hum_offset
        
        return reading
    
    def apply_cross_correction(self, chamber_id, temp, hum):
        """应用交叉校准修正"""
        corr = self.cross_corrections.get(chamber_id)
        if corr:
            temp += corr.get("temp_offset", 0.0)
            hum += corr.get("hum_offset", 0.0)
        return temp, hum

"""
quality_and_audit.py - 数据质量标记与审计
"""

from enum import IntFlag
from dataclasses import dataclass, field

class DataQuality(IntFlag):
    GOOD = 0x00
    CLOCK_DRIFT = 0x01
    DEVICE_UNSYNCED = 0x02
    CALIBRATION_EXPIRED = 0x04
    SENSOR_FAULT = 0x08
    OUT_OF_RANGE = 0x10
    INTERPOLATED = 0x20

@dataclass
class AuditedReading:
    dev_id: str
    chamber_id: str
    temperature: float
    humidity: float
    timestamp: float
    quality: DataQuality
    calibration_offset: float
    calibration_date: str
    operator: str = ""

class AuditLogger:
    def __init__(self, log_path="audit.log"):
        self.log_path = log_path
    
    def log_calibration(self, dev_id, chamber_id, operator, result, cert_ref):
        """记录校准操作"""
        import json
        entry = {
            "event": "calibration",
            "timestamp": time.time(),
            "dev_id": dev_id,
            "chamber_id": chamber_id,
            "operator": operator,
            "offset": result.offset,
            "slope": result.slope,
            "rmse": result.rmse,
            "sample_count": result.sample_count,
            "cert_ref": cert_ref
        }
        with open(self.log_path, "a") as f:
            f.write(json.dumps(entry) + "\n")
    
    def log_quality_change(self, dev_id, old_quality, new_quality, reason):
        """记录数据质量变化"""
        import json
        entry = {
            "event": "quality_change",
            "timestamp": time.time(),
            "dev_id": dev_id,
            "old_quality": str(old_quality),
            "new_quality": str(new_quality),
            "reason": reason
        }
        with open(self.log_path, "a") as f:
            f.write(json.dumps(entry) + "\n")