帮我优化代码

Prompt

**【角色设定】** 你是一位精通GIS空间分析、InSAR时序数据处理（PS-InSAR）以及Python自动化制图的专家级算法工程师。请帮我编写完整的Python实验代码，以复现铁路沉降监测论文的实验部分。不需要复现按高程筛选桥梁 **【项目目录结构】** 项目根目录为 `E:/360MoveData/Users/Lenovo/Desktop/olso`，包含以下两个结构： - olso/psgpt2/：用于存放可执行的Python脚本。 - olso/psgpt2/data:用于存放输出的图片和数据。 - `data/`：用于存放所有输入数据，包含以下文件： - tron_ps_fullpoints.shp：铁路线附近PS点矢量数据 - `tron_ps_fullpoints.csv`：铁路线附近PS点数据，相干性在大于0.6，包含7年76景的时间序列位移属性。用于计算**贪婪启发式分段回归**和**RANSAC（随机抽样一致）算法**。内容示例： | velocity | coherence | MuSigma | shp nbr | Scatterer | Hprecision | Vprecision | range | azimuth | SubArea ID | lon | lat | xpos | ypos | zpos | Z | ALOS | ILOS | Hcorrectio | D_20150512 | ...(共79景) | D_20210901 | D_20210913 | D_20210925 | | -------- | --------- | ------- | ------- | --------- | ---------- | ---------- | ----- | ------- | ---------- | ------- | -------- | ------- | -------- | -------- | -------- | ------- | ------ | ---------- | ---------- | --------- | ---------- | ---------- | ---------- | | 0.007 | 0.612 | 1.922 | 0 | PS | 5.032 | 0.263 | 2078 | 5 | 27-30-33 | 10.4294 | 63.44662 | 10.4294 | 63.44662 | 60.72296 | 20.84884 | -78.469 | 37.054 | 9.765 | 0 | ... | 0.5 | 1.1 | -0.9 | validation_sources.csv：论文中测量车和事故日志报告的异常 | source | start_km | end_km | rate_mm_year | anomaly_type | label | | ----------------- | -------- | ------ | ------------ | ------------ | ------------------- | | Event-Log | 506.3 | 506.3 | | disturbance | ILIAS 506.3 | | Measurement-Train | 508.29 | 508.46 | | unknown | MTIAS 508.29-508.46 | - `Støren501.8kmandTrondheim553.8kmdot.shp`：铁路线起点`Støren501.8km和终点Trondheim553.8km大概位置shp数据 - `Støren501.8km-Trondheim553.8kmrailwayline.shp`：铁路线矢量shp数据。利用 **反距离权重法（IDW）** 插值到铁路线上，将连续的铁路线“离散化为一系列离散点”，计算铁路沉降速率。 【方法复现——逐步实现】 Step 1 预处理 preprocess.py a. PS 点与铁路线统一转 EPSG:32632；coherence ≥ 0.6 复核过滤； b. 999999/NaN 统一置 NaN；有效观测 < MIN_OBS(20) 景的点剔除并记录数量； c. LOS→垂直投影：d_vert = d_LOS / cos(ILOS·π/180)（开关 LOS_TO_VERTICAL）； 与 velocity 列做相关性与量纲校核（自动判别 velocity 单位并打印偏差统计），若整体反号可用 FLIP_SIGN 纠正； d. 输出清洗后的 N×M 时序矩阵与 PS 元数据表（CSV 落盘）。 Step 2 时间域：贪婪启发式分段回归 + RANSAC temporal_processing.py（对应论文 3.1.1 与 Fig.2） 对每个 PS 点的 (t_i, d_i) 序列： 1. 初始分段数 k = ceil(n/2)（n 为有效观测数，即每段约 2 个点，末段可能 3 个）； 2. 对每对相邻段分别做 OLS 直线拟合（各段直线互不强制连续、允许跳跃），记录合并该段对后的 总 RMSE 增量 Δ = RMSE(合并段) − RMSE(段a) − RMSE(段b)； 3. 迭代合并 Δ 最小的相邻段对，直到"最优合并仍使总 RMSE 相对增加 > 3%"（MERGE_RMSE_TOL=0.03） 或仅剩 1 段； 4. 对分段结果施加 RANSAC（sklearn RANSACRegressor + LinearRegression，random_state=42， 残差阈值默认 2×MuSigma，与论文 coherence=0.6 对应的理论噪声色散 4.5 mm 同量级，可配置）， 稳健估计位移速率并剔除离群点； 5. 每点输出：rate_overall（全程 RANSAC 斜率，mm/yr，默认采用）、rate_last_segment（末段斜率，可选）、 分段数、断点位置、内外点掩膜、R²_分段 与 R²_全程线性； 6. 挑 3–4 个代表性 PS 点（分段线性型 / 准线性型 / 含离群点型）绘图：原始散点 + 普通 OLS 直线 + 分段直线 + RANSAC 直线，标注两种拟合的 R² 对比（复现 Fig.2）。 Step 3 空间域：铁路离散化与里程系统 spatial_processing.py a. 铁路线按弧长参数化，以 Støren 端为里程零点：km(s) = 501.8 + (s / S_total) × 52.0， 并校核 S_total ≈ 52 km；打印首末中心点里程与 dots 点位置做一致性检查； b. 以步长 = L（默认 100 m）沿線生成截断窗中心离散点（对应 Fig.11(a) 紫点），记录每点 km； c. 将全部 PS 点垂直投影到线上（shapely project，向量化）：得到每点里程 km_j 与垂距 d_j。 Step 4 IDW 截断窗插值 spatial_processing.py（对应论文 3.1.2 与 Fig.6） 对每个窗口中心 i（沿线矩形：长 L × 宽 W，默认 L=100 m、W=40 m，均可在 config 修改）： - 窗内 PS 点条件：|km_j − km_i| ≤ L/2 且 d_j ≤ W/2； - 插值公式：R_i = Σ( ż_j / D_j^p ) / Σ( 1 / D_j^p )；其中 ż_j = |速率|（绝对值参与插值， 避免沉降/隆起正负抵消），D_j 默认取 PS 点到铁路中心线的垂距（可配置改为到窗口中心的距离）， p 为幂指数（默认 0.5；p=0 退化为算术平均，p→+∞ 退化为最近邻）； - 同步记录每窗 PS 点数 n_p 与密度 n_p/(L·W)； - 绘 Fig.6 风格示意图：截断窗矩形 + 窗内 PS 点按权重 1/D^p 着色 + 沿线纵向权重均匀的说明。 Step 5 异常识别 anomaly_detection.py（对应 3.1.3 与 Fig.11(a)） a. μ+3σ 准则：对所有窗口中心的 |R_i| 计算 μ 与 σ，|R_i| > μ + 3σ 判为异常（正态下约 0.3% 超限）； b. 将相邻异常窗口合并为 ISIAS 段，输出 km 区间（如 "508.29–508.46"）； c. 绘图：横轴 DOVB 里程（501.8→553.8）、纵轴 |R|（mm/yr），紫点=窗口中心，红色虚线=μ+3σ， 灰色虚线=铁路平均速率 μ（论文参考值 4.3 mm/yr），并叠加 ILIAS/MTIAS 标记； d. 各异常段内分别统计强沉降（红）与强隆起（蓝）PS 点数并出图（论文 3.1.2 末段要求）。 Step 6 验证与匹配 validation.py（对应 4.2 节与 Table 3） a. ISIAS 与 ILIAS/MTIAS 匹配：km 区间重叠，或区间最小里程差 ≤ 400 m 判为对应； 另统计 ≤ 300 m 的"定位精度"比例； b. 输出 Table 3 风格对照表 CSV（ILIAS 位置/描述、MTIAS 位置、ISIAS 位置、偏差米数、匹配结论） 与汇总统计：n_ISIAS、与 ILIAS 匹配数、其余与 MTIAS 匹配数、未匹配（疑似误判）段及其窗内 PS 点数 （论文锚点：约 7 个 ISIAS，4 个对应事故日志，余 3 个中 2 个被测量车捕获，1 个疑似误判且位于 PS 稀疏区）； c. 对未匹配段回查其窗内 PS 点数与密度，验证"PS 稀疏 → 误判"机制（论文 5.1 节， 如 536.9 km 处仅 8 个 PS 点）。 Step 7 传统方法对照 traditional_baseline.py（对应 Fig.12） 仅取垂距 ≤ BALLAST_HALF_WIDTH(10 m，道砟带半宽，可配置) 的 PS 点，按最近铁路离散点直接投影取平均 （不做 IDW、不用截断窗），同样施加 μ+3σ；输出与所提方法的对比图与统计（传统法应呈现： 无数据段更多、误判更多），佐证所提方法更优。 Step 8 参数研究与 PS 密度分析 parametric_study.py（对应 5.1–5.2 节与 Fig.17–20） a. 幂指数 p ∈ {0, 0.5, 1, 2, +∞}（固定 L=100, W=40）； b. 窗长 L ∈ {20, 40, 100, 200, 500} m（固定 p=0.5, W=40）； c. 窗宽 W ∈ {10, 20, 40, 100, 200} m（固定 p=0.5, L=100）； 每组记录 n_ISIAS、μ+3σ 阈值、平均 |R|、异常段列表 → 三张对比图 + parametric_summary.csv； 趋势须与论文一致：p↑ → 异常数量与幅值↑（平均效应减弱）；L/W 过小 → PS 稀疏、离群点主导； 过大 → 过度平滑、定位精度下降； d. PS 密度分析（Fig.17）：沿线各窗 PS 点数与密度 n_p/(L·W) 分布图，标注低密度区段。 【关键默认参数（全部集中在 config.py）】 P_POWER=0.5; WIN_L=100; WIN_W=40; WIN_STEP=WIN_L; COHERENCE_MIN=0.6; MERGE_RMSE_TOL=0.03; SIGMA_K=3; MATCH_TOL_M=400; PRECISION_TOL_M=300; KM_START=501.8; KM_END=553.8; RANSAC_SEED=42; LOS_TO_VERTICAL=True; FLIP_SIGN=False; BALLAST_HALF_WIDTH=10; MIN_OBS=20; RATE_MODE="overall" # 可选 "last_segment" 【工程与输出要求】 - 模块划分：config.py、utils_io.py、preprocess.py、temporal_processing.py、spatial_processing.py、 anomaly_detection.py、validation.py、traditional_baseline.py、parametric_study.py、visualization.py、 main.py（一条命令跑通全流程）；允许合理合并文件，但必须保留 main.py 总入口； - 所有中间结果落盘 psgpt2/data/（CSV），支持断点续跑；窗口归属计算用"里程排序 + np.searchsorted" 向量化实现，避免逐窗空间连接；随机过程固定种子； - 输出文件名（全 ASCII，图内文字用英文避免中文乱码，dpi=300）： fig02_piecewise_regression_examples.png、fig06_idw_weight_distribution.png、 fig11a_rate_profile_with_criteria.png、fig11b_anomaly_map.png （铁路线底图 + ISIAS 异常段高亮 + ILIAS/MTIAS 标记）、fig12_traditional_comparison.png、 fig17_ps_density.png、fig18_param_power_index.png、fig19_param_window_length.png、 fig20_param_window_width.png、 ps_points_rate.csv、railway_displacement_profile.csv、ISIAS_segments.csv、 table3_validation_match.csv、parametric_summary.csv、summary_report.md、run.log； - summary_report.md 汇总：清洗后 PS 点数、铁路平均 |R|、μ+3σ 阈值、ISIAS 数与 km 区间、 匹配统计（x/y、300 m 精度比例）、与传统法对比、参数研究结论，并用 [PASS]/[WARN] 软校核论文参考值； - 环境：Python ≥ 3.10；numpy、pandas、geopandas、shapely≥2.0、pyproj、scikit-learn、scipy、matplotlib， 生成 requirements.txt 固化版本；函数带 docstring，关键步骤写日志。 【验收标准】 1. main.py 一键跑通、无手工干预，全部输出落盘； 2. 铁路平均 |R| 处于 2–6 mm/yr 量级（论文 4.3 mm/yr）； 3. ISIAS 约 5–9 段（论文 7 段）；≥4 段在 400 m 内匹配 ILIAS，其余多数能匹配 MTIAS； 4. 若存在未匹配段，其窗内 PS 点数应显著偏低（复现 536.9 km 附近仅 8 点的机制）； 5. 参数研究趋势与论文 5.2 节文字描述一致； 6. 图表要素（坐标轴、单位、参考线颜色语义、图层符号）与论文对应插图一致。