# -*- coding: utf-8 -*- """ ICM 2023 E — Light Pollution（光污染风险指标体系与干预优化） 真实赛题：COMAP 照明控制任务（ICM）要求 (1) 开发一个广泛适用的指标，识别一个地点的光污染风险水平（LPRI）； (2) 把该指标应用到四类地点（受保护土地 / 乡村社区 / 郊区社区 / 城市社区）并解释结果； (3) 描述三种干预策略，讨论实施行动及其对光污染总体影响的潜在影响； (4) 选取两个地点，用指标确定各自最有效的干预策略，并讨论其对风险水平的影响； (5) 为某一地点及其最有效干预策略制作一页宣传传单。 真源模型（纯标准库，确定性，无随机噪声；SEED 仅作版本锚定）： 1. 每个地点由 6 个归一化指标刻画（均取自公开文献常用的光污染驱动维度）： L 夜间平均辐亮度（VIIRS 类，越高越亮） G 照明增长速率（近年灯光扩张趋势，越高未来风险越大） P 人口暴露度（人类受光照影响的人数规模） B 生态/保护敏感性（生物多样性与保护地权重，越高越脆弱） S 安全照明需求（犯罪/出行安全对光照的依赖，越高越"需要光"） R 治理成熟度（规划/管制能力，越高越能抑制风险，作减项） 2. 两个分量风险（均 clip 到 [0,1]）： 人类风险 HumanRisk = 0.45·L + 0.35·P + 0.20·(L·(1−S)) # 过亮且安全需求低 → 眩光/浪费 生态风险 EcoRisk = 0.40·L + 0.45·B + 0.15·G # 亮度×生态敏感，叠加增长 3. 综合光污染风险指数 LPRI = w_h·HumanRisk + w_e·EcoRisk，默认 w_h=w_e=0.5 （赛题强调"人类与非人类关切并重"，等权最易辩护；权重敏感性另做 tornado）。 4. 风险分级：LPRI<0.33 低 / 0.33–0.60 中 / >0.60 高。 5. 三干预策略（对指标的结构化削减）： I1 定向遮光 + 暖色 LED 替换：L → 0.65·L（降眩光与生态扰动，成本 0.6） I2 分时调光 + 宵禁：L → 0.75·L，G → 0.50·G（成本 0.4） I3 功能分区 + 暗夜廊道：B → 0.60·B，L → 0.85·L（成本 0.9） 对每地点重算 LPRI，取 ΔLPRI 最大者为"最有效"；ΔLPRI/成本 为成本效益。 6. 面板：16 个代表地点（每类 4 个），值经手选以覆盖低/中/高三档， 用于方法演示与可复现分析；真实建模应接入 VIIRS/WorldPop/WDPA 实测。 输出：assets/problems/data/icm2023e.csv（16 行地点风险表 + 表头） 作者说明：指标权重与干预削减系数为方法论演示值，真实参赛应以文献/实测校准； 本模型聚焦"指标构建 → 分类 → 干预优化"的可复现链路。 """ import os import math SEED = 2023 HERE = os.path.dirname(os.path.abspath(__file__)) # --------------------------------------------------------------------------- # 1. 代表地点面板（4 类 × 4 个；指标均为 [0,1] 归一化） # --------------------------------------------------------------------------- # (id, 名称, 类型, L, G, P, B, S, R) LOCATIONS = [ # 受保护土地 (1, "暗夜保护区 A", "protected", 0.06, 0.10, 0.05, 0.92, 0.25, 0.75), (2, "内陆 reserve B", "protected", 0.12, 0.20, 0.10, 0.85, 0.30, 0.65), (3, "城郊 reserve C", "protected", 0.38, 0.30, 0.25, 0.80, 0.35, 0.55), (4, "海岸湿地 D", "protected", 0.20, 0.25, 0.15, 0.78, 0.30, 0.60), # 乡村社区 (5, "山村 E", "rural", 0.25, 0.30, 0.30, 0.55, 0.45, 0.50), (6, "农田 F", "rural", 0.32, 0.35, 0.35, 0.50, 0.45, 0.48), (7, "平原 G", "rural", 0.40, 0.40, 0.42, 0.45, 0.50, 0.45), (8, "河谷 H", "rural", 0.48, 0.45, 0.48, 0.42, 0.52, 0.42), # 郊区社区 (9, "卫星城 I", "suburban", 0.55, 0.45, 0.55, 0.35, 0.55, 0.45), (10, "近郊 J", "suburban", 0.62, 0.50, 0.62, 0.30, 0.60, 0.42), (11, "都会边缘 K", "suburban", 0.70, 0.55, 0.70, 0.25, 0.65, 0.40), (12, "城缘 L", "suburban", 0.66, 0.52, 0.66, 0.28, 0.62, 0.41), # 城市社区 (13, "城市核心 M", "urban", 0.85, 0.60, 0.90, 0.10, 0.80, 0.35), (14, " downtown N", "urban", 0.92, 0.65, 0.95, 0.08, 0.85, 0.32), (15, "港口 O", "urban", 0.78, 0.58, 0.85, 0.12, 0.78, 0.38), (16, "老城 P", "urban", 0.72, 0.55, 0.82, 0.15, 0.75, 0.40), ] TYPES = ["protected", "rural", "suburban", "urban"] TYPE_CN = {"protected": "受保护土地", "rural": "乡村社区", "suburban": "郊区社区", "urban": "城市社区"} # 权重与分级阈值 W_H = 0.50 W_E = 0.50 THR_LO = 0.33 THR_HI = 0.60 # 干预策略：对 (L,G,P,B,S,R) 的削减 + 相对成本 INTERVENTIONS = { "I1": {"name": "定向遮光+暖色LED替换", "cost": 0.6, "f": lambda L,G,P,B,S,R: (0.65*L, G, P, B, S, R)}, "I2": {"name": "分时调光+宵禁", "cost": 0.4, "f": lambda L,G,P,B,S,R: (0.75*L, 0.50*G, P, B, S, R)}, "I3": {"name": "功能分区+暗夜廊道", "cost": 0.9, "f": lambda L,G,P,B,S,R: (0.85*L, G, P, 0.60*B, S, R)}, } # --------------------------------------------------------------------------- # 2. 核心计算 # --------------------------------------------------------------------------- def _clip(v, lo=0.0, hi=1.0): return max(lo, min(hi, v)) def human_risk(L, G, P, B, S, R): return _clip(0.45*L + 0.35*P + 0.20*(L*(1.0 - S))) def eco_risk(L, G, P, B, S, R): return _clip(0.40*L + 0.45*B + 0.15*G) def lpri(L, G, P, B, S, R, wh=W_H, we=W_E): return _clip(wh*human_risk(L, G, P, B, S, R) + we*eco_risk(L, G, P, B, S, R)) def level(v): if v < THR_LO: return "Low" if v <= THR_HI: return "Moderate" return "High" def _minmax(vals): lo, hi = min(vals), max(vals) if hi - lo < 1e-12: return [0.5]*len(vals) return [(v - lo)/(hi - lo) for v in vals] # --------------------------------------------------------------------------- # 3. 主生成 # --------------------------------------------------------------------------- def gen_icm2023e(seed=SEED, out_csv="icm2023e.csv"): # 基础风险 base = [] for (i, name, t, L, G, P, B, S, R) in LOCATIONS: hr = human_risk(L, G, P, B, S, R) er = eco_risk(L, G, P, B, S, R) v = lpri(L, G, P, B, S, R) base.append({ "id": i, "name": name, "type": t, "L": L, "G": G, "P": P, "B": B, "S": S, "R": R, "hr": hr, "er": er, "lpri": v, "level": level(v), }) # 排名（风险降序） order = sorted(range(len(base)), key=lambda k: base[k]["lpri"], reverse=True) for r, k in enumerate(order): base[k]["rank"] = r + 1 # 干预模拟 for loc in base: L, G, P, B, S, R = loc["L"], loc["G"], loc["P"], loc["B"], loc["S"], loc["R"] loc["interv"] = {} best = None for key, iv in INTERVENTIONS.items(): L2, G2, P2, B2, S2, R2 = iv["f"](L, G, P, B, S, R) v2 = lpri(L2, G2, P2, B2, S2, R2) d = loc["lpri"] - v2 ce = d / iv["cost"] loc["interv"][key] = {"lpri2": v2, "delta": d, "ce": ce, "cost": iv["cost"]} if best is None or d > best[1]: best = (key, d) loc["best"] = best[0] loc["best_delta"] = best[1] # 每类汇总 by_type = {} for t in TYPES: rows = [loc for loc in base if loc["type"] == t] by_type[t] = { "n": len(rows), "mean_lpri": sum(r["lpri"] for r in rows)/len(rows), "mean_hr": sum(r["hr"] for r in rows)/len(rows), "mean_er": sum(r["er"] for r in rows)/len(rows), "counts": {"Low": 0, "Moderate": 0, "High": 0}, } for r in rows: by_type[t]["counts"][r["level"]] += 1 # 全样本统计 n = len(base) mean_lpri = sum(r["lpri"] for r in base)/n cnt = {"Low": 0, "Moderate": 0, "High": 0} for r in base: cnt[r["level"]] += 1 # 权重敏感性：w_e 从 0.30 到 0.70，观察分级是否翻转 sens = [] for we in [0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70]: wh = 1.0 - we flips = 0 for loc in base: v2 = lpri(loc["L"], loc["G"], loc["P"], loc["B"], loc["S"], loc["R"], wh, we) if level(v2) != loc["level"]: flips += 1 sens.append({"we": we, "wh": wh, "flips": flips}) # 指标与 LPRI 的相关性（Pearson） def _corr(a, b): m = len(a) ma = sum(a)/m; mb = sum(b)/m cov = sum((a[i]-ma)*(b[i]-mb) for i in range(m)) sa = math.sqrt(sum((x-ma)**2 for x in a)) sb = math.sqrt(sum((x-mb)**2 for x in b)) return cov/(sa*sb) if sa*sb > 0 else 0.0 Ls = [r["L"] for r in base]; Gs=[r["G"] for r in base]; Ps=[r["P"] for r in base] Bs = [r["B"] for r in base]; Ss=[r["S"] for r in base]; Rs=[r["R"] for r in base] V = [r["lpri"] for r in base] corr = { "L": _corr(Ls, V), "G": _corr(Gs, V), "P": _corr(Ps, V), "B": _corr(Bs, V), "S": _corr(Ss, V), "R": _corr(Rs, V), } # 选点：城郊 reserve C(id3) 与 城市核心 M(id13) pick = [3, 13] picked = [loc for loc in base if loc["id"] in pick] D = { "N": n, "base": base, "order": order, "by_type": by_type, "mean_lpri": mean_lpri, "counts": cnt, "sens": sens, "corr": corr, "interventions": INTERVENTIONS, "picked": picked, "W_H": W_H, "W_E": W_E, "THR_LO": THR_LO, "THR_HI": THR_HI, "TYPES": TYPES, "TYPE_CN": TYPE_CN, } # 写 CSV（16 行地点风险表） out = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data", out_csv)) os.makedirs(os.path.dirname(out), exist_ok=True) lines = ["id,name,type,L,G,P,B,S,R,human_risk,eco_risk,lpri,risk_level,rank,best_interv,best_delta"] for r in base: lines.append("%d,%s,%s,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f,%s,%d,%s,%.4f" % (r["id"], r["name"], r["type"], r["L"], r["G"], r["P"], r["B"], r["S"], r["R"], r["hr"], r["er"], r["lpri"], r["level"], r["rank"], r["best"], r["best_delta"])) with open(out, "w", encoding="utf-8-sig") as f: f.write("\n".join(lines) + "\n") D["csv"] = out return D if __name__ == "__main__": d = gen_icm2023e() print("ICM2023E 光污染风险 生成完成 ->", d["csv"]) print("面板 N=%d 平均LPRI=%.4f" % (d["N"], d["mean_lpri"])) print("风险分级:", d["counts"]) print("四类平均LPRI:") for t in d["TYPES"]: bt = d["by_type"][t] print(" %s: %.4f (Low=%d Mod=%d High=%d)" % ( d["TYPE_CN"][t], bt["mean_lpri"], bt["counts"]["Low"], bt["counts"]["Moderate"], bt["counts"]["High"])) print("风险排名 TOP5:") for k in d["order"][:5]: r = d["base"][k] print(" #%d %s [%s] LPRI=%.4f" % (r["rank"], r["name"], r["level"], r["lpri"])) print("指标-LPRI 相关性:", {k: round(v, 4) for k, v in d["corr"].items()}) print("权重敏感性(翻转数):", [(s["we"], s["flips"]) for s in d["sens"]]) for loc in d["picked"]: print("选点 %s (LPRI=%.4f) 最有效干预=%s Δ=%.4f" % ( loc["name"], loc["lpri"], loc["best"], loc["best_delta"])) for key in d["interventions"]: iv = loc["interv"][key] print(" %s: LPRI->%.4f Δ=%.4f CE=%.3f" % ( key, iv["lpri2"], iv["delta"], iv["ce"]))