印度尼西亚Morowali工业园区镍加工扩张对沿海水体透明度下降的因果归因 Causal Attribution of Coastal Water Clarity Degradation to Nickel Processing Expansion at the Indonesia Morowali Industrial Park, Sulawesi

方法步骤的完整描述包括四个主要阶段。第一阶段是LULC变化分析：从Esri 10m年度土地利用土地覆盖产品提取2017-2024年影响区的年度地图，排除水体和云掩膜像素，保留6个陆地类别。对每对连续年度地图构建像素级交叉制表矩阵$C(t)=[C^{(t)}_{ij}]$，其中$C^{(t)}_{ij}$表示在时间$Y_t$被分类为类别$i$、在$Y_{t+1}$为类别$j$的像素数。应用Aldwaik和Pontius的三层次强度分析框架，在间隔、类别和转换层次评估变化的系统性特征，并使用QES分解量化变化性质。第二阶段是海洋时间序列探索性分析：从CMEMS GlobColour多传感器合并产品提取1998-2024年每月Kd(490)场，从GLORYS12V1全球海洋再分析提取SST和SSS。计算区域加权的空间均值，使用纬度余弦权重。将记录划分为四个政策对齐的时期：冶炼前基线（2015年4月之前）、初始IMIP运营（2015年4月至2019年12月）、出口禁令后超扩张（2020年1月之后）和完整记录。对每个时期和区域计算鲁棒的位置和尺度度量，使用Theil-Sen估计器和Kendall's $\tau$估计单调趋势。第三阶段是结构断点检测：对每个区域的月度Kd(490)时间序列独立应用多算法共识断点框架。采用径向基函数核成本$C(x_{a+1:b})=\frac{1}{n_{ab}}-\frac{1}{n_{ab}^2}\sum_{i=a+1}^b\sum_{j=a+1}^b\exp(-\gamma(x_i-x_j)^2)$作为段同质性度量。三种算法解决惩罚分割问题：PELT采用线性惩罚$f(k)=\beta k$，Binary Segmentation和基于窗口的检测。通过敏感度扫描选择PELT的惩罚参数$\beta$，通过BIC准则$BIC(k)=n\ln(RSS(k)/n)+pk\ln n$独立选择断点数。共识规则要求至少两个算法在±5个月容差窗口内达成一致。对每个共识断点通过排列检验评估显著性，并使用Welch's t检验、Mann-Whitney U检验、Levene检验和KS检验量化相邻制度间的分布变化。第四阶段是BSTS因果影响估计：将干预日期设为$\hat{\tau}_1$（最早的共识断点），将记录划分为干预前期$T_{pre}=\{1,\dots,t^*-1\}$和干预后期$T_{post}=\{t^*,\dots,n\}$。观测方程为$y_t=\mu_t+\gamma_t+\beta^\top z_t+\varepsilon_t$，其中$\mu_t$是局部线性趋势，$\gamma_t$是周期为$s=12$个月的随机季节分量，$\beta$是静态回归系数向量，$\varepsilon_t \sim \mathcal{N}(0,\sigma^2_\varepsilon)$。趋势演化为$\mu_{t+1}=\mu_t+\nu_t+\eta_{\mu,t}$，$\nu_{t+1}=\nu_t+\eta_{\nu,t}$，其中$\eta_{\mu,t} \sim \mathcal{N}(0,\sigma^2_\mu)$，$\eta_{\nu,t} \sim \mathcal{N}(0,\sigma^2_\nu)$。季节分量满足$\sum_{j=0}^{s-1}\gamma_{t-j}=\eta_{\gamma,t}$，$\eta_{\gamma,t} \sim \mathcal{N}(0,\sigma^2_\gamma)$。通过最大似然估计超参数，使用L-BFGS-B优化器（最多1000次迭代），通过卡尔曼滤波递归评估对数似然。模型在$T_{pre}$上拟合，使用观测的$z_t$在$T_{post}$上向前投影产生反事实。点态和累积因果效应为$\delta_t=y_t-\hat{y}_t$和$\Delta_T=\sum_{t=t^*}^T\delta_t$。通过双侧z检验评估显著性。进行三个稳健性检验：安慰剂测试、留一协变量敏感性分析和前期残差诊断。通过方程$Z_{eu} \approx 4.6/Kd(490)$将因果效应转换为真光层深度变化。

技术新颖性体现在四个方面。第一，方法整合了贝叶斯结构时间序列因果推断、多算法共识断点检测和卫星土地利用强度分析，形成了一个多重证据链的环境影响评估框架。这种组合在沿海环境影响研究中是新颖的，能够穿透由季风、ENSO和印度洋偶极子驱动的巨大自然变异性。第二，方法使用长达27年的多传感器合并海洋颜色数据（SeaWiFS、MODIS-Aqua、MERIS、VIIRS-SNPP、VIIRS-JPSS1和Sentinel-3A/3B上的OLCI），构建了深度干预前记录以实现稳定的反事实估计。第三，方法采用了分布无关的安慰剂秩检验，不依赖于高斯假设，为因果效应提供了独立于模型假设的显著性验证。第四，方法在4km的GlobColour分辨率下操作，整合了远大于单个羽流区域的光学性质，使得估计是保守下界，同时保持了因果推断所需的深度时间记录。

独立分析的弱点包括：第一，空间分辨率限制。4km的GlobColour产品整合了远大于单个羽流区域的光学性质，可能低估峰值浊度和羽流传播的空间细节。改进方向：融合更高分辨率传感器（如Sentinel-3 OLCI的300m数据）进行近期时段的亚像素梯度分析，虽然这需要处理不同传感器间的光谱和几何差异。第二，水质代理的特异性限制。Kd(490)不能区分悬浮沉积、有色溶解有机物质和浮游植物色素的具体贡献。改进方向：结合额外的水质参数如总悬浮物(TSM)、有色溶解有机物质(CDOM)和叶绿素a浓度，使用多变量光学模型分解Kd(490)的组成贡献。第三，缺乏原位验证。没有同步的原位测量排除了直接地面验证。改进方向：与当地机构合作部署锚系观测站，测量近岸浊度、沉积物通量和痕量金属浓度，建立卫星反演的地面真值。第四，控制区识别假设。方法依赖于控制区不受干预影响的假设，虽然零结果支持这一点。改进方向：扩展控制区到多个独立地点，进行敏感性分析，评估控制区选择对结果的影响。第五，静态回归系数假设。BSTS模型假设静态回归系数，可能无法捕获干预后时期的潜在非线性动态。改进方向：使用时变系数BSTS模型或机器学习方法（如高斯过程）捕获非平稳关系。

未来研究方向包括作者提出的和基于成果可延伸的。作者提出的方向：复制IMIP模式于Weda Bay（北马鲁古）和Virtue Dragon综合体（东南苏拉威西），研究这些影响是否在类似设施位于贫营养水域的地方被复制。我们的延伸方向：第一，扩展框架到其他沿海工业类型。将BSTS因果推断框架应用于沿海采矿、港口建设、水产养殖和石油天然气平台等其他沿海工业，评估它们的海洋环境影响。第二，集成更高分辨率数据源。融合PlanetScope（3m）、Sentinel-2 MSI（10m）和Landsat 8/9（30m）的数据，实现亚羽流尺度的空间分析，同时保持深度时间记录。第三，多变量因果分析。扩展BSTS框架到多变量设置，同时建模多个水质参数（Kd(490)、TSM、Chl-a、CDOM）及其相互依赖关系，提供更全面的环境影响评估。第四，生态系统建模集成。将BSTS估计的因果效应集成到生态系统模型中，预测对珊瑚礁、海草床和鱼类种群的长期生态影响。第五，政策影响评估。评估印尼镍下游化政策的环境成本，为环境监管和可持续发展提供科学依据。第六，机器学习方法探索。使用深度学习、图神经网络和因果机器学习方法，捕获更复杂的非线性关系和空间依赖性。第七，实时监测系统。开发基于卫星的实时监测和预警系统，为环境管理和应急响应提供及时信息。

复现评估：开源情况良好，所有原始卫星和海洋数据集可从其各自公共存储库获得：海洋颜色和物理海洋变量通过CMEMS（https://marine.copernicus.eu），LULC复合物通过Esri（https://livingatlas.arcgis.com/landcoverexplorer），水深通过SRTM15+V2（https://topex.ucsd.edu/WWW_html/srtm15_plus.html）。处理后的数据集、脚本和补充材料可在https://github.com/sandyherho/supplMorowaliOcean获得。代码使用Python实现，依赖NumPy、SciPy、netCDF4、Matplotlib、xarray、pandas和statsmodels库，都是标准开源工具。数据量适中：324个每月Kd(490)场，8个年度LULC地图，SST和SSS时间序列。算力需求不高：统计分析可在标准笔记本电脑或台式机上运行，断点检测和BSTS建模在几分钟到几小时内完成。难度中等到高：需要理解时间序列分析、因果推断和海洋光学原理，以及Python编程和数据处理技能。方法的可转移性强：框架可应用于其他沿海工业场所和热带环境，数据可获取性和代码开源性支持复现和扩展。