IoT设备群体分析实战：用DolphinDB实现从异常诊断到健康评分

// 创建数据库：日期VALUE分区 + 设备HASH分区
dbName = "dfs://device_health"
if (existsDatabase(dbName)) dropDatabase(dbName)
db1 = database("", VALUE, 2024.01.01..2024.12.31)
db2 = database("", HASH, [SYMBOL, 10])
db = database(dbName, COMPO, [db1, db2])

schema = table(1:0,
    `ts`deviceId`current`vibration`temperature`rpm`,
    [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, DOUBLE, DOUBLE])

db.createPartitionedTable(schema, "sensor", `ts`deviceId)

def generateDeviceData(deviceCount, days, anomalyDevices, anomalyLevel) {
    // anomalyDevices: 异常设备列表，如 ["DEV_003", "DEV_017"]
    // anomalyLevel: 劣化程度，如 1.3 表示比正常高出30%
    devices = "DEV_" + lpad(string(1..deviceCount), 3, "0")
    totalSecs = days * 86400
    t = table(totalSecs * deviceCount:0,
        `ts`deviceId`current`vibration`temperature`rpm`,
        [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, DOUBLE, DOUBLE])

    startTs = tempParse(timestamp, "2024.06.01 00:00:00.000")

    for (sec in 0..(totalSecs - 1)) {
        ts = startTs + sec * 1000
        for (dev in devices) {
            // 正常基线
            baseCurrent = 12.5 + 2.0 * sin(double(sec % 86400) / 86400.0 * 2 * pi)
            baseVibration = 1.2 + norm(0, 0.15)
            baseTemp = 55.0 + 5.0 * sin(double(sec % 86400) / 86400.0 * 2 * pi)
            baseRpm = 1480 + norm(0, 10)

            // 异常设备：各项指标偏高
            isAnomaly = anomalyDevices.find(dev) >= 0
            factor = iif(isAnomaly, anomalyLevel, 1.0)

            insert into t values(
                ts, dev,
                baseCurrent * factor + norm(0, 0.5),
                baseVibration * factor + norm(0, 0.05),
                baseTemp * factor + norm(0, 1.0),
                baseRpm / factor + norm(0, 5)
            )
        }
    }
    return t
}

// 生成7天数据：50台设备，DEV_003和DEV_017有轻度劣化
data = generateDeviceData(50, 7, ["DEV_003", "DEV_017"], 1.3)
loadTable("dfs://device_health", "sensor").append!(data)

// 提取每台设备7天的统计特征
select deviceId,
    // 电流特征
    avg(current) as avgCurrent,
    std(current) as stdCurrent,
    percentile(current, 95) as p95Current,
    max(current) - min(current) as rangeCurrent,

    // 振动特征
    avg(vibration) as avgVib,
    std(vibration) as stdVib,
    percentile(vibration, 99) as p99Vib,
    max(vibration) as maxVib,

    // 温度特征
    avg(temperature) as avgTemp,
    std(temperature) as stdTemp,
    max(temperature) as maxTemp,

    // 转速特征
    avg(rpm) as avgRpm,
    std(rpm) as stdRpm,
    min(rpm) as minRpm
from loadTable("dfs://device_health", "sensor")
group by deviceId

// 每台设备的日内小时统计
select deviceId,
       hour(ts) as hour,
       avg(vibration) as hourlyAvgVib,
       avg(current) as hourlyAvgCurrent,
       avg(temperature) as hourlyAvgTemp
from loadTable("dfs://device_health", "sensor")
group by deviceId, hour(ts)

// 变异系数 = 标准差 / 均值
select deviceId,
    std(current) / avg(current) as cvCurrent,
    std(vibration) / avg(vibration) as cvVibration,
    std(temperature) / avg(temperature) as cvTemperature,
    std(rpm) / avg(rpm) as cvRpm
from loadTable("dfs://device_health", "sensor")
group by deviceId

// 各维度排名：1=最好，50=最差
select deviceId,
    avgVib, rank(avgVib, true) as vibRank,
    avgCurrent, rank(avgCurrent, true) as currentRank,
    maxTemp, rank(maxTemp, true) as tempRank,
    avgRpm, rank(avgRpm, false) as rpmRank
from (
    select deviceId,
        avg(vibration) as avgVib,
        avg(current) as avgCurrent,
        max(temperature) as maxTemp,
        avg(rpm) as avgRpm
    from loadTable("dfs://device_health", "sensor")
    group by deviceId
)
order by vibRank desc

// 每台设备每小时的平均振动值，转为矩阵（行=小时，列=设备）
exec avg(vibration) as avgVib
from loadTable("dfs://device_health", "sensor")
pivot by hour(ts), deviceId

// 提取每台设备的分钟级振动序列
vibMatrix = exec avg(vibration)
    from loadTable("dfs://device_health", "sensor")
    pivot by bar(ts, 1m), deviceId

// 计算设备间的相关系数矩阵
corrMatrix = cross(corr, vibMatrix, vibMatrix)

// 计算每个特征在截面上的Z-Score
select deviceId,
    avgVib,
    (avgVib - avg(avgVib)) / std(avgVib) as zVib,
    avgCurrent,
    (avgCurrent - avg(avgCurrent)) / std(avgCurrent) as zCurrent,
    maxTemp,
    (maxTemp - avg(maxTemp)) / std(maxTemp) as zTemp,
    avgRpm,
    (avgRpm - avg(avgRpm)) / std(avgRpm) as zRpm
from (
    select deviceId,
        avg(vibration) as avgVib,
        avg(current) as avgCurrent,
        max(temperature) as maxTemp,
        avg(rpm) as avgRpm
    from loadTable("dfs://device_health", "sensor")
    group by deviceId
)
context by NULL

// 用IQR方法识别异常设备
select deviceId, avgVib,
    iqr(avgVib) as vibIQR,
    percentile(avgVib, 25) - 1.5 * iqr(avgVib) as vibLower,
    percentile(avgVib, 75) + 1.5 * iqr(avgVib) as vibUpper,
    iif(avgVib > percentile(avgVib, 75) + 1.5 * iqr(avgVib), "异常偏高",
        iif(avgVib < percentile(avgVib, 25) - 1.5 * iqr(avgVib), "异常偏低", "正常")) as vibStatus
from (
    select deviceId, avg(vibration) as avgVib
    from loadTable("dfs://device_health", "sensor")
    group by deviceId
)
context by NULL

// 多指标联合判定：至少2个指标Z-Score超过阈值
select deviceId,
    zVib, zCurrent, zTemp, zRpm,
    (abs(zVib) > 2 ? 1 : 0) +
    (abs(zCurrent) > 2 ? 1 : 0) +
    (abs(zTemp) > 2 ? 1 : 0) +
    (abs(zRpm) > 2 ? 1 : 0) as anomalyCount,
    iif(
        (abs(zVib) > 2 ? 1 : 0) +
        (abs(zCurrent) > 2 ? 1 : 0) +
        (abs(zTemp) > 2 ? 1 : 0) +
        (abs(zRpm) > 2 ? 1 : 0) >= 2,
        "疑似异常", "正常"
    ) as status
from (
    select deviceId,
        avgVib, (avgVib - avg(avgVib)) / std(avgVib) as zVib,
        avgCurrent, (avgCurrent - avg(avgCurrent)) / std(avgCurrent) as zCurrent,
        maxTemp, (maxTemp - avg(maxTemp)) / std(maxTemp) as zTemp,
        avgRpm, (avgRpm - avg(avgRpm)) / std(avgRpm) as zRpm
    from (
        select deviceId,
            avg(vibration) as avgVib,
            avg(current) as avgCurrent,
            max(temperature) as maxTemp,
            avg(rpm) as avgRpm
        from loadTable("dfs://device_health", "sensor")
        group by deviceId
    )
    context by NULL
)
where (abs(zVib) > 2 ? 1 : 0) +
      (abs(zCurrent) > 2 ? 1 : 0) +
      (abs(zTemp) > 2 ? 1 : 0) +
      (abs(zRpm) > 2 ? 1 : 0) >= 2
order by anomalyCount desc

// Step 1: 计算原始指标
rawMetrics = select deviceId,
    avg(vibration) as avgVib,
    max(vibration) as maxVib,
    avg(temperature) as avgTemp,
    max(temperature) as maxTemp,
    avg(current) as avgCurrent,
    avg(rpm) as avgRpm,
    avg(current) / avg(rpm) * 1000 as currentPerRpm,
    std(vibration) / avg(vibration) as cvVib,
    std(current) / avg(current) as cvCurrent,
    std(temperature) / avg(temperature) as cvTemp
from loadTable("dfs://device_health", "sensor")
group by deviceId

// Step 2: 截面排名归一化（rank越靠前=越好=分数越高）
rankedMetrics = select deviceId,
    (1.0 - rank(avgVib, true) / count(avgVib)) * 100 as vibScore,
    (1.0 - rank(maxTemp, true) / count(maxTemp)) * 100 as tempScore,
    (1.0 - rank(currentPerRpm, true) / count(currentPerRpm)) * 100 as effScore,
    (1.0 - rank(cvVib + cvCurrent + cvTemp, true) / count(cvVib)) * 100 as stabScore
from rawMetrics
context by NULL

// Step 3: 加权求和
healthScore = select deviceId,
    vibScore, tempScore, effScore, stabScore,
    vibScore * 0.30 + tempScore * 0.25 + effScore * 0.25 + stabScore * 0.20 as healthScore,
    iif(vibScore * 0.30 + tempScore * 0.25 + effScore * 0.25 + stabScore * 0.20 >= 80, "优秀",
        iif(vibScore * 0.30 + tempScore * 0.25 + effScore * 0.25 + stabScore * 0.20 >= 60, "良好",
            iif(vibScore * 0.30 + tempScore * 0.25 + effScore * 0.25 + stabScore * 0.20 >= 40, "一般", "警告"))) as grade
from rankedMetrics
order by healthScore asc

// 按天计算每台设备的健康评分趋势
dailyScore = select deviceId, date(ts) as day,
    avg(vibration) as avgVib,
    max(temperature) as maxTemp,
    avg(current) / avg(rpm) * 1000 as currentPerRpm,
    std(vibration) / avg(vibration) as cvVib
from loadTable("dfs://device_health", "sensor")
group by deviceId, date(ts)

// 计算每天的截面排名并生成评分
select deviceId, day,
    avgVib,
    (1.0 - rank(avgVib, true) / count(avgVib)) * 100 as dailyVibScore,
    (1.0 - rank(maxTemp, true) / count(maxTemp)) * 100 as dailyTempScore,
    (1.0 - rank(currentPerRpm, true) / count(currentPerRpm)) * 100 as dailyEffScore
from dailyScore
context by day

// 7天滑动平均健康分（窗口覆盖全部数据）
select deviceId, day, dailyVibScore,
    mavg(dailyVibScore, 3) as vibScoreTrend,
    dailyVibScore - move(dailyVibScore, 1) as vibScoreChange
from (
    select deviceId, date(ts) as day,
        avg(vibration) as avgVib,
        (1.0 - rank(avg(vibration), true) / count(avg(vibration))) * 100 as dailyVibScore
    from loadTable("dfs://device_health", "sensor")
    group by deviceId, date(ts)
    context by day
)
context by deviceId csort day

// 定义劣化检测函数
@state
def degradeAlert(vibScore, threshold, declineDays) {
    // 如果当前评分低于阈值且连续下降N天，触发预警
    decline = deltas(vibScore, 1) < 0
    consecutiveDecline = segment(decline, 0)
    return iif(vibScore < threshold and decline, 1, 0)
}

// 策略A：日期VALUE + 设备HASH
db1 = database("", VALUE, 2024.01.01..2024.12.31)
db2 = database("", HASH, [SYMBOL, 10])
db = database("dfs://health_a", COMPO, [db1, db2])

// 策略B：仅日期VALUE
db1 = database("", VALUE, 2024.01.01..2024.12.31)
db = database("dfs://health_b", db1)

// 策略C：日期RANGE + 设备VALUE
db1 = database("", RANGE, 2024.01.01 2024.06.01 2025.01.01)
db2 = database("", VALUE, "DEV_001" "DEV_010" "DEV_020" "DEV_030" "DEV_040" "DEV_050")
db = database("dfs://health_c", COMPO, [db1, db2])

// 好：指定日期范围，只扫描相关分区
select ... from sensor
where date(ts) between 2024.06.01 and 2024.06.07
group by deviceId

// 差：不带时间条件，全表扫描
select ... from sensor
group by deviceId

// 正确：csort 保证时间顺序
select deviceId, day, mavg(score, 3) as trend
from dailyScore
context by deviceId csort day

// 错误：缺少csort，mavg结果不可靠
select deviceId, day, mavg(score, 3) as trend
from dailyScore
context by deviceId

// 定义每日健康评分计算函数
def calcDailyHealth() {
    // 计算昨天的设备特征
    yesterday = date(now()) - 1

    rawMetrics = select deviceId,
        avg(vibration) as avgVib,
        max(temperature) as maxTemp,
        avg(current) / avg(rpm) * 1000 as currentPerRpm,
        std(vibration) / avg(vibration) as cvVib
    from loadTable("dfs://device_health", "sensor")
    where date(ts) = yesterday
    group by deviceId

    // 截面排名并计算评分
    result = select deviceId, yesterday as scoreDate,
        (1.0 - rank(avgVib, true) / count(avgVib)) * 100 as vibScore,
        (1.0 - rank(maxTemp, true) / count(maxTemp)) * 100 as tempScore,
        (1.0 - rank(currentPerRpm, true) / count(currentPerRpm)) * 100 as effScore
    from rawMetrics
    context by NULL

    // 写入评分结果表
    loadTable("dfs://device_health", "health_score").append!(result)
}

// 注册每天凌晨2点执行的定时任务
scheduleJob("dailyHealth", "Calculate daily health score",
    calcDailyHealth, 0 2 * * *)

// 健康评分结果表
scoreDb = database("", VALUE, 2024.01.01..2025.12.31)
scoreSchema = table(1:0,
    `deviceId`scoreDate`vibScore`tempScore`effScore`,
    [SYMBOL, DATE, DOUBLE, DOUBLE, DOUBLE])
scoreDb.createPartitionedTable(scoreSchema, "health_score", "scoreDate")