简介这份资源面向古生物学、地球科学与气候学领域的科研人员、学生及数据分析爱好者围绕有孔虫化石形态测量数据的气候重建与群落结构评估展开。内容基于377个地质时间点、每点约2000个个体的测量数据涵盖大小、面积、形状因子、伸长率、球形度、周长与灰度值等参数提供完整的Python代码及解释涉及数据加载与预处理、基本统计、时间序列分解、多变量分析、PCA降维及环境关联分析并扩展了年龄模型验证、自动化参数提取与分类单元变化分析。资源包为1个PDF文件约516KB便于携带与查阅。已有59人学习。读者可借此掌握Pandas、NumPy、Matplotlib、Seaborn、statsmodels与scikit-learn在科研场景中的实际用法理解形态变化与环境变化的关联并获得可复用的分析框架与排错思路。1. 有孔虫化石形态测量与气候重建一份能跑通的 Python 复现包如果你手头有一批有孔虫化石的形态测量数据——面积、周长、形状因子、伸长率、球形度、灰度值——却卡在“怎么把它们和古气候指标串起来”这一步这份复现包值得拆开看。它对应的是论文Morphometric Analysis of Foraminifera Fossils for Climate Reconstruction and Community Structure Assessment核心思路是用比“大小分布偏度”更细的形态参数去追踪群落结构、物种数量、选择压力差异、数据周期性以及形态变化与环境变化之间的关联。数据集覆盖 377 个地质时间点每个点约 2000 个个体的测量记录量级不算小。适合古生物学、地球科学、气候学方向的研究生和科研人员也适合想拿真实科学数据练 pandas、numpy、scipy、seaborn 的人。下面按“数据怎么进 → 统计怎么做 → 时间序列怎么拆 → 多变量怎么降维 → 环境怎么关联 → 坑在哪”的顺序走一遍。2. 数据加载与预处理把 377 个 Excel 拼成一张可分析的表2.1 为什么不能直接pd.read_excel一把梭原始数据是“一个样本一个 Excel 文件”的散装结构外加一份样本概览表。直接读单个文件只能看到某个时间点的个体测量值拿不到年龄、环境指标这些跨样本信息。常见做法是先把所有样本文件纵向拼成一张长表再按Sample_ID左连接概览表把Age (Ma)、d18O这类字段挂上去。这里有两个硬约束一是样本文件名必须能解析出Sample_ID二是概览表里的Sample列要和样本 ID 对得上否则 merge 之后全是 NaN后面所有分析都会静默跑出一堆空图。import pandas as pd import numpy as np import os def load_and_preprocess_data(sample_file_path, overview_file_path): 加载并预处理有孔虫形态测量数据 sample_file_path: 样本数据文件夹路径 overview_file_path: 样本概览文件路径 overview pd.read_excel(overview_file_path) combined_data pd.DataFrame() for file in os.listdir(sample_file_path): if file.endswith(.xlsx): sample_id file.split(.)[0] # 文件名即样本 ID sample_data pd.read_excel(os.path.join(sample_file_path, file)) sample_data[Sample_ID] sample_id combined_data pd.concat([combined_data, sample_data], ignore_indexTrue) # 左连接概览表保留所有个体测量记录 combined_data combined_data.merge( overview, left_onSample_ID, right_onSample, howleft ) # 清洗去掉没有年龄的样本去掉面积为 0 或负的无效测量 combined_data combined_data.dropna(subset[Age (Ma)]) combined_data combined_data[combined_data[Area] 0] return combined_data逻辑说明os.listdir遍历文件夹file.split(.)[0]取文件名前缀作为样本 ID这一步假设文件名就是纯 ID没有多余后缀。pd.concat用ignore_indexTrue重置索引避免后续 groupby 出现重复索引。merge 用howleft而不是 inner是为了保留那些在概览表里找不到匹配的样本方便后面排查数据对齐问题。清洗阶段dropna(subset[Age (Ma)])是必须的因为时间序列分析依赖年龄轴Area 0过滤掉分割失败或测量异常产生的零面积对象。参数上sample_file_path指向存放所有.xlsx样本文件的目录overview_file_path是概览表路径。如果概览表里年龄列名不是Age (Ma)要么改代码要么在读入后统一 rename。我一般会在 merge 之后加一句print(combined_data[Age (Ma)].isna().sum())确认没有大量样本因为 ID 对不上而丢失年龄。2.2 列名标准化与形态参数提取不同批次的测量软件导出的列名可能大小写不一比如area、Area、AREA混着来。在进入统计分析前最好做一次列名映射把object_id、sample_id、shape_factor这类变体统一成标准名。同时把形态参数单独拎出来做标准化后面 PCA 和相关性分析都要用。from sklearn.preprocessing import StandardScaler def automated_parameter_extraction(data): 统一列名并标准化形态参数 param_order [ Object_ID, Sample_ID, Area, Perimeter, Shape_Factor, Elongation, Sphericity, Gray_Value, Age (Ma) ] processed_data data.copy() column_map { object_id: Object_ID, sample_id: Sample_ID, area: Area, perimeter: Perimeter, shape_factor: Shape_Factor, elongation: Elongation, sphericity: Sphericity, gray_value: Gray_Value, age: Age (Ma) } processed_data.rename(columnscolumn_map, inplaceTrue) missing_cols [c for c in param_order if c not in processed_data.columns] if missing_cols: print(f警告: 缺少以下列: {missing_cols}) available_cols [c for c in param_order if c in processed_data.columns] processed_data processed_data[available_cols] morph_params [Area, Perimeter, Shape_Factor, Elongation, Sphericity, Gray_Value] morph_params [p for p in morph_params if p in processed_data.columns] scaler StandardScaler() processed_data[morph_params] scaler.fit_transform(processed_data[morph_params]) return processed_data逻辑说明rename用字典做批量列名替换只替换存在的键不会因为某个变体不存在而报错。param_order定义了标准列顺序最后用列表推导筛出实际存在的列保证输出表结构一致。StandardScaler对形态参数做 z-score 标准化均值 0、标准差 1这一步对 PCA 是必须的因为面积和灰度值的量纲差了几个数量级不标准化的话 PCA 会被面积主导。参数上morph_params列表可以根据实际数据增减比如有些数据集没有Sphericity代码会自动跳过。标准化后的值会覆盖原列如果还想保留原始值可以先 copy 一份到_raw后缀列。3. 基本统计与时间序列从偏度到周期性的完整链路3.1 按样本聚合形态统计量拿到长表之后第一步是按Sample_ID分组算每个样本的均值、中位数、标准差和偏度。偏度是论文特别强调的——传统研究只看大小分布的偏度这里把偏度扩展到面积、周长、形状因子等多个参数看不同参数的选择压力是否同步。聚合之后把多级列索引拍平方便后续按列名取用。def basic_statistical_analysis(data): grouped data.groupby(Sample_ID) stats_df grouped.agg({ Area: [mean, median, std, skew], Perimeter: [mean, median, std, skew], Shape_Factor: [mean, median, std, skew], Age (Ma): first }) stats_df.columns [_.join(col).strip() for col in stats_df.columns.values] stats_df.reset_index(inplaceTrue) # 面积均值随时间变化颜色映射偏度 import seaborn as sns import matplotlib.pyplot as plt plt.figure(figsize(12, 6)) sns.scatterplot(datastats_df, xAge (Ma)_first, yArea_mean, hueArea_skew, paletteviridis) plt.title(有孔虫平均面积随时间变化(颜色表示偏度)) plt.xlabel(年龄(百万年)) plt.ylabel(平均面积) plt.colorbar(label偏度) plt.show() return stats_df逻辑说明groupby(Sample_ID)后agg对每个形态参数同时算四个统计量Age (Ma)取first因为同一样本内年龄是常量。列名拍平用_.join(col)把(Area, mean)变成Area_mean避免多级索引在后续操作中带来的麻烦。散点图用hueArea_skew把偏度映射到颜色一眼能看出哪些时间点的分布不对称。参数上skew是 scipy 的偏度计算对样本量敏感如果某个样本个体数太少比如少于 30偏度会很不稳定建议在聚合前先过滤掉个体数过少的样本。paletteviridis是连续色板适合偏度这种连续变量。3.2 时间序列分解与自相关按年龄分组后把Area_mean当作时间序列做分解看趋势、周期和残差。seasonal_decompose的period参数需要根据地质时间尺度来定代码里默认 10意思是每 10 个百万年一个周期这个值不能照搬得结合具体数据的时间分辨率和地质背景调整。自相关图用来判断形态变化是否存在周期性记忆。from statsmodels.tsa.seasonal import seasonal_decompose def time_series_analysis(data): time_stats data.groupby(Age (Ma)).agg({ Area: [mean, std, skew, lambda x: np.percentile(x, 95)], Shape_Factor: [mean, std], Elongation: [mean, std] }) time_stats.columns [Area_mean, Area_std, Area_skew, Area_95th, Shape_Factor_mean, Shape_Factor_std, Elongation_mean, Elongation_std] time_stats time_stats.sort_index() result seasonal_decompose(time_stats[Area_mean], modeladditive, period10) fig result.plot() fig.set_size_inches(12, 8) plt.suptitle(有孔虫面积时间序列分解) plt.show() plt.figure(figsize(10, 5)) pd.plotting.autocorrelation_plot(time_stats[Area_mean]) plt.title(有孔虫平均面积自相关图) plt.show() return time_stats逻辑说明groupby(Age (Ma))按年龄聚合sort_index()保证时间轴单调递增否则分解会乱序。seasonal_decompose用加法模型假设趋势、周期、残差是相加关系如果数据波动幅度随均值增大而增大应该换modelmultiplicative。lambda x: np.percentile(x, 95)算 95 分位数用来捕捉大个体的尾部变化比只看均值更敏感。参数上period10是最大的坑之一。如果年龄分辨率是 0.1 Ma那 10 个百万年就是 100 个点一个周期如果分辨率是 1 Ma就是 10 个点。设错了分解出来的周期成分没有地质意义。自相关图里横轴是滞后阶数纵轴是相关系数超出置信区间的滞后阶数说明存在显著自相关。4. 多变量分析与环境关联PCA 降维和 δ18O 相关4.1 PCA 降维与载荷解读六个形态参数之间往往高度相关直接做回归会共线性爆炸。PCA 把六个参数压成两个主成分看样本在 PC1-PC2 平面上的分布颜色映射年龄能直观判断形态是否随地质时间发生系统性漂移。载荷矩阵告诉你每个主成分主要代表哪些参数。from sklearn.decomposition import PCA def multivariate_analysis(data): features [Area, Perimeter, Shape_Factor, Elongation, Sphericity, Gray_Value] scaler StandardScaler() scaled_data scaler.fit_transform(data[features]) pca PCA(n_components2) principal_components pca.fit_transform(scaled_data) plt.figure(figsize(10, 8)) sns.scatterplot(xprincipal_components[:, 0], yprincipal_components[:, 1], huedata[Age (Ma)], paletteviridis, alpha0.6) plt.title(有孔虫形态数据PCA分析(颜色表示年龄)) plt.xlabel(主成分1 (解释方差: {:.2f}%).format( pca.explained_variance_ratio_[0] * 100)) plt.ylabel(主成分2 (解释方差: {:.2f}%).format( pca.explained_variance_ratio_[1] * 100)) plt.colorbar(label年龄(百万年)) plt.show() loadings pd.DataFrame(pca.components_.T, columns[PC1, PC2], indexfeatures) print(\n主成分载荷矩阵:) print(loadings) return pca, loadings逻辑说明StandardScaler在 PCA 前必须做否则量纲大的参数会主导方差。PCA(n_components2)只保留前两个主成分方便可视化如果解释方差太低比如 PC1PC2 不到 60%说明形态变异不是低维结构需要考虑更多主成分或换方法。载荷矩阵pca.components_.T的每一行对应一个原始参数值越大说明该参数对主成分贡献越大。参数上alpha0.6让散点半透明避免密集区域糊成一团。explained_variance_ratio_给出每个主成分解释的方差比例PC1 通常代表“整体大小”PC2 可能代表“形状扁平度”或“伸长程度”具体看载荷符号。4.2 形态参数与 δ18O 的关联δ18O 是经典的古温度代用指标把它和形态参数放一起看能回答“形态变化是否跟环境变化同步”这个核心问题。代码里先按样本聚合形态均值和 δ18O再画散点图必要时加回归线和相关系数。def environmental_correlation_analysis(data, overview): if d18O not in overview.columns: print(概览表中没有 d18O 列跳过环境关联分析) return merged data.merge(overview[[Sample, d18O]], left_onSample_ID, right_onSample, howleft) env_stats merged.groupby(Sample_ID).agg({ Area: mean, Shape_Factor: mean, d18O: first, Age (Ma): first }) fig, axes plt.subplots(1, 2, figsize(15, 6)) sns.scatterplot(dataenv_stats, xd18O, yArea_mean, axaxes[0]) axes[0].set_title(有孔虫平均面积与δ18O的关系) axes[0].set_xlabel(δ18O (‰)) axes[0].set_ylabel(平均面积) sns.scatterplot(dataenv_stats, xd18O, yShape_Factor_mean, axaxes[1]) axes[1].set_title(有孔虫平均形状因子与δ18O的关系) axes[1].set_xlabel(δ18O (‰)) axes[1].set_ylabel(平均形状因子) plt.tight_layout() plt.show() # 补充计算 Pearson 相关系数和 p 值 from scipy import stats r_area, p_area stats.pearsonr(env_stats[d18O], env_stats[Area_mean]) r_shape, p_shape stats.pearsonr(env_stats[d18O], env_stats[Shape_Factor_mean]) print(f面积 vs δ18O: r{r_area:.3f}, p{p_area:.4f}) print(f形状因子 vs δ18O: r{r_shape:.3f}, p{p_shape:.4f}) return env_stats逻辑说明merge把 δ18O 挂到个体测量表上再groupby(Sample_ID)聚合保证每个样本只有一个 δ18O 值取first。散点图看线性趋势pearsonr给出相关系数和显著性。如果 p 值大于 0.05说明线性相关不显著可能需要考虑非线性关系或滞后效应。参数上d18O列名要和概览表一致有些数据集用delta18O或d18O_permil需要提前 rename。pearsonr假设线性关系和正态分布如果散点图明显非线性应该换 Spearman 相关或加多项式拟合。5. 避坑与排查数据对齐、周期参数和内存这三处最容易翻车5.1 样本 ID 对不上导致年龄全空现象merge 之后Age (Ma)列大量 NaN时间序列图空白。原因样本文件名解析出的Sample_ID和概览表Sample列的格式不一致比如一个是S001另一个是Sample_001或者有前后空格。解决merge 前先overview[Sample] overview[Sample].astype(str).str.strip()样本 ID 也做同样处理再打印两边 ID 集合的交集大小确认匹配率。如果匹配率低于 90%先别往下跑回去查数据源。5.2seasonal_decompose的 period 设错现象分解出的周期成分是一条直线或剧烈震荡没有地质意义。原因period参数和实际时间分辨率不匹配。比如年龄间隔是 0.5 Maperiod10意味着 20 个点一个周期但如果数据只有 50 个时间点周期成分根本拟合不出来。解决先算time_stats.index.to_series().diff().median()得到中位时间间隔再用period 目标周期长度 / 中位间隔反推。目标周期长度要参考地质背景比如轨道周期约 100 ka 或 400 ka。5.3 全量数据做 PCA 内存爆掉现象StandardScaler.fit_transform或PCA.fit_transform报 MemoryError。原因377 个样本 × 2000 个体 ≈ 75 万行六个 float64 列约 36 MB本身不大但如果中间做了多次 copy 或者 merge 产生了笛卡尔积行数会膨胀。解决在 merge 前检查概览表是否有重复Sample值重复会导致一对多合并用overview.drop_duplicates(subsetSample)去重。另外PCA 之前可以按样本抽样比如每个样本随机抽 200 个个体既降内存又降个体测量噪声。5.4 δ18O 列名不匹配静默跳过现象环境关联分析没有任何输出也没有报错。原因代码里用if d18O not in overview.columns做了静默跳过如果列名是d18O_permil或delta_O18直接 return 了。解决把列名检查改成模糊匹配比如[c for c in overview.columns if 18O in c or d18O in c]找到后 rename 成标准名。或者在函数开头打印overview.columns.tolist()确认实际列名。5.5 中文显示方块现象图表标题和轴标签中文变成方块。原因matplotlib 默认字体不支持中文。解决在脚本开头加plt.rcParams[font.sans-serif] [SimHei]和plt.rcParams[axes.unicode_minus] False。如果系统没有 SimHei换成Microsoft YaHei或Noto Sans CJK SC。这个坑不影响计算但影响出图质量血泪经验是每次新建脚本都先设这两行。6. 年龄模型验证与分类单元变化两个容易被忽略的进阶检查年龄模型验证这一步很多人直接跳过但它是整个时间序列分析的地基。如果原始年龄和模型年龄偏差太大后面所有“随时间变化”的结论都站不住。代码里用散点图对比Age (Ma)和model_age加一条 yx 参考线再算平均差异和最大绝对差异。我一般会额外算一个 Bland-Altman 图看差异是否随年龄均值变化如果差异有系统性趋势说明年龄模型存在非线性偏差。def validate_age_model(data, age_model_path): age_model pd.read_csv(age_model_path) if age_model_path.endswith(.csv) \ else pd.read_excel(age_model_path) merged data.merge(age_model, left_onSample_ID, right_onsample_id, howleft) merged[age_diff] merged[Age (Ma)] - merged[model_age] plt.figure(figsize(10, 6)) sns.scatterplot(datamerged, xAge (Ma), ymodel_age) plt.plot([merged[Age (Ma)].min(), merged[Age (Ma)].max()], [merged[Age (Ma)].min(), merged[Age (Ma)].max()], r--) plt.title(样本年龄与年龄模型比较) plt.xlabel(原始年龄(百万年)) plt.ylabel(模型年龄(百万年)) plt.show() print(f平均年龄差异: {merged[age_diff].mean():.2f} 百万年) print(f最大年龄差异: {merged[age_diff].abs().max():.2f} 百万年) return merged逻辑说明merge用left_onSample_ID和right_onsample_id注意两边列名大小写可能不同。age_diff为正说明原始年龄偏老为负说明偏年轻。参考线是 yx点越靠近红线说明模型越准。平均差异接近 0 不代表没问题还要看最大差异和差异的分布形态。分类单元变化分析是另一个进阶点。论文提到“物种数量变化”但代码框架里没有直接体现。常见做法是按样本统计不同形态聚类比如用 KMeans 对形态参数聚类的簇数量或者按形状因子阈值划分形态类型再看类型数量随年龄的变化。这一步没有标准答案取决于你对“物种”的定义——是有孔虫的生物学种还是形态种。我一般会先用 PCA 前两个主成分做 KMeans聚类数用肘部法或轮廓系数定然后画簇数量随年龄的折线图。环境数据整合那块代码里用了xarray读 NetCDF 和pd.merge_asof做最近邻匹配。merge_asof要求两边都按匹配键排序directionnearest表示找最近的年龄点。这里有个坑如果环境数据的时间分辨率和样本年龄分辨率差很多最近邻匹配会引入插值误差。更稳妥的做法是先对环境数据做线性插值到样本年龄再合并。另外merge_asof的tolerance参数可以限制最大匹配距离超过就置 NaN避免把相差几十个百万年的环境值硬匹配上去。从那以后我每次跑这类“形态 环境”的关联分析都强制先走一遍年龄模型验证和 ID 匹配率检查确认地基没问题再往下做统计。希望帮到你。本文还有配套的精品资源点击获取