在生态环境、自然地理及可持续发展研究中，面对海量的散点数据，传统的线性回归模型通常难以全面反映变量间的真实关系。由于真实生态系统的高度复杂性，生态系统服务与气象、地形、社会经济等驱动因子之间，往往存在显著的非线性交互作用（Nonlinear Interactions）和阈值效应（Threshold Effects）。
为了准确刻画这些非线性特征，“双约束线”（Dual-Constraint Lines）分析框架被提出并广泛应用。本文将系统梳理该方法的底层统计原理与典型应用场景，并提供完整的 Python 实现代码，以期为相关领域的数据分析与可视化工作提供客观参考。
替换自变量与因变量即可使用！
在散点图中，数据点往往呈现出一团不规则的“云”状分布。传统的回归分析（如 OLS）关注的是数据的“平均趋势”，直接穿透点云的中心。但这忽略了一个致命问题：在生态系统中，起决定作用的往往不是“平均水平”，而是“极限边界”。
“约束线”方法（Constraint Line Method）的核心思想正是去寻找数据分布的边界包络线。而“双约束线”框架，顾名思义，同时关注了上限和下限：
双约束线的提取并不是简单地把最外围的点连起来（那样极其容易受到极端异常值的干扰），而是基于严谨的统计学步骤：
分段分位数回归（Segmented Quantile Regression）： 将 X 轴划分为多个等距区间（Bins），在每个区间内提取高分位点（如 99.9%）作为上限候选点，提取低分位点（如 0.1%）作为下限候选点。
Tukey 栅栏法去噪： 利用四分位距（IQR）对提取出的边界点进行过滤，剔除偶发的极端异常值。
受限三次样条拟合（Restricted Cubic Splines, RCS）： 这是最关键的一步。普通的平滑曲线在数据两端容易出现剧烈震荡（过拟合），而 RCS 强制规定了曲线在两端节点外呈线性，从而在保证曲线整体平滑且连续的同时，极大地提高了边界拟合的稳健性。
掌握了双约束线，你能用它在论文里写出什么样的好故事？以下是三大经典应用：
1. 精准识别“生态阈值”（Threshold Identification）上约束线通常呈现“倒 U 型”或“单调饱和型”。倒 U 型曲线的顶点（导数为 0 的位置）即为关键阈值点（Tipping Point）。例如，随着植被覆盖度（NDVI）增加，产水量（Water Yield）最初会协同上升；但当 NDVI 超过某个阈值（如 0.57）后，过度的蒸腾作用反而会导致产水量下降。找到这个阈值，就能为“退耕还林”等生态工程提供科学的造林密度指导。
2. 揭示隐藏的“生态风险”（Ecological Risk Warning）许多研究只看上限，却忽略了下限。双约束线的巧妙之处在于，将拟合出的下约束线与政策规定的“生态红线”（Baseline）进行比对。如果某个区域的下约束线跌破了基线（例如最低生态基流、最低碳汇标准），即使它的上限表现很好，这里也是一个亟待干预的高风险脆弱区。
3. 洞察人类活动对自然机制的“扭曲”（Anthropogenic Disturbance）通过对比自然保护区（无干扰）与人类高强度活动区（如农田、城镇）的双约束线，可以清晰地看到人类活动如何导致阈值延迟（Threshold Postponement）或上限衰减（Attenuation）。这为协调经济发展与生态安全提供了极具说服力的直观证据。
下面这段 Python 代码为您打通了从数据加载、区间切片、异常值剔除、RCS平滑拟合，到极值点（阈值）计算与最终出图的全流程。代码严格遵守顶刊规范：设置了 Times New Roman（英文）与 宋体（中文）的字体组合，带有置信区间阴影，并自动标注出阈值坐标与拟合优度 R2。
请确保您的运行环境中已安装：pandas, numpy, matplotlib, statsmodels, patsy。
import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport statsmodels.formula.api as smffrom patsy import crimport warningswarnings.filterwarnings('ignore')# 地数空间-地学小助手-太阳花🌻# ==========================================# 1. 字体与绘图全局设置# ==========================================# 设置全局字体：英文为 Times New Roman，中文为宋体plt.rcParams['font.sans-serif'] = ['SimSun'] # 用来正常显示中文标签plt.rcParams['font.family'] = 'sans-serif'plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号plt.rcParams['font.serif'] = ['Times New Roman']# ==========================================# 2. 数据加载与参数设置# ==========================================# 请确认文件路径正确file_path = r"E:\data_study1双约束线\data1.xlsx"df = pd.read_excel(file_path)# 定义变量target_y = 'HQ' # 因变量features_x = ['DEM', 'dist_waterway', 'TEMP', 'PRCP', 'Population', 'LHGI', 'GPP', 'FTI', 'LE', 'ETP', 'Nightlight', 'RX5day', 'CDD', 'TNn', 'TXx']# ==========================================# 3. 核心算法函数 (新增统计指标提取)# ==========================================def tukey_fences_filter(x, y): """使用 Tukey's fences (IQR) 剔除异常值""" q1 = np.percentile(y, 25) q3 = np.percentile(y, 75) iqr = q3 - q1 lower_bound = q1 - 1.5 * iqr upper_bound = q3 + 1.5 * iqr mask = (y >= lower_bound) & (y <= upper_bound) return x[mask], y[mask]def fit_restricted_cubic_spline(x, y, df_spline=4): """对过滤后的散点进行受限三次样条拟合(RCS)，并返回完整统计指标""" data = pd.DataFrame({'x': x, 'y': y}) data = data.dropna().sort_values('x') # 使用 patsy 创建受限三次样条，df=4 通常设置3个内部节点 formula = "y ~ cr(x, df={})".format(df_spline) try: model = smf.ols(formula, data=data).fit() # 生成平滑的预测线 x_pred = np.linspace(data['x'].min(), data['x'].max(), 300) pred_data = pd.DataFrame({'x': x_pred}) predictions = model.get_prediction(pred_data) y_pred = predictions.predicted_mean y_ci_lower = predictions.conf_int()[:, 0] y_ci_upper = predictions.conf_int()[:, 1] # 计算极值点 (导数为0处) - 即阈值点 dy = np.gradient(y_pred, x_pred) # 寻找平稳点(斜率变号) sign_dy = np.sign(dy) sign_changes = ((np.roll(sign_dy, 1) - sign_dy) != 0).astype(int) sign_changes[0] = 0 peak_idx = np.argmax(y_pred) # 上约束通常找最高点 valley_idx = np.argmin(y_pred) # 下约束找最低点 # ★ 新增：提取丰富的统计指标 ★ stats = { 'R_squared': model.rsquared, 'Adj_R_squared': model.rsquared_adj, 'F_statistic': model.fvalue, 'P_value': model.f_pvalue, 'AIC': model.aic, 'BIC': model.bic, 'RMSE': np.sqrt(model.mse_resid) # 残差均方差的平方根 } return x_pred, y_pred, y_ci_lower, y_ci_upper, x_pred[peak_idx], y_pred[peak_idx], x_pred[valley_idx], y_pred[valley_idx], stats except Exception as e: return None, None, None, None, None, None, None, None, None# ==========================================# 4. 绘图与统计结果收集循环# ==========================================fig = plt.figure(figsize=(20, 25))fig.subplots_adjust(hspace=0.4, wspace=0.3)# 创建一个空列表，用于收集所有的统计信息stats_results_list = []for idx, x_var in enumerate(features_x): ax = fig.add_subplot(5, 3, idx + 1) # 提取当前数据并剔除缺失值 temp_df = df[[x_var, target_y]].dropna() x_data = temp_df[x_var].values y_data = temp_df[target_y].values # 原始数据散点（浅灰色） ax.scatter(x_data, y_data, color='lightgray', alpha=0.5, s=10, label='Original Data') # 1. 将 X 分为 100 个等宽的区间（Bins） bins = np.linspace(np.min(x_data), np.max(x_data), 101) bin_indices = np.digitize(x_data, bins) x_binned = [] y_upper_999 = [] y_lower_001 = [] for b in range(1, 101): in_bin = y_data[bin_indices == b] in_bin_x = x_data[bin_indices == b] if len(in_bin) > 0: x_binned.append(np.mean(in_bin_x)) # 该区间X的均值 y_upper_999.append(np.percentile(in_bin, 99.9)) y_lower_001.append(np.percentile(in_bin, 0.1)) x_binned = np.array(x_binned) y_upper_999 = np.array(y_upper_999) y_lower_001 = np.array(y_lower_001) # 2. Tukey's Fences 异常值剔除 x_up_clean, y_up_clean = tukey_fences_filter(x_binned, y_upper_999) x_low_clean, y_low_clean = tukey_fences_filter(x_binned, y_lower_001) # 绘制上约束过滤后的散点 ax.scatter(x_up_clean, y_up_clean, color='crimson', alpha=0.8, s=20, edgecolors='k', linewidth=0.5) # 绘制下约束过滤后的散点 ax.scatter(x_low_clean, y_low_clean, color='dodgerblue', alpha=0.8, s=20, edgecolors='k', linewidth=0.5) # 3. 受限三次样条拟合 (接收新的统计字典参数) x_pred_up, y_pred_up, ci_low_up, ci_hi_up, peak_x_up, peak_y_up, _, _, stats_up = fit_restricted_cubic_spline(x_up_clean, y_up_clean) x_pred_low, y_pred_low, ci_low_low, ci_hi_low, _, _, valley_x_low, valley_y_low, stats_low = fit_restricted_cubic_spline(x_low_clean, y_low_clean) # 绘制上约束拟合线与置信区间 if x_pred_up is not None: r2_up = stats_up['R_squared'] # 获取R²用于图表展示 ax.plot(x_pred_up, y_pred_up, color='crimson', linewidth=2.5, label='Upper Constraint') ax.fill_between(x_pred_up, ci_low_up, ci_hi_up, color='crimson', alpha=0.15) # 标注阈值线 (峰值) ax.axvline(peak_x_up, color='black', linestyle='--', alpha=0.7) ax.scatter(peak_x_up, peak_y_up, color='white', edgecolors='crimson', s=60, zorder=5) ax.text(peak_x_up, peak_y_up * 1.02, f"({peak_x_up:.2f}, {peak_y_up:.2f})", fontname='Times New Roman', fontsize=10, ha='center', color='crimson') # 标出 R^2 ax.text(0.95, 0.95, f'$R^2$={r2_up:.2f}', transform=ax.transAxes, ha='right', va='top', color='crimson', fontname='Times New Roman', fontsize=10, bbox=dict(facecolor='white', alpha=0.6, edgecolor='crimson')) # 记录上约束统计指标 up_info = {'Variable': x_var, 'Constraint_Type': 'Upper', 'Threshold_X': peak_x_up, 'Threshold_Y': peak_y_up} up_info.update(stats_up) stats_results_list.append(up_info) # 绘制下约束拟合线与置信区间 if x_pred_low is not None: r2_low = stats_low['R_squared'] # 获取R²用于图表展示 ax.plot(x_pred_low, y_pred_low, color='dodgerblue', linewidth=2.5, label='Lower Constraint') ax.fill_between(x_pred_low, ci_low_low, ci_hi_low, color='dodgerblue', alpha=0.15) # 标出 R^2 ax.text(0.95, 0.05, f'$R^2$={r2_low:.2f}', transform=ax.transAxes, ha='right', va='bottom', color='dodgerblue', fontname='Times New Roman', fontsize=10, bbox=dict(facecolor='white', alpha=0.6, edgecolor='dodgerblue')) # 记录下约束统计指标 low_info = {'Variable': x_var, 'Constraint_Type': 'Lower', 'Threshold_X': valley_x_low, 'Threshold_Y': valley_y_low} low_info.update(stats_low) stats_results_list.append(low_info) # 设置网格与外观 ax.grid(True, linestyle='--', alpha=0.4) ax.set_title(f'HQ = f({x_var})', fontname='Times New Roman', fontsize=14, fontweight='bold') ax.set_xlabel(f'{x_var}', fontname='Times New Roman', fontsize=12) ax.set_ylabel('HQ', fontname='Times New Roman', fontsize=12) # 设置刻度字体 for label in ax.get_xticklabels() + ax.get_yticklabels(): label.set_fontname('Times New Roman') label.set_fontsize(11)# ==========================================# 5. 全局排版、导出图片及统计表格# ==========================================# 添加全局图例handles, labels = ax.get_legend_handles_labels()fig.legend(handles, labels, loc='lower center', ncol=3, fontsize=14, prop={'family': 'Times New Roman'}, bbox_to_anchor=(0.5, 0.08))plt.suptitle("Dual-Constraint Framework for Ecosystem Service (HQ) Optimization", fontsize=20, fontname='Times New Roman', y=0.92, fontweight='bold')# ★ 导出图片：同时保存为 PNG(预览用) 和 PDF(发文章用矢量图) ★output_png = r"E:\data_study1双约束线\Dual_Constraint_Plot.png"output_pdf = r"E:\data_study1双约束线\Dual_Constraint_Plot.pdf"plt.savefig(output_png, dpi=300, bbox_inches='tight')plt.savefig(output_pdf, format='pdf', dpi=300, bbox_inches='tight')plt.show()# ★ 导出统计指标至 Excel ★output_excel = r"E:\data_study1双约束线\Dual_Constraint_Statistics.xlsx"stats_df = pd.DataFrame(stats_results_list)stats_df.to_excel(output_excel, index=False, float_format="%.4f")
