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列，矩陣元素是（經

度，緯度）這種形式的二維數組，給出結果的范圍為