STR混合图谱解析四问：分类建模、非负最小二乘与小波去噪

简介：面向2025深圳杯数学建模竞赛D题参赛者与数学建模学习者的完整分析论文PDF，围绕赛题给出从问题重述、问题分析、模型假设与符号定义到多问建模求解的系统框架，适合赛前梳理思路、赛后复盘及机器学习建模方法参考。资源包共1个PDF文件，约2.71MB，正文以论文形式整合模型推导、可运行代码与相关数据说明，目录覆盖特征工程设计、分类模型结构与数学表达、模型总体性能对比、不同贡献者数量表现，以及问题三的贡献者兼容性筛选、混合比例优化估计和组合匹配评分函数等模块。已有330人学习下载。读者可据此理清D题完整解题链路，理解分类模型构建、评估指标与优化匹配的实现细节，并借助论文中的代码和数据复现实验、对照结果，提升把建模思路转化为可验证方案的能力。

1. 半夜对着一个准确率下降的分类器时，这份 D 题论文能帮你解决什么问题

凌晨两点，我把一份混合样本的峰高表丢进分类器，2 人混合样本的 F1 冲到 0.99，换成 4 人混合样本后直接掉到 0.34，特征没改、模型没改、随机种子也没改。这是我第一次真正意识到，STR 混合图谱解析的核心难点不在于是选 GBDT 还是随机森林，而在于人数一多、等位基因一开始重叠，特征空间就失去了区分度。

这份某数学建模竞赛 D 题的完整分析论文，讨论的正是这个场景：从混合 STR 图谱里判断有几个人、各自占多少比例、还原每个人的基因型、再把噪声和拖尾峰清理掉。四个问题从判别到回归再到组合搜索，最后一问落在信号去噪上，形成一条完整链路。论文带模型推导、可运行代码和数据，适合正在备赛建模的同学，也适合想练一遍「多分类 + 非负最小二乘 + 小波阈值」这条技术路线的开发者。

2. 四问一条流水线：模型分层、技术栈与可执行环境

2.1 四个问题不是四篇独立论文，而是一条前后强耦合的流水线

很多人拿到这类论文会按章节单独阅读，读完四章仍不知道代码怎么串起来。这四个问题的实际依赖关系是单向的：问题一输出贡献者人数 N，问题二在 N 已知的前提下才能确定比例类别的维度，问题三需要 N 才能确定枚举组合的规模，问题四的去噪结果又反哺前三问的输入质量。任何一环判断错误，误差都会向后传播。

论文的结构设计恰好对应这条链路。问题一被建模成多类分类，标签是 2/3/4/5 人；特征层从等位基因结构、峰高统计、片段长度分布、基因座 one-hot 四个方向提取，模型层对比了随机森林、GBDT、MLP、SVM 四套方案。问题二把连续的比例估计离散成类别分类，按人数分别训练分类器——2 人用 GBDT、3 人用 MLP、4 到 5 人用随机森林。问题三转成组合搜索：先按兼容性把候选个体池筛到能枚举的规模，再用加权最小二乘解混合比例，最后用综合评分排序。问题四引入小波阈值去噪、拖尾峰规则过滤、自适应峰值筛选三段处理。

值得注意的是问题二的建模思路。严格来说比例是连续量，直接做回归更自然，但论文选择了离散类别，理由是样本量小、类别设计本身就是竞赛给定的网格（1:1、1:4、1:9 这类），当成分类问题能直接复用问题一的特征管线和评估指标，工程上更省事。这个取舍在数据量有限的场景下是合理的，代价是泛化到网格外的比例时会失效——这一点在二次开发部分值得再处理。

2.2 技术栈选型：为什么最终是 sklearn 配 PyWavelets 这一套

整个方案没有引入深度学习框架，模型层全部落在 scikit-learn 上。论文里 GBDT 的准确率是 0.9608，交叉验证 0.9584±0.0136，随机森林 0.7255，MLP 0.7549 但交叉验证只有 0.3603，SVM 0.5931。这组数字本身就说明了选型逻辑：样本量在千级，特征维度在二十上下，树模型在这个区间几乎总是赢家，而 MLP 训练不充分时交叉验证会大幅下滑，这正是它 0.36 那个数字的来源。

问题四的信号处理部分依赖 PyWavelets。小波分解、软阈值、重构这三个动作在 pywt 里各一个函数就能完成，不需要自己写滤波器组。数据处理和可视化走 pandas、numpy、matplotlib，交互调试用 Jupyter。这套组合的安装成本很低，一台普通笔记本十几分钟能配完。

提示：论文正文里的公式在 OCR 后有一定丢失，符号定义那一章尤其明显。拿代码对照公式读，比硬啃 PDF 快得多。

2.3 把 PDF 里的代码变成能跑的工程

论文给的是分析文档加代码片段，直接复制粘贴跑不通，得先补一层工程结构。我一般会按照「数据、特征、模型、评估」四个目录切开，让每一问对应一个 notebook，公共逻辑抽到 src 下面。这样改一个问题不会牵连其他三个。

# 建议的环境准备，Python 3.10 上下都行 python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy pandas scikit-learn scipy PyWavelets matplotlib jupyter

目录按下面这个形状铺开，能让后面每一章的代码都有明确去处：

str-mix/ ├── data/ │ ├── raw/ # 题目给的混合图谱峰表，一个样本一行或一个基因座一行 │ └── reference/ # 候选个体基因型库，问题三要用 ├── src/ │ ├── features.py # 特征构造函数，问题一、二共用 │ ├── train_q1.py # 人数判别训练脚本 │ ├── train_q2.py # 比例分类训练脚本 │ ├── search_q3.py # 兼容性筛选 + 最小二乘配比 │ └── denoise_q4.py # 小波阈值去噪 + 拖尾峰过滤 ├── notebooks/ # 四个问题各一个，负责画图和对比实验 └── outputs/ # 模型文件、评估报告、去噪前后对照图

把公共特征函数放进features.py而不是复制到两个 notebook 里，是我踩过坑之后养成的习惯。问题一和问题二的特征高度重叠，一旦两个地方各改一版，后面做消融实验时你根本说不清是哪份代码在跑。

3. 特征工程拆解：把峰高表变成特征向量的每一步

3.1 图谱原始字段长什么样

原始数据的最小单元是「一个样本 × 一个基因座」下的一组峰，每个峰带四个字段：等位基因编号（allele）、峰高（height）、片段长度（size）、是否 OL 标记（is_ol）。OL 指的是无法与标准梯度比对的片段，这类峰在论文里被单独计数，因为它们的数量本身就和混合复杂度相关。

论文的观察是，单个基因座下正常情况最多两个等位基因，超过两个就基本可以判定为多人混合。这个判断构成了整个特征工程的地基——allele_count和expected_allele_ratio之所以在特征重要性排序里排到前两位，就是因为它们直接刻画了这个结构信号。峰高则服从明显的幂律分布，大多数峰集中在低位，少数主峰高度突出，主次峰比值（top1/top2）能反映混合是否均衡。

数据分布上还有一个细节值得记住：论文统计的样本里，2 人和 3 人混合各占 240 和 256 个，4 人 176 个，5 人 144 个。类别本身就不均衡，所以评估时必须看加权 F1 和混淆矩阵，只看 accuracy 会被多数类带偏。

3.2 特征构造函数的逐段拆解

特征构造是整个方案里最值得抄的一段代码，因为它同时服务问题一和问题二，改动一处就能影响两个任务的输入。下面这个版本按论文的特征清单重写，把「训练时可用的特征」和「预测时可用的特征」分开处理。

import numpy as np import pandas as pd # 16 个核心基因座，按论文顺序编号，避免直接写死具体位点名 CORE_LOCI = [f"Locus{i:02d}" for i in range(1, 17)] def locus_features(peaks: pd.DataFrame, n_contrib: int | None = None) -> dict: """ peaks: 单个基因座下的峰表，必须含 allele / height / size / is_ol 四列 n_contrib: 真实贡献者人数。只在离线造训练集时传入； 线上预测时必须传 None，否则会造出标签泄漏特征 """ h = peaks["height"].to_numpy(dtype=float) h = h[h > 0] # 零峰视为未检出，先剔掉 h_sorted = np.sort(h)[::-1] # 降序，方便取 top-k feats = {} # 1) 结构特征：等位基因数量类 feats["allele_count"] = len(peaks) feats["non_ol_count"] = int((peaks["is_ol"] == 0).sum()) feats["ol_ratio"] = 1.0 - feats["non_ol_count"] / max(len(peaks), 1) # 2) 峰高统计特征：均值、离散度、极值 feats["height_mean"] = float(h.mean()) if len(h) else 0.0 feats["height_std"] = float(h.std(ddof=0)) if len(h) else 0.0 feats["height_cv"] = feats["height_std"] / max(feats["height_mean"], 1e-6) feats["height_min"] = float(h.min()) if len(h) else 0.0 feats["height_max"] = float(h.max()) if len(h) else 0.0 feats["height_range"] = feats["height_max"] - feats["height_min"] # 3) 主次峰比值：判断是否存在单方主导 top = np.pad(h_sorted[:4], (0, max(0, 4 - len(h_sorted))), constant_values=0.0) feats["height_top1"] = float(top[0]) feats["height_top4"] = float(top[3]) feats["top1_top2_ratio"] = float(top[0] / top[1]) if top[1] > 0 else 10.0 feats["height_ratio_mean"] = float(top[:3].mean() / max(top[0], 1e-6)) # 4) 片段长度分布：混合越复杂，size 展布通常越宽 s = peaks["size"].to_numpy(dtype=float) feats["size_mean"] = float(s.mean()) if len(s) else 0.0 feats["size_std"] = float(s.std(ddof=0)) if len(s) else 0.0 feats["size_range"] = float(s.max() - s.min()) if len(s) else 0.0 # 5) 这个特征只在训练期存在，线上拿不到 n_contrib if n_contrib is not None: feats["expected_allele_ratio"] = feats["allele_count"] / (2.0 * n_contrib) return feats

这段代码的关键点在最后那个if。expected_allele_ratio的计算式是「观测等位基因数 ÷ 理论等位基因数」，而理论值2 × n_contrib里直接含了标签 N。训练时它能让模型准确率冲到 0.96，上线时这个输入根本不存在，模型性能会大幅下降。论文的特征重要性图里它权重超过 0.5，排在第一，恰恰说明模型学到的主要是这个泄漏信号，而非真实的图谱规律。

参数上还有两处要留意。top1_top2_ratio在次高峰为 0 时会给一个默认值 10，这个哨兵值如果设得太小，会和真实的高比值样本混在一起；设成 10 在多数 STR 数据里够用，但换数据集时最好重新看一下分布。height_cv的分母做了 1e-6 的保护，防止全零峰导致 NaN，这类 NaN 一旦漏进管线，树模型会直接报错。

3.3 问题二的比例标签：把连续占比离散成类别

问题二把比例估计做成分类，标签得自己从数据里还原。论文给出的类别是按人数变化的网格：2 人样本有 1:1、1:4、1:9 三类，3 人有 6 类，4 人和 5 人各 7 到 8 类。实际数据里比例可能带小数误差，需要按最近邻映射回设计网格。

def ratio_label(shares, grid): """ shares: 一串占比，例如 [0.2, 0.8] grid: 该人数下所有合法比例模板，形如 [[1,1],[1,4],[1,9]] 返回: 最接近的模板对应的字符串标签 """ x = np.asarray(shares, dtype=float) x = x / x.sum() # 先归一化，避免采集误差 best, best_err = None, np.inf for g in grid: t = np.asarray(g, dtype=float) t = t / t.sum() err = np.abs(np.sort(x) - np.sort(t)).sum() # 排序后比较，规避排列歧义 if err < best_err: best, best_err = g, err return ":".join(str(v) for v in best)

这里把两个向量都排序后再比，是为了规避「谁排第一位」带来的歧义。混合样本里贡献者本身没有顺序，1:4 和 4:1 是同一个类别，如果不排序就生成标签，同一类会被拆成两个，类别数直接翻倍。

标签生成完之后，训练要按人数分组做。论文分别给 2 人、3 人、4/5 人挑了三套模型，理由是不同人数下的特征分布差异太大，混在一起训一个多类分类器会互相干扰。这个做法在小样本场景下确实有效，但每个分组只有一百到两百条样本，交叉验证的方差会很大，报告指标时最好把折间标准差一起写出来。

4. 问题三与问题四的实现：兼容性搜索、最小二乘配比与小波去噪

4.1 兼容性评分：先把候选组合砍到能枚举的规模

问题三的输入是候选个体基因型库，输出是「哪几个人组成了这个混合样本」。如果库里有一百人、要找 3 人组合，枚举量是十六万量级，逐个算峰高拟合不现实。论文的做法是先算兼容性，把明显不可能的个体剔掉，再在小集合上枚举。

兼容性评分的定义很直接：某个候选人的基因型，在该基因座上能否覆盖观测到的等位基因。完全覆盖记满分，部分覆盖扣分，完全对不上直接出局。把所有基因座的得分汇总，取前若干名进入组合搜索。

def compatibility(cand_geno: dict, observed: dict) -> float: """ cand_geno: {基因座: (等位基因a, 等位基因b)}，纯合则 a == b observed: {基因座: set(观测到的等位基因)} 返回: 0~1 之间的兼容度，按基因座取平均 """ scores = [] for loc, obs in observed.items(): g = cand_geno.get(loc) if g is None: scores.append(0.0) continue allow = set(g) # 该候选人可能贡献的等位基因集合 hit = len(allow & obs) scores.append(hit / max(len(allow), 1)) return float(np.mean(scores)) if scores else 0.0

这段代码的判分逻辑是「分母用候选人自己的基因型规模」，不是用观测集合的规模。这一点容易被写反：如果用观测集合做分母，等位基因多的混合样本会把所有候选人的分数都压低，排序结果被稀释；用候选人自己的基因型做分母，纯合子（两个相同等位基因）和杂合子就能在同一个尺度上比较。

阈值怎么定是个经验值。我一般会按分数排序取前 15 到 20 个候选进入枚举，3 人组合大约一千组，算下来可接受。如果库里人数更多，可以改成贪心：先定兼容度最高的一个，再在剩下的里面找与之互补度最高的，牺牲一点最优性换速度。

4.2 用非负最小二乘解混合比例

候选组合确定之后，剩下的问题是「这组人各自占多少」。论文的假设是峰高与 DNA 含量近似线性相关，于是可以写成线性方程组：观测峰高 ≈ 各贡献者在该位点的等位基因拷贝数 × 各自比例 × 一个缩放系数。用最小二乘解比例，再归一化。

from scipy.optimize import nnls import numpy as np def design_matrix(alleles, genotypes): """ alleles: 该基因座下观测到的等位基因列表 genotypes: 候选人的基因型列表，每个形如 (a, b) 返回: 矩阵 G，形状 (len(alleles), n_candidates) G[i][k] = 第 k 个候选人在 alleles[i] 上的拷贝数（0/1/2） """ G = np.zeros((len(alleles), len(genotypes))) idx = {a: i for i, a in enumerate(alleles)} for k, (a, b) in enumerate(genotypes): if a in idx: G[idx[a], k] += 1 if b in idx: G[idx[b], k] += 1 return G def fit_shares(heights, alleles, genotypes): """heights 与 alleles 一一对应；返回归一化后的比例向量""" G = design_matrix(alleles, genotypes) w, _ = nnls(G, np.asarray(heights, dtype=float)) # 非负约束，杜绝负比例 return w / w.sum() if w.sum() > 0 else w

之所以用nnls而不用普通的lstsq，是因为普通最小二乘在噪声存在时会解出负的比例值，物理上无意义，还得手工截断，截断后又破坏了残差最小的性质。非负最小二乘直接把这个约束加进去，代价是求解慢一点，但在这个规模下可以忽略。

解出比例后，论文还设计了一个综合匹配评分，把「等位基因解释率」和「峰高分布相关性」两件事合到一起：解释率衡量这个组合能不能覆盖所有观测峰，相关性衡量拟合出来的峰高曲线和实际曲线像不像。前者防止漏掉贡献者，后者防止比例解偏。两者加权求和排序，取最高分作为最终答案。

注意：解释率高但相关性低的组合很常见，典型情况是把一个次要贡献者的等位基因也算进去了，覆盖率上去了，峰高却对不上。加权时相关性那一项别设太小，我的经验是两项接近等权比较稳。

4.3 小波阈值去噪与拖尾峰过滤

问题四处理的是输入质量。混合图谱里的噪声有两类：一类是高频随机噪声，体现在峰高的小幅抖动上；另一类是 PCR 扩增产生的拖尾峰（stutter），特征是紧挨着主峰、片段长度差 2 到 4 bp、高度明显偏低。

小波阈值去噪负责第一类。做法是把信号做多层小波分解，对高频细节系数做软阈值收缩，再重构回去。阈值取 σ√(2lnN) 这个经典形式，σ 用第一层细节系数的中位数绝对偏差估计。

import pywt import numpy as np def wavelet_denoise(y, wavelet="db4", level=3): """ y: 一条峰高序列（按片段长度排序） wavelet: 小波基，db4 在保留尖峰形状上比 haar 好 level: 分解层数，信号短于 2**level 时要往下调 """ coeffs = pywt.wavedec(y, wavelet, level=level) sigma = np.median(np.abs(coeffs[-1])) / 0.6745 # 用最细层估计噪声水平 thr = sigma * np.sqrt(2 * np.log(len(y))) # 第 0 个是近似系数，不动；只收缩细节层 coeffs[1:] = [pywt.threshold(c, thr, mode="soft") for c in coeffs[1:]] recon = pywt.waverec(coeffs, wavelet) return recon[: len(y)] # 重构后长度可能多出 1~2 点

三个参数里最容易出问题的是level。分解层数过深，低频近似系数被过度压缩，主峰高度会被削平，后面的比例估计全跟着偏；层数过浅则噪声去不干净。论文的经验是 3 层，这个值跟信号长度有关，序列点数是 2 的幂次以下时，层数别超过 log2(长度)-1。

拖尾峰的过滤走规则而不是滤波。论文给的判据是峰高比值加片段长度差：与主峰的片段长度差在 2 到 4 bp 之内，且峰高比低于某个阈值（常见取 0.15 上下），就判为 stutter 剔掉。这个规则简单但有效，前提是同一基因座内的峰已经按片段长度排好序。

最后一步是自适应峰值筛选，对候选峰做质量评分——峰高归一化后加上与相邻峰之间的间距，低于阈值的峰不进后续匹配。三段处理串起来之后，论文报告在多数样本上贡献者识别准确率能到 0.8 至 1.0，极端样本上的表现也优于只做预处理或传统去噪的方案。

5. 避坑与排查：这套代码跑起来最容易出问题的五个地方

5.1 训练阶段的两个隐性问题

现象一：离线准确率 0.96，换个数据集直接掉到 0.6 以下。原因是expected_allele_ratio这个特征在训练期用到了真实人数。它本质上是把标签做了个除法再喂回模型，树模型只要在这个维度上切一刀就能分类。排查方法很简单：把所有含n_contrib的特征列删掉重训一遍，如果准确率大幅下降，说明模型之前主要靠泄漏特征。解决方式是在locus_features里强制n_contrib=None生成线上特征，训练集也只保留这一版，泄漏特征单独拿去做错误分析，不进模型。

现象二：MLP 的训练准确率 0.75，交叉验证只有 0.36。这是论文里真实出现的数字，典型的小样本过拟合。MLP 参数量相对千级样本太多，默认配置下几轮就把训练集背下来了。排查时看训练损失和验证损失的曲线是否分叉，分叉就是过拟合。解决方式是砍隐藏层宽度、加 dropout、加 L2，或者干脆换成树模型——在特征维度只有二十上下、样本量不大的情况下，GBDT 的性价比确实更高。

现象三：混淆矩阵里 4 人样本大量被判成 3 人。这是等位基因重叠导致的，不是代码 bug。人数越多，重复等位基因越多，观测到的allele_count增长会饱和，模型分不清 3 人和 4 人。排查方式是画出不同人数下allele_count的分布箱线图，如果 3 人和 4 人的箱子大面积重叠，就说明这个特征本身已经到上限了。解决方向是引入峰高比例的细粒度特征（各峰占该位点总峰高的比例），或者引入相邻基因座的联合特征，让模型有更多区分依据。

5.2 数据与工程层面的三个坑

现象四：某几个基因座的特征全是 NaN，训练直接报错。常见原因是该基因座在所有样本里都没检出峰，height_mean除零后变成 NaN，height_cv再一除就变成 inf。排查时对特征矩阵跑一遍df.isna().sum()和np.isinf检查。解决方式是在locus_features里对空峰表直接返回全零字典，并在管线最前面加一层SimpleImputer(strategy="constant", fill_value=0)保底。

现象五：小波去噪后主峰高度普遍变矮，比例估计整体偏移。原因是分解层数超过了信号实际能支撑的长度，或者软阈值把主峰的细节系数一起收缩了。排查方式是画出去噪前后的对照曲线，重点看最高的那几个峰有没有被削。解决方式是把level降到 2 或者换mode="garrote"，同时限定只对最细的两层做阈值处理，近似系数层完全不动。

提示：去噪前后一定要存图。这个环节的失败是静默的——程序不报错，指标略微变差，你不看曲线根本发现不了。

除了上面五条，还有一个容易忽视的点：问题三的候选个体池里如果混进了亲缘关系较近的个体，兼容性评分会集体偏高，排序失效。这属于模型假设被打破的情况，论文里明确假设了贡献者之间基因型独立、不存在亲缘关系，超出这个假设的数据不该硬套。

6. 二次开发：把 0.96 的准确率拆开重装一遍

拿到这套代码之后，我做的第一件事是把特征表按「是否用到标签」分成两栏，逐列过一遍。这一步不改模型，只改输入，但决定了后面所有指标是否可信。

删掉泄漏特征后重训，准确率肯定会下降，这个过程不必追求高指标。指标下降说明模型之前确实靠了那条捷径，现在学到的东西才对应真实的图谱规律。

第二个改造方向是把问题二的离散分类换成回归加约束。原来按人数分组、每组挑一个分类器，本质上是用分组绕开了「比例」这个连续量的建模难度。改成回归之后，输出一个 N 维的占比向量，约束非负、和为 1，样本利用率比分组分类高得多，也能处理网格外的比例。代价是需要重新设计评估指标，MAE 比 accuracy 更能反映实际误差。

第三个方向是验证流程。原有方案用的是普通交叉验证，同一个样本的不同基因座可能同时出现在训练集和验证集里，导致指标虚高。改成按样本分组、按基因座分层，即用GroupKFold保证同一样本不跨折，评估结果才贴近真实场景。

from sklearn.model_selection import GroupKFold, cross_val_score from sklearn.ensemble import GradientBoostingClassifier # groups 是样本 ID，同一样本的所有基因座行共享同一个 group cv = GroupKFold(n_splits=5) scores = cross_val_score( GradientBoostingClassifier(n_estimators=200, max_depth=3, random_state=0), X_clean, y, groups=groups, cv=cv, scoring="f1_weighted" ) print(scores.mean(), scores.std()) # 这个数字才是能写进报告的

从那以后，我每次拿到一份带「准确率 0.96」的分析文档，都强制先做一遍泄漏排查：把所有特征列按「预测时能不能拿到」过一遍，删掉拿不到的，再重跑指标。这一步通常只花二十分钟，却能避免后面两周基于不真实指标做出来的错误设计。这套 D 题论文的价值不在于那个 0.96，而在于它把四个问题的完整链路、特征设计、去噪思路和失败案例都留了下来，改起来有据可依。希望帮到你。

题外一句，论文正文里还有几处 OCR 断字，公式部分对不上的时候，直接对着代码里的实现倒推，比反复读文字快。

本文还有配套的精品资源，点击获取