2012研究生数学建模题目

DNA序列内含子(Intron)

ADNA序列内含子(Intron)

基因识别问题及其算法实现

一、背景介绍DNA是生物遗传信息的载体，其化学名称为脱氧核糖核酸（Deoxyribonucleicacid，缩

写为DNA）。DNA分子是一种长链聚合物，DNA序列由腺嘌呤（Adenine,A），鸟嘌呤（Guanine,

G），胞嘧啶（Cytosine,C），胸腺嘧啶（Thymine,T）这四种核苷酸（nucleotide）符号按一定

的顺序连接而成。其中带有遗传讯息的DNA片段称为基因（Gene）（见图1第一行）。其他的

DNA序列片段，有些直接以自身构造发挥作用，有些则参与调控遗传讯息的表现。

在真核生物的DNA序列中，基因通常被划分为许多间隔的片段（见图1第二行），其中编

码蛋白质的部分，即编码序列（CodingSequence）片段，称为外显子（Exon），不编码的部

分称为内含子（Intron）。外显子在DNA序列剪接（Splicing）后仍然会被保存下来，并可在

基因（Gene）

o

外显子(Exon)

图1真核生物DNA序列（基因序列）结构示意图

蛋白质合成过程中被转录（transcription）、复制（replication）而合成为蛋白质（见图2）。

DNA序列通过遗传编码来储存信息，指导蛋白质的合成，把遗传信息准确无误地传递到蛋

白质（protein）上去并实现各种生命功能。

基因（Gene）

DNA序列

剪接、转录、复制蛋白质序列

图2蛋白质结构示意图

对大量、复杂的基因序列的分析，传统生物学解决问题的方式是基于分子实验的方法，

其代价高昂。诺贝尔奖获得者W.吉尔伯特（WalterGilbert，1932—；【美】，第一个制备

出混合脱氧核糖核酸的科学家）1991年曾经指出：现在，基于全部基因序列都将知晓，并

以电子可操作的方式驻留在数据库中，新的生物学研究模式的出发点应是理论的。一个科学

家将从理论推测出发，然后再回到实验中去，追踪或验证这些理论假设。随着世界人类基

1

[4]。I{A,T,G,C}，长度（即核苷酸符号个数，又称碱基对（BasePair1]。现对于任意确定的bI1,S[n]b[0],ub(bI[4]。I{A,T,G,C}，长度（即核苷酸符号个数，又称碱基对（BasePair1]。现对于任意确定的bI1,S[n]b[0],ub(bI)。S[n]b，n0,1,2,N1

息，对生物学、医学、药学等诸多方面都具有重要的理论意义和实际价值，也是目前生物信

息学领域的一个研究热点。

二、数字序列映射与频谱3-周期性：

对给定的DNA序列，怎么去识别出其中的编码序列（即外显子），也称为基因预测，是

一个尚未完全解决的问题，也是当前生物信息学的一个最基础、最首要的问题。

基因预测问题的一类方法是基于统计学的[1]。很多国际生物数据网站上也有“基因识别”

的算法。比如知名的数据网站/GENSCAN.html提供的基因识别软件

GENSCAN（由斯坦福大学研究人员研发的、可免费使用的基因预测软件）,主要就是基于隐马

尔科夫链（HMM）方法。但是，它预测人的基因组中有45000个基因，相当于现在普遍认可

数目的两倍。另外，统计预测方法通常需要将编码序列信息已知的DNA序列作为训练数据

集来确定模型中的参数，从而提高模型的预测水平。但在对基因信息了解不多的情况下，基

因识别的准确率会明显下降。

因此在目前基因预测研究中，采用信号处理与分析方法来发现基因编码序列也受到广泛

重视

1.数字序列映射

在DNA序列研究中，首先需要把A、T、G、C四种核苷酸的符号序列，根据一定的规

则映射成相应的数值序列，以便于对其作数字处理。

令）长度，单位

记为bp）为N的任意DNA序列，可表达为

S{S[n]|S[n]I,n0,1,2,N1}

即A、T、G、C的符号序列S：S[0],S[1],,S[N，令

ub[n]0,

称之为Voss映射[5]，于是生成相应的0-1序列（即二进制序列）{ub[n]}：u

u[N1]b

例如，假设给定的一段DNA序列片段为S=ATCGTACTG，则所生成的四个0-1序列分

别为：

2

；；Nb[T5000050000k{u[n]}：{0,0,0,1,0,0,0,0,1{u[n]}：{0,1,0,0,1,0,0,1,0}1uG100100N3n]eC200200处，具有较大的频谱峰值(Peakj2nkN(2)300300,400400k0,1,,N1500500(1)；；Nb[T5000050000k{u[n]}：{0,0,0,1,0,0,0,0,1{u[n]}：{0,1,0,0,1,0,0,1,0}1uG100100N3n]eC200200处，具有较大的频谱峰值(Peakj2nkN(2)300300,400400k0,1,,N1500500(1)600600A

{u[n]}：{0,0,1,0,0,0,1,0,0}C这样产生的四个数字序列又称为DNA序列的指示序列(indicatorSequence)。2.频谱3-周期性

为研究DNA编码序列（外显子）的特性，对指示序列分别做离散Fourier变换(DFT)

Ub[k]n0

以此可得到四个长度均为N的复数序列{Ub[k]}，bI。计算每个复序列{Ub[k]}的平方

功率谱，并相加则得到整个DNA序列S的功率谱序列{P[k]}：

P[k]U[k]2U[k]2U[k]2U[k]2,k0,1,N1A对于同一段DNA序列，其外显子与内含子序列片段的功率谱通常表现出不同的特性

10000

(Pk)

0

k

10000

P()k

0

k

图3编号为BK006948.2的酵母基因DNA序列的功率谱（因为对称性，实际这里只给出了功率谱图

的一半）。(a)上图是基因上一段外显子(区间为[81787，82920]，长1134bp)对应的指示序列映射的功率

谱，它具有3-周期性；(b)下图是基因上一段内含子（区间为[96361,97551]，长1191bp）的指示序列的功

率谱，它不具有3-周期性。

可以看到：外显子序列的功率谱曲线在频率

Value)，而内含子则没有类似的峰值。这种统计现象被称为碱基的3-周期(3-basePeriodicity)

[2][3]。

记DNA序列S的总功率谱的平均值为

3

k0ER0n0,1,2,I{A,T,G,C}出现在该序列的0,3,6,...N－3与1,4,7,…N－2以及2,5,8,…N－1等位置上,yzbbIbIn0bbIx,y,zbbN3NN3(4)2），是N1}中，若N为3的倍数，将核苷酸符号N3xyN3(3)U[]k0ER0n0,1,2,I{A,T,G,C}出现在该序列的0,3,6,...N－3与1,4,7,…N－2以及2,5,8,…N－1等位置上,yzbbIbIn0bbIx,y,zbbN3NN3(4)2），是N1}中，若N为3的倍数，将核苷酸符号N3xyN3(3)U[]uej23zN3j2bIn0232

bbbbbbbbN1N(xyzxyxzyz2nN3u[n]ej23b2N122n2E

而将DNA序列在特定位置，即k处的功率谱值，与整个序列S的总功率谱的平均值的

比率称为DNA序列的“信噪比”（SignalNoiseRatio，SNR），即

P[N3]R

DNA序列的信噪比值的大小，既表示频谱峰值(PeakValue)的相对高度，也反映编码或

非编码序列3-周期性的强弱。

信噪比大于某个适当选定的阈值R（比如RDNA

显子）通常满足的特性，而内含子则一般不具有该性质[6]。

在DNA序列{S[n],

b

的频数分别记为x和，则

2P[N3]bb

bbI

易见，当四种核苷酸符号b（bI）在序列的上述第一、第二、第三个子序列上出现的频

数越接近相等时，

曲线，在频率处具有较大的频谱峰值(PeakValue)，反映了在基因外显子片段上，四种核

苷酸符号在序列的三个子序列上分布的“非均衡性”。通常认为这种现象源于编码基因序列

“密码子”（coden）使用的偏向性（bias）。虽然目前对此现象产生的“机理”还不是十分

地清楚，但是频谱的3-周期性被普遍认为是可用于识别基因编码序列（外显子）的一个重要

的特征信息。

3.基因识别

频谱峰值特征的发现，或者频谱与信噪比概念的引入，其最终目的是要探测、预报一个

4

值数值化变换0,1,2,N1。nN1），在以ub[i]eM33M3p(n;)0,1,2,N1，经过标准化处理（即除以最大频谱值M3DFT信噪比计算取长度M（通M12inM1p(n;]U[TM3，并画出其频谱曲线功率谱或判别分类M122M23外显子]j2ikM)M

G预测结果,，即]值数值化变换0,1,2,N1。nN1），在以ub[i]eM33M3p(n;)0,1,2,N1，经过标准化处理（即除以最大频谱值M3DFT信噪比计算取长度M（通M12inM1p(n;]U[TM3，并画出其频谱曲线功率谱或判别分类M122M23外显子]j2ikM)M

G预测结果,，即]U[k0,1,,M123MC]U[23MM]p(n;)23

阈

DNA序列映射

图4基于序列频谱3—周期性的的基因预测方法流程图

已经有一些研究者提出了识别基因的算法（如参见[6]及其后面的文献）。目前利用信噪

比的基因识别算法通常有两种：一是固定长度窗口滑动法[2][3]；另一是移动信噪比曲线识别

法[6]。

基于固定长度滑动窗口上频谱曲线的基因识别方法:

对一个DNA序列S和它的指示序列{ub[n]}，bI，n

常取为3的倍数，例如M=99,129,255,513等）作为固定窗口长度。

对任意n（0n为中心的长度为M的序列片段[n，n

上（当n接近序列的两端时，窗口实际有效长度可能会小于M），作四个指示序列的离散

Fourier变换(DFT)

Ub[k]

inM12

并求出它在处总频谱

P[M3]U[A

把这样得到的频谱值，n

max{p(n;)}）0nN1

5

0.335004000450050005500600065007000

pp(n0.335004000450050005500600065007000

pp(n;)

p(n;)[n]}，bI0,1,2,N1。对任意nN1），通常in10,1,,NU[n]U[]U[]U[ATGCE[n]E[n]P[k]E[n]为移动子序列SE[n]0~n1M3M3ub[i]e22223333]k0j2ikMnnn

，0n,。在坐标系中画出移动序k0,1,,n1

0.9

0.8

0.7n)p(umectrspNAD

0.2

0.1

03000

nucleotidepositionn

图5固定长度滑动窗口的频谱曲线（人类线粒体基因，NC_012920_1.fasta）

图中红色水平细线条是DNA序列实际的基因外显子的区间。滑动窗口频谱曲线的

峰与基因外显子区间具有“对应”关系。

基于DNA序列上“移动序列”信噪比曲线的基因识别方法：

设已知DNA序列S和它的指示序列{u，n

（0n取3的倍数并逐渐增大。在n的左边一个长度为n的序列片段[0，

n-1]上，相应的子序列S称为DNA序列S

个指示序列的离散Fourier变换(DFT)

Ub[k]i0

并求出移动子序列S，n

P[n3]R[n]nN1

n1其中的功率谱的平均值

列S的信噪比曲线R[n]（称为信噪比移动曲线（SNRwalkcurve），见图

6

8(NR3500DNA移动序列其指示序列的信噪比曲线。（人类线粒体基因，NC_012920_1.fasta）A0,C1,G2,T8(NR3500DNA移动序列其指示序列的信噪比曲线。（人类线粒体基因，NC_012920_1.fasta）A0,C1,G2,T3，也给出功率谱与信噪比264000450050005500600065007000

12

10

n)R

S

4

2

30000nucleotidepositionn

图6

图中红色水平细线条是DNA序列实际的基因外显子的区间。DNA序列的信噪比移动曲线

的峰、谷与基因外显子区间的端点也具有较“明显的”的对应关系。

三、请研究的几个问题：

1.功率谱与信噪比的快速算法

对于很长的DNA序列，在计算其功率谱或信噪比时，离散Fourier变换(DFT)的总体计

算量仍然很大，会影响到所设计的基因识别算法的效率。大家能否对Voss映射，探求功率

谱与信噪比的某种快速计算方法？

在基因识别研究中，为了通过引入更好的数值映射而获取DNA序列更多的信息，除了

上面介绍的Voss映射外，实际上人们还研究过许多不同的数值映射方法。例如，著名的

Z-curve映射（参见[5]或者附件1）。试探讨Z-curve映射的频谱与信噪比和Voss映射下的频

谱与信噪比之间的关系；

此外，能否对实数映射，如：

的快速计算公式？

2.对不同物种类型基因的阈值确定

对特定的基因类型的DNA序列，将其信噪比R的判别阈值取为R0，带有一定的主

观性、经验性。对不同的基因类型，所选取的判别阈值也许应该是不同的。附件中给出了来

自于著名的生物数据网站：/guide/的几个基因序列数据，另外也

给出了带有编码外显子信息的100个人和鼠类的，以及200个哺乳动物类的基因序列的样本

数据集合。大家还可以从生物数据库下载更多的数据，找你们认为具有代表性的基因序列，

并对每类基因研究其阈值确定方法和阈值结果。此外，对按照频谱或信噪比特征将编码与非

7

编码区间分类的有效性，以及分类识别时所产生的分类错误作适当分析。

3.基因识别算法的实现

我们的目的是要探测、预报尚未被注释的、完整的DNA序列的所有基因编码序列（外

显子）。目前基因识别方面的多数算法结果还不是很充分。例如前面所列举的某些基因识别

算法，由于DNA序列随机噪声的影响等原因，还很难“精确地”确定基因外显子区间的两

个端点。

对此，你的建模团队有没有更好的解决方法？请对你们所设计的基因识别算法的准确率

做出适当评估，并将算法用于对附件中给出的6个未被注释的DNA序列（gene6）的编码

区域的预测。

4.延展性研究

在基因识别研究中，还有很多问题有待深入探讨。比如

（1）采用频谱或信噪比这样单一的判别特征，也许是影响、限制基因识别正确率的一

个重要原因。人们发现，对某些DNA序列而言，其部分编码序列（外显子），尤其是短的（长

度小于100bp）的编码序列，就可能不具有频谱或者信噪比显著性。你们团队能否总结，甚

至独自提出一些识别基因编码序列的其它特征指数，并对此做相关的分析？

（2）“基因突变是生物医学等方面的一个关注热点。基因突变包括DNA序列中单个核

苷酸的替换，删除或者插入等。那么，能否利用频谱或信噪比方法去发现基因编码序列可能

存在的突变呢？

上面提出的基于频谱3-周期性的基因预测四个方面问题中，“快速算法”与“阈值确定”

是为设计基因预测算法做准备的。此外，在最后的延展性研究中，各队也可以对你们自己认

为有价值的其它相关问题展开探讨。

参考文献：

【1】Burge,C.,Karlin,S.,1997.PredictionofcompletegenestructuresinhumangenomicDNA.

J.Mol.Biol.268,78–94.

【2】Anastassiou,D.,2000.Frequency-domainanalysisofbiomolecularsequences.

Bioinformatics16,1073–1081.

【3】Kotlar,D.,Lavner,Y.,2003.Genepredictionbyspectralrotationmeasure:anewmethodfor

identifyingprotein-codingregions.GenomeRes.13,1930–1937.

8

th

【4】Berryman,M.J.,Allison,A.,2005.Reviewofsignalprocessingingenetics.Fluctuationandth

NoiseLetters.5(4),13-35.

【5】Sharma,S.D.,Shakya，K.,Sharma，S.N.,2011.EvaluationofDNAMappingSchemesfor

ExonDetection.InternationalConferenceonComputer,CommunicationandElectrical

Technology–ICCCET.2011,18th&19

【6】Yin,C.,Yau,S.S.-T.2007.Predictionofproteincodingregionsbythe3-baseperiodicity

analysisofaDNAsequence.JournalofTheoreticalBiology.247,687–694.

B【2012第九届全国研究生数学建模竞赛B题】基于卫星无源探测的空间飞行器主动段轨道估计与误差分析

有些国家会发射特殊目的的空间飞行器，如弹道式导弹、侦察卫星等。对他国发射具有敌意的空间飞行器实施监控并作出快速反应，对于维护国家安全具有重要的战略意义。发现发射和探测其轨道参数是实现监控和作出反应的第一步，没有观测，后续的判断与反应都无从谈起。卫星居高临下，是当今探测空间飞行器发射与轨道参数的重要平台。观测卫星按轨道特点，可分为高轨地球同步轨道卫星和中低轨近圆轨道卫星。其中同步轨道距地球表面约3.6万千米，轨道平面与地球赤道平面重合，理论上用3颗间隔120度分布的同步轨道卫星可覆盖地球绝大部分表面。中低轨近圆轨道距地球表面数百到几千千米不等，根据观测要求，其轨道平面与赤道平面交成一定角度，且常由若干颗卫星实现组网探测。装置于卫星上的探测器包括有源和无源两类：有源探测器采用主动方式（如雷达，激光）搜寻目标，同时具备定向和测距两种能力；无源探测器则被动接收目标辐射。采用无源探测器的观测卫星常采用红外光学探测器，只接收目标的红外辐射信息，可定向但不能测距。对于火箭尾部喷焰的高度敏感性是红外技术的长处，但易受气候影响与云层干扰则是其缺点。探测的目的是为了推断空间飞行器的轨道参数，推断是基于观测数据并通过数学模型与计算方法作出的。当观测卫星飞行一段时间，探测器测得目标相对于运动卫星的观测数据，以观测卫星和空间飞行器的运动模型和观测模型为基础，对空间飞行器的轨道参数（包括轨道位置、速度初值和其他模型参数）进行数学推断，为飞行器类别、飞行意图的判断提供信息基础。

9

vBAYEQ经线x

yCD赤道r(t)P椭圆轨道

vBAYEQ经线x

yCD赤道r(t)P椭圆轨道惯性飞行段和再入大气层后的攻击段。主动段通常由多级火箭相继推进，前一级火箭完成推进后脱落，由后一级火箭接力。惯性飞行段在空气阻力极小的大气层外，靠末级火箭关机前获得的速度在椭圆轨道上作无动力惯性飞行。攻击段则根据任务需求，受控制后再入大气层，飞向目标。对于卫星而言，在其寿命结束前一直绕地飞行，故无攻击段。图1是空间飞行器的主动段示意图（未按实际比例）。主动段又可细分为若干子段：垂直上升段，程序拐弯段和重力斜飞段。按最优轨道设计，为节约燃料，箭体应尽快穿过稠密大气层，故火箭一般先垂直发射。设A点为地面发射点，AB为垂直上升段，BC弧段为程序拐弯段，CD弧段为重力斜飞段，DE弧段为椭圆轨道。程序拐弯段连接垂直上升段与重力斜飞段，在外力矩控制下使箭体转过一定角度，该段完成后外加力矩撤销，进入斜飞状态。第一级火箭通常负担“垂直段+程序拐弯段（加外力矩）+重力斜飞段的前段”的推进（视发动机的特性），重力斜飞段的后程则靠第二、第三级火箭相继完成。由于斜飞状态下地球引力与推力不在同一直线，所以箭体质心的运动轨迹为带一定弧度的光滑曲线。

垂直段ZC

北极vr(t)

纬线C

o

赤道平面

XC

南极

图1空间飞行器主动段轨道的示意图为描述观测卫星和空间飞行器的运动，需要建立适当的坐标系。本题基础坐标系为随地心平移的坐标系，取地球中心O为原点，地球自转轴取为z轴，指向北极为正向，轴c由O指向零时刻的0经度线，再按右手系确定轴，建立直角坐标系OXYZ。地心Occcccc在绕日椭圆轨道上运动，所以理论上OXYZ系是非惯性系。但地球公转周期远大于空cccc间飞行器的观测弧段时长，故本题在短时间内认定该系为惯性坐标系，该基础坐标系不随地球旋转。

10

观测坐标系示意图FF|r(tm)|r(t)v

m(t)为

对时间t观测坐标系示意图FF|r(tm)|r(t)v

m(t)为

对时间t的二阶导数，即加速度；为地球引力常数（本题中地球引力常数取

2

F|r(tm)|vyzxsxm(t)

m

GsGm(t)

（2）（3）（1）第二个坐标系是随卫星运动的观测坐标系OXYZ，见图2，原点取为卫星中心Os，ssssX轴沿OO连线，离开地球方向为正，Z轴与X垂直指向正北，Y轴按右手系确定。由scssss于一般测量卫星的轨道都不会严格经过南北极上空，所以这种坐标系的定义是明确的。如此定义的观测坐标系也叫做UEN坐标系，因为三个坐标轴分别指向上（UP）、东（EAST）和北（NORTH）三个方向。根据变质量质点的动力学，空间飞行器在基础坐标系下的主动段的简化运动方程如下：

r(t)eTr

eT表示火箭产生的推力加速度，瞬时质量；m(t)是质量变化率；r(t)为空间飞行器在基础坐标系下的位置矢量；r(t)表示r(t)G

G3.9860051014m3/s为了更明确地表示推力加速度的方向，vm于火箭尾部喷口的喷射速度的逆矢量。方程（1）中如果只保留右侧第一项，则可以表示观测卫星的简化运动方程：

r(t)e

在给定基础坐标系下的位置和速度初值情况下，可以利用常微分方程组数值解方法计

算空间飞行器的运动轨迹。不同空间飞行器的本质差异就在于m(t)一般而言应为严格单调递减的非负函数。v

或相同，其大小一般较为稳定。观测卫星对于空间飞行器的观测数据通过化简可以由观测坐标系下的两个无量纲比值确定：

;

ss其中x,y,z为空间飞行器在观测坐标系中的坐标。sss

11

来表示，分别表示第一观测量的平移量、第二观测量的平移量以及观测量在z,x,来表示，分别表示第一观测量的平移量、第二观测量的平移量以及观测量在z,x,y,z。卫星编号从上到下递增并从0开观测卫星在任意时刻的位置计算是估计的前提，请根据satinfo.txt和观测卫星的在本题给定的仿真数据下，06号和09号观测卫星对0号空间飞行器形成了立体机误差为直接叠加在观测数据上的白噪声，可能产生于背景辐射干扰与信息处理等多个方面。系统误差也包括多种来源，如卫星定位误差、指向机构误差、图像校准误差、传感器安装误差等等。在本题框架内，我们假定只考虑与卫星平台相关的系统误差，即不同观测卫星的系统误差相互没有关联，同一观测卫星对于不同空间飞行器的系统误差是一样的。经由适当的简化模型，各种系统误差最终可以折合为观测坐标系的原点位置误差和三轴指向误差。根据工程经验，原点位置误差影响较小，而三轴指向误差影响较大，对三轴指向误差进行估计对于提高估计精度很有帮助，本题只考虑三轴指向误差。三轴指向误差在二维观测数据平面上表现为两个平移误差和一个旋转误差，具体可以用三个常值小量d,d,d

平面内的旋转量。

单个红外光学探测器不具备测距能力，但借助多颗（含两颗）观测卫星的同步观测能够进行逐点定位，再结合空间飞行器的运动模型，可以进行轨道参数估计。在单星观测条件下，利用空间飞行器轨道的特殊性，结合较强的模型约束也可得到一定精度轨道参数估计。由于受大气影响，垂直上升段的火箭尾焰不易观测，程序拐弯段的运动方程又较为复杂，所以本题重点关注重力斜飞段的后程段，本题所附仿真数据也集中于此段。本题以中低轨近圆轨道卫星为观测星座对假想的空间飞行器进行仿真观测，生成仿真观测数据，要求利用仿真观测数据，对假想空间飞行器的轨道参数进行估计。本题所附文件包括：参数文件satinfo.txt用来存储观测卫星信息，每行表示一颗卫星，包含六列，分别表示零时刻卫星在基础坐标系下的位置和速度x,y,

始。仿真数据文件meadata_i_j.txt用来存储仿真观测数据信息。i、j为占位符，表示编号为i的卫星对编号为j的飞行器的仿真观测数据信息，按照时间顺序分行，每行分三列，分别是观测时刻t以及对应观测数据,。

本题所涉及的数据与结果，均应采用国际标准单位，即：时间单位为秒、距离单位为米、速度单位为米每秒等等；所有位置和速度均指基础坐标系下的位置和速度。在仅考虑随机误差的条件下，请你们团队研究下列问题：1.简化运动方程（2），计算09号观测卫星在50.0s、100.0s、150.0s、200.0s、250.0s五个时刻的三维位置。结果保留6位有效数字。2.交叠观测，请结合立体几何知识按照逐点交汇定位的思路，给出0号空间飞行器

12

若06和若06和09号两颗观测卫星均有可能带有一定的系统误差，对系统误差进行正确对只有09号观测卫星单星观测的01号空间飞行器进行轨道估计，结果形式要求

从50.0s到170.0s间隔10.0s进行采样，计算并列表给出0号空间飞行器在各个采样点的位置和速度，并给出估计残差。结果保留6位有效数字。同时绘制0号空间飞行器的三个位置t-x、t-y、t-z和三个速度t-vx、t-vy、t-vz曲线示意图。在同时考虑系统误差的条件下，进一步研究下列问题：3.的估计能够有效提高精度。利用上述的逐点交汇方法能否同时对系统误差进行估计？若不能，是否还有其他的思路能够同时估计系统误差与轨道？给出你的解决方案与估计结果。在报告中除给出与第二问要求相同的结果外，还应分别给出两颗观测卫星的系统误差估计结果，共六个数值，分别是两颗卫星的d,d,d。

如果你们还有时间和兴趣，还可考虑下列：4.同第三问，注意参考第三问的系统误差估计结果。并进一步考虑在同时有多颗观测卫星观测多个空间飞行器的情况下能否联合进行系统误差估计？本题要求提供可计算出所提交报告中答案的计算程序，所使用的语言和工具不限，但推荐使用C\C++、Fortran、Matlab、Mathematica、...。

参考文献[1].王志刚,施志佳.远程火箭与卫星轨道力学基础[M].西北工业大学出版社，2006.[2].张毅,肖龙旭,王顺宏.弹道导弹弹道学[M].国防科技大学出版社,2005.[3].中国人民解放军总装备部军事训练教材编辑工作委员会著.外弹道测量数据处理[M].国防工业出版社,2002.[4].王正明,易东云著.测量数据建模与参数估计[M].国防科技大学出版社,1996.[5].科普托夫编著.弹道式导弹设计和试验[M].国防工业出版社.

有杆抽油系统的数学建模及诊断C

目前，开采原油广泛使用的是有杆抽油系统（垂直井，如图1）。电机旋转运动转化为抽油杆上下往返周期运动，带动设置在杆下端的泵的两个阀的相继开闭，从而将地下上千米深处蕴藏的原油抽到地面上来。钢制抽油杆由很多节连接而成，具有相同直径的归为同一级，级数从上到下按1,2„进行编号，可多达5级，从上端点到下端点可能长达上千米。描述抽油

13

ttttt=0t信息的通用方法是示功图：它是该点随时间而变化的荷载（合力，向下为正）数据作为纵坐标，以该点垂直方向上随时间而变化的位置相对于时刻该点位置的位移数据作为横坐标构成的图形。函数关系表现为位移-荷载关于时间的参数方程。一个冲程（冲程的说明见附录）中示功图是一条封闭的曲线。构成示功图的数据称为示功数据。抽油杆上端点称为悬点，图4示意了悬点E的运动过程。在一个冲程期间，仪器以一系列固定的时间间隔测得悬点E处的一系列位移数据和荷载数据，据此建立悬点E的示功图称为悬点示功图。附件1、2中的位移-荷载数据是某油田某井采油工作时采集的悬点处原始示功数据。“泵”是由柱塞、游动阀、固定阀、部分油管等几个部件构成的抽象概念（见图2），泵中柱塞处的示功图称为泵功图。因为受到诸多因素的影响，在同一时刻t，悬点处的受力（荷载）与柱塞的受力是不相同的；同样，在同一时刻t，悬点处的相对位移与柱塞的相对位移也不相同。因此悬点示功图与泵功图是不同的。图5给出了理论悬点示功图和理论泵功图。示功图包含了很多信息，其中就有有效冲程，泵的有效冲程是指泵中柱塞在一个运动周期内真正实现从出油口排油的那段冲程。工程上一般根据示功图形状与理论示功图进行对比来判断抽油机工作状态。通过悬点示功图可以初步诊断该井的工作状况，如产量、气体影响、阀门漏液、沙堵等等。要精确诊断油井的工作状况，最好采用泵功图。然而，泵在地下深处，使用仪器测试其示功数据实现困难大、成本高。因此，通过数学建模，把悬点示功图转化为杆上任意点的示功图（统称为地下示功图）并最终确定泵功图，以准确诊断该井的工作状况，是一个很有价值的实际问题。

请解决以下问题：问题一：光杆悬点运动规律电机旋转运动通过四连杆机构转变为抽油杆的垂直运动。假设驴头外轮廓线为部分圆弧、电机匀速运动，悬点E下只挂光杆（光杆下不接其它杆，不抽油，通常用来调试设备）。请按附录4给出四连杆各段尺寸，利用附件1的参数，求出悬点E的一个冲程的运动规律：位移函数、速度函数、加速度函数。并与有荷载的附件1的悬点位移数据进行比较。问题二：泵功图计算1966年，Gibbs给出了悬点示功图转化为地下示功图的模型[3],[4]，由于受

14

2u2uut2u2uutxt22悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；附件1是只有一级杆的某油井参数和悬点示功数据，附件2是有三级杆的另一油井参数和悬点示功数据，利用它们分别计算出这两口油井的泵功图数据；并分别绘制出两油井的悬点示功图和泵功图（每口井绘一张图，同一井的悬点示功图与泵功图绘在同一张图上，请标明坐标数据）。问题三：泵功图的应用（下面2小问选作一问。鼓励全做）1）建立2个不同的由泵功图估计油井产量的模型，其中至少一个要利用“有效冲程”；并利用附件1和附件2的数据分别估算两口油井一天（24小时）的产液量。（单位：吨，这里所指的液体是指从井里抽出来的混合液体）2）如图5（C）形式的泵功图表示泵内有气体，导致泵没充满。请建立模型或算法，以由计算机自动判别某泵功图数据是否属于泵内有气体的情况。并对附件1、附件2对应的泵功图进行计算机诊断是否属于泵内充气这种情况。问题四：深入研究的问题（下面2小问选作一问。鼓励全做）1）请对Gibbs模型进行原理分析，发现它的不足。在合理的假设下，重新建立抽油系统模型或对现有模型进行改进；并给出由悬点示功图转化为泵功图的详细计算过程，包括：原始数据的处理、边界条件、初始条件、求解算法；利用附件1、附件2的数据重新进行计算；对计算结果与问题二的计算结果进行比较，分析你的模型的优缺点。

2）Gibbs模型在数学上可简化为“波动方程”：a其中a为已知常数，c称为阻尼系数，鉴于大多数的阻尼系数公式[1][2]是作了诸多假设后推出的，并不能完整地反应实际情况。如果能从方程本身和某些数据出发用数学方法估计参数c，贡献是很大的。对此，请你进行研究，详细给出计算c的理论推导过程并尽可能求出c。如果需要题目之外的数据，请用字母表示之并给出计算c的推导过程。

参考文献[1].王鸿勋张琪，《采油工业原理》，石油行业出版社，1985年4月,第二章[2].万仁溥，《采油工程手册》，石油行业出版社，2000年8月，第五章第二节[3].Gibbs.S.G,Neely,A.B,ComputerDiagnosisofDownholeConditioninSuckerRodPumpingWells,J.Pet.Tech.,Jan.1966.[4].Gibbs.S.G,MethodofDeterminingSuckerRodPumpPerformance,UnitedStatesPatentOffice,Sep.1967.

15

附录1：有杆抽油系统实图及名词解释

图1.有杆抽油系统地面实图说明：前臂、后臂是对应的是同一个部件称为游梁，游梁与驴头固定连接；钢缆固定在驴头上部；悬点是钢缆与光杆的连接点；光杆是第一节抽油杆，长度比系统的理论冲程稍长。光杆上接钢缆，下接其它抽油杆，由于光杆有时与空气接触有时与油接触，环境较恶劣，所以材质较好，做得也比较光滑。光杆与第一级抽油杆粗细相同，计算时把光杆与第一级抽油杆同等看待，长度也计入了第一级。

16

出油口抽油杆

附录2：有杆抽油系统的工作原理出油口抽油杆

套管压力

油管压力

动液面地面

套管

柱塞油管游动阀

固定阀

油层

图2：有杆抽油系统原理图上冲程：抽油杆带着柱塞向上运动，柱塞上的游动阀受管内液柱压力而关闭。此时，柱塞下面的泵腔容积增大，压力降低，固定阀在其上下压差下打开，原油吸入泵中。如果油管内被液体充满，上冲程将在出油口排出相当于柱塞冲程长度的一段液体。原来作用在固定阀上的油管内的液柱压力将从油管转移到柱塞上，从而引起抽油杆柱的伸长和油管的缩短。上冲程是泵内吸入液体，而井口排出液体的过程。下冲程：抽油杆带着柱塞向下运动，柱塞压缩固定阀和游动阀之间的液体。当泵内压力增大到一定程度时，固定阀先关闭，当泵内压力增大到大于柱塞以上的液柱压力时，游动阀被顶开，柱塞下面的液体通过游动阀进入柱塞上部。原来作用在柱塞以上的液体重力转移到固定阀上，因此引起抽油杆柱的缩短和油管的伸长。柱塞向上向下活动一次叫一个冲程，冲程还表示物体在一个运动周期内的最大位移。每分钟完成冲程的次数成为冲次。

17

(b)(b)—柱塞上行接近上死点(d)—柱塞下行接近下死点O(c)(d)

附录3：泵的抽汲循环及阀门开闭(b)(b)—柱塞上行接近上死点(d)—柱塞下行接近下死点O(c)(d)

游动阀柱塞固定阀

(a)

(a)—柱塞接近下死点处上行(c)—柱塞接近上死点处下行

图3：泵的抽汲循环及阀开闭示意图附录4：悬点运动过程

A

B

E

D

φ

O’

图4：抽油机器四连杆机构简图

18

荷载C（B）理论泵功图形状AD（C）理论泵内充气泵功图形状B

“悬点”E的运动过程：t=0时刻，曲柄滑块D位于上顶点（=0），AB荷载C（B）理论泵功图形状AD（C）理论泵内充气泵功图形状B平行于水平面，E对应坐标原点（称为E的下死点），E的位移为0；D运动到下顶点（=）时，E的位移到达最大（称为E的上死点）；D接着运动到上顶点（=2）时，E又回到位移为0的位置，完成一个周期（即一个冲程）。前臂AO=4315mm，后臂BO=2495mm，连杆BD=3675mm，曲柄半径O’D=R=950mm

附录5：理论悬点示功图与理论泵功图

荷载

位移位移

（A）理论悬点示功图形状

图5：理论悬点示功图与理论泵功图

附录6：几点说明1、问题二得到的泵功图数据（位移，荷载）请在论文中以适当的方式表达；同时，按照附件所给数据格式，填入到“提交数据.xls”中，并将该文件名换为：“题号+报名号.xls”。2、本题所有抽油杆均为钢制，密度为：8456（单位：kg/m3），弹性模量为：2.1*1011（单位：pa）。3、因为悬点功图数据是自动测试的，附件1、2所给数据的第一对并不一定刚好是一个冲程的起点。上行和下行用的时间也不一定完全相等，请自行判断那些数据属于上冲程，那些数据属于下冲程。4、本题所有抽油机的油管是锚定的，因此本题不必考虑抽油管的长度变化（伸缩）。

19

D基于卫星云图的风矢场（云导风）度量模型与算法探讨

卫星云图在掌握大气环流、中长期天气预报以及灾害性天气学的

研究中有重要作用。它由地球同步卫星上的红外探测仪探测地球上空

的温度数据再转换成灰度数据制作而成。附件中定标数据文件

k_temp.txt给出了灰度数据与温度数据的转换关系，k_temp.txt内有

1024个实型数，依次是图象灰度数据为0到1023所对应的K氏温度

值，灰度值为-1时对应的是地球以外的探测点。[注：地球是被探测

温度的唯一来源，如果天空无云，探测到的温度可以看成是地球表面

的温度；在有云层的地方，探测到的温度相对较低，且云层越高越厚

温度就越低，探测到的温度可看成云层所在区域的温度]。红外探测

仪扫描采样时，按步进角(南北方向)和行扫描角（东西方向）均为140

微弧（1弧度＝1000000微弧）采样。在卫星与地球中心的连线和地

球表面的交点（称为星下点）处的分辨率大约是5公里。本题提供的

卫星探测数据文件都是2288×2288的灰度值矩阵，矩阵的每个元素

都对应地球上或地球外的一个探测点（或称采样点）。同步卫星离地

球中心的高度为42164000米,星下点在东经86.5度,北纬0度，星

下点对应的矩阵元素位于矩阵的第1145行和第1145列相交处。

为解答本题，首先要确定灰度矩阵中每个元素对应的采样点在地

球上的经纬度。地球可视为理想椭球，这个理想椭球可以由地球的一

个经过南北极的椭圆截面绕南北极的连线旋转而得到，椭圆截面的长

半轴（赤道半径）=6378136.5m,短半轴（极半径）=6356751.8m；据

20

此就可以将灰度矩阵中非负元素的行列号按上北下南、左西右东的地

图规则换算成地球上经纬度坐标，此结果既可用于估算各探测点之间

的距离，还可用于在云图上依据海岸线经纬度坐标标出海岸线以方便

看图。

观测大气环流情况的一个方法是在卫星云图上标出风矢。风矢的

大小和方向由云块移动的速度决定。风矢与风的速度有所不同，如某

个台风中一些区域的风速可达每秒五、六十米，而台风（看作云块）

中心的移动速度可能仅每小时十多公里。没有云或云块不稳定处的风

矢规定为零风矢，这种用云块的移动所定义的风矢被称为云迹风。气

象部门已经有一些方法根据变化的卫星云图计算云迹风，这类方法称

为云导风方法。计算云迹风时通常将云块大小限定为16×16个像素，

搜索范围限定为64×64个像素。

本题的主要目的是希望大家充分利用卫星图像数据及其特点建

立尽可能准确地描述实际风矢场的度量模型和算法。

题目提供了我国风云2号卫星获得的三个灰度矩阵，

IR1_2030.mat,IR1_2100.mat,IR1_2130.mat，分别表示某天的20:30,

21:00,21:30时刻红外探测仪探测到的地球上空的温度数据对应的灰

度值。又给出了海岸线经纬度坐标数据文件coastline0.txt，文件的第

1列为经度(东经),第2列是纬度(北纬),每一行2个数据对应海岸线

上一点，而特大数据（99999.99，99999.99）表示前一曲线已结束，

将要开始下一曲线。

具体要求解决如下问题：

21

1.换算视场坐标。给出灰度矩阵元素行列号对应于经纬度坐标的

换算公式，建立矩阵形式的经纬度坐标文本文件，这里矩阵的第i行

与第j列，分别对应灰度矩阵的450+i行与450+j列，矩阵元素是（经

度，纬度）这种形式的二维数组，给出结果的范围为