提高低信噪比核磁共振測(cè)井弛豫信息真實(shí)性回波處理方法

王鵬，譚茂金，2

（1.中國地質(zhì)大學(xué)地球物理與信息技術(shù)學(xué)院，北京100083；2.地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室，北京100083）

0 引 言

核磁共振（NMR）測(cè)井回波信號(hào)（CPMG Echo Trains）的采集不可避免地要受到儀器設(shè)備和井下環(huán)境噪聲的影響，記錄到的回波信號(hào)帶有很強(qiáng)的噪聲，給反演解釋工作帶來極大的不確定性。當(dāng)信噪比小于5時(shí)，直接反演便不能得出可靠的T2譜分布。Hansen等［1］針對(duì)噪聲影響的強(qiáng)度給出了離散的Picard條件，說明當(dāng)信噪比低于5時(shí)，回波信號(hào)與弛豫譜滿足的第1類Fredholm積分方程所反映的模型已經(jīng)完全不滿足離散Picard條件，也就是說因?yàn)樵肼暤挠绊憜螐臄?shù)學(xué)角度不能反演出真實(shí)可靠信息，必須對(duì)數(shù)據(jù)濾波或者重新采樣［1－2］。

實(shí)際生產(chǎn)中，數(shù)據(jù)采集采用不同相位信號(hào)疊加，壓制了很大一部分噪聲，對(duì)于解釋處理，疊加后的數(shù)據(jù)仍具有相當(dāng)高的噪聲。低信噪比核磁共振數(shù)據(jù)處理有翁愛華等［3］提出在反演中將基線偏移作為獨(dú)立的參數(shù)參與反演，并結(jié)合時(shí)間依賴濾波技術(shù)提高了長弛豫信息反演的準(zhǔn)確性。該方法將噪聲影響合理地抽象成反演參數(shù)并參與反演與校正，但沒能在小孔隙信息恢復(fù)上得到有效校正。Ma S等［4］提出了改進(jìn)的分窗小波降噪方法，并與中值濾波、普通小波濾波和FIR濾波的降噪效果進(jìn)行了對(duì)比，得出改進(jìn)的分窗小波方法能夠更好地提高帶噪信號(hào)的信噪比，但在反演信息準(zhǔn)確恢復(fù)上缺乏進(jìn)一步驗(yàn)證。Ahmed O A W［5］提出了SLFT（Stable Linear Time－Frequence Transforms）方法，通過限制和平均變換系數(shù)實(shí)現(xiàn)了對(duì)分子藥物核磁共振數(shù)據(jù)成像上的降噪。對(duì)于低磁場(chǎng)高噪聲核磁共振測(cè)井?dāng)?shù)據(jù)，該方法目前僅具有借鑒意義。

本文針對(duì)低信噪比的核磁共振測(cè)井?dāng)?shù)據(jù)提出了基于奇異值分解法SVD（Singlar Value Decomposition）濾波的核磁共振數(shù)據(jù)變換反演方法BRD（Butler Reeds Dawson）。該方法利用Hankel矩陣與回波信號(hào)的相互構(gòu)建關(guān)系，把回波信號(hào)的Hankel矩陣作為SVD濾波處理的信號(hào)矩陣，實(shí)現(xiàn)了對(duì)回波信號(hào)的分解和濾波。結(jié)合BRD反演算法，利用SVD處理后的高信噪比全波列數(shù)據(jù)反演出流體的T2分布。在隨后的模擬數(shù)據(jù)試驗(yàn)和實(shí)際數(shù)據(jù)處理中，可以看到這樣的分解和濾波在確?？蓜?dòng)流體成分和提高束縛流體成分在反演信息準(zhǔn)確恢復(fù)上是有效的。

1 處理方法與原理

整個(gè)數(shù)據(jù)的處理流程包括回波信號(hào)分解為高信噪比部分和低信噪比部分；低信噪比數(shù)據(jù)濾波；低、高信噪比數(shù)據(jù)反演；反演譜疊加。

1.1 回波信號(hào)分解

設(shè)核磁共振采集的回波信號(hào)為X，對(duì)應(yīng)的T2譜為F。X可以分解為多項(xiàng)疊加的形式，如式（1）所示分解為X1和X2，相應(yīng)的T2譜為F1和F2，與F的關(guān)系如式（2）所示，有

式（1）稱為回波信號(hào)分解，這樣的分解處理意義在于，對(duì)于低信噪比回波信號(hào)，回波信號(hào)的有效分解能夠保證更多的有效信息與噪聲信號(hào)分離，提高主要信息在反演處理中恢復(fù)的準(zhǔn)確性。

圖1（a）給出了含有100個(gè)樣本的回波信號(hào)X，X分解為X1和X2的2個(gè)信號(hào)，它們之間滿足式（1）。圖1（b）給出了回波信號(hào)的T2譜，從左往右依次是X、X1和X2對(duì)應(yīng)T2分布，它們之間的關(guān)系滿足式（2）。

圖1 回波信號(hào)分解

1.2 SVD濾波處理

1.2.1 濾波原理

對(duì)于一個(gè)m×n維信號(hào)矩陣A，對(duì)其進(jìn)行SVD分解［6－7］，有

式中，U為m×m維正交矩陣；S為m×n維似對(duì)角矩陣；V為n×n維正交矩陣。當(dāng)m＝n，S＝Σ；m＜n，S＝［Σm×mOm×（n－m）］；m＞n，S＝［Σn×n O（m－n）×n］T。Om×（n－m）是m×（n－m）維 零 矩 陣，O（m－n）×n是（m－n）×n維零矩陣。Σ＝diag（σ1，σ2，…，σp，…，σr）為非負(fù)對(duì)角矩陣，它的對(duì)角元素是由從大到小排列的奇異值構(gòu)成。

通常，沒有噪聲的信號(hào)矩陣的秩很小，而含噪聲的信號(hào)矩陣是一個(gè)滿秩矩陣。受噪聲的影響，原來排列靠后很小的奇異值會(huì)變得很大，SVD濾波原理就是將排列在某一奇異值之后所有奇異值置0，通過重構(gòu)信號(hào)矩陣達(dá)到降噪的目的［8－9］。

為了描述置0奇異值分界點(diǎn)的選取，定義重構(gòu)信號(hào)矩陣的能量比為

式中，p作為置0奇異值的分界點(diǎn)，奇異值σp以后的奇異值將被置0，輸出信號(hào)占總能量比例為η。在實(shí)際處理中，為獲取最優(yōu)濾波信號(hào)，往往根據(jù)信號(hào)矩陣奇異值分布曲線的拐點(diǎn)選取p的位置，此時(shí)得到的能量比稱為最優(yōu)能量比。

1.2.2 信號(hào)矩陣

對(duì)于核磁共振回波數(shù)據(jù)X，它由N個(gè)等間隔采樣點(diǎn)構(gòu)成，表示為

考慮噪聲影響，令S表示期望信號(hào)，N0表示噪聲，S和N0與X有相同的采樣形式，式（5）可寫為

根據(jù)回波信號(hào)與Hankel矩陣的構(gòu)建關(guān)系［10］，令回波信號(hào)的Hankel矩陣作為反映該回波信息的矩陣，稱為信號(hào)矩陣。回波信號(hào)X的信號(hào)矩陣H如式（7）所示，信號(hào)S和噪聲N0也可以寫出自己的信號(hào)矩陣，它們之間的關(guān)系H＝HS＋HN；HS和HN是S和N0對(duì)應(yīng)的信號(hào)矩陣，有

信號(hào)矩陣大小的選取應(yīng)滿足min（R，L）＞rank（H），對(duì)于數(shù)據(jù)大樣本信號(hào)，構(gòu)建信號(hào)矩陣應(yīng)盡可能減少行或列的維數(shù)，以減少計(jì)算量。

信號(hào)矩陣通過式（8）計(jì)算可恢復(fù)相應(yīng)回波信號(hào)

式中，g＝max（1，i－R＋1）；k＝min（i，N－R＋1）。

1.2.3 算法實(shí)現(xiàn)

回波信號(hào)的SVD濾波處理包含3個(gè)步驟。①構(gòu)建信號(hào)矩陣，利用式（7）構(gòu)建待處理回波信號(hào)的信號(hào)矩陣；②重構(gòu)矩陣，根據(jù)濾波目的選擇合適的能量比，得到降噪后的信號(hào)矩陣；③重構(gòu)信號(hào)，由回波信號(hào)濾波后的矩陣經(jīng)過式（8）得到重構(gòu)信號(hào)，此時(shí)的重構(gòu)信號(hào)即為經(jīng)過SVD濾波處理后的回波信號(hào)。

1.3 反演弛豫信息

核磁共振測(cè)井依據(jù)地層流體在外加磁場(chǎng)作用下的弛豫特性反映地層流體孔隙度和區(qū)分地層中可動(dòng)流體與束縛流體。根據(jù)核磁共振測(cè)井原理，從核磁共振測(cè)井采集的CPMG回波信號(hào)反演T2分布是典型的病態(tài)問題［2，11］，該問題被描述為式（9）及考慮噪聲影響的離散形式（10），有

式中，xi為記錄的第i個(gè)回波幅度的回波信號(hào)；ti（其中i＝1，2，…，N）為第i個(gè)相等回波間隔的衰減時(shí)間；Tj對(duì)應(yīng)第j種流體的弛豫特性（其中j＝1，2，…，M），即弛豫時(shí)間T2；fj為第j種流體的T2幅度。TW為等待時(shí)間，R為T1／T2，方程（1）為第1類Fredholm方程。

測(cè)量信號(hào)xi不可避免地要受到噪聲的影響。小的噪聲干擾可能造成fj很大程度的不穩(wěn)定，使得T2的分布強(qiáng)度fj不能有效、唯一地確定。目前提出的很多反演算法和正則化方法旨在抑制解的不穩(wěn)定性和加強(qiáng)解的可靠性［12－15］，其中較為有效的算法是Venkataramanan等［13］提出BRD反演算法，該算法在目標(biāo)方程中增加一項(xiàng)范數(shù)比例懲罰項(xiàng)抑制噪聲影響。如目標(biāo)方程（11）所示，方程中增加了帶有平滑因子α的懲罰項(xiàng)

具體計(jì)算中，解的可行域范圍內(nèi)預(yù)設(shè)弛豫基確保解的合理性以及計(jì)算效率。例如，均勻分布在0.1～104ms范圍內(nèi)的30個(gè)離散值構(gòu)成弛豫基。為方便起見，方程（12）給出了矩陣形式，待求參數(shù)為弛豫時(shí)間T2的分布強(qiáng)度向量F

1.4 處理流程

核磁共振測(cè)井采集的回波信號(hào)受噪聲的影響主要表現(xiàn)，①在大、中孔隙可動(dòng)流體成分造成基線偏移，在反演的弛豫譜上基線偏離表現(xiàn)為異常的長弛豫組分，這種影響在核磁數(shù)據(jù)信噪比較低時(shí)尤其明顯；②在微、小孔隙束縛流體成分產(chǎn)生畸變，在反演的弛豫譜上表現(xiàn)為跳變的尖銳假峰。提高低信噪比回波信號(hào)反演弛豫信息的準(zhǔn)確性，就是確?？蓜?dòng)流體成分和提高束縛流體成分在反演信息準(zhǔn)確恢復(fù)。

首先，應(yīng)用SVD濾波高效的信噪分離能力進(jìn)行回波信號(hào)分解，得到高信噪比的回波信號(hào)X1和低信噪比的回波信號(hào)X2，高信噪比的回波信號(hào)X1含有大量的T2弛豫信息可以直接反演得到相應(yīng)的譜分布。其次，低信噪比的回波信號(hào)X2嚴(yán)重受到噪聲影響，需要再次經(jīng)過SVD濾波提取有用信息，然后反演得到相應(yīng)譜分布。最后，根據(jù)回波信號(hào)分解關(guān)系，疊加2部分反演的T2譜作為最終要求取的T2分布。

反演數(shù)據(jù)的處理過程中2次用到SVD濾波。

（1）反演的回波信號(hào)經(jīng)過SVD濾波實(shí)現(xiàn)回波信號(hào)的分解。SVD濾波算法采用低于最優(yōu)能量比的設(shè)置，保證分解得到高信噪比數(shù)據(jù)X1。

（2）由回波信號(hào)分解理論，由X2＝X－X1計(jì)算得到相應(yīng)的低信噪比回波信號(hào)X2。

（3）低信噪比數(shù)據(jù)X2的SVD濾波，提取有用信息?；夭ㄐ盘?hào)X2含有較高噪聲，SVD濾波采用高于最優(yōu)能量比的設(shè)置，保證有用信息不被丟失。

（4）回波信號(hào)反演。2部分經(jīng)SVD濾波后的數(shù)據(jù)，反演得到相應(yīng)T2分布、F2和F1。

（5）反演譜的疊加。由F＝F1＋F2計(jì)算回波信號(hào)X的T2分布、F。

2 方法驗(yàn)證與分析

根據(jù)鄧克俊等［11］提出的 FCD（Fluid Component Decomposition）回波串構(gòu)建弛豫信號(hào)的方法，取Gaussian函數(shù)作為基函數(shù)，高斯寬度（Full Width at Half Maximum）為0.9，考慮不同強(qiáng)度的噪聲影響，正演得到不同信噪比的回波信號(hào)。對(duì)于含有束縛水和中等黏度油的雙峰特征T2譜模型，束縛水的T2峰值在20ms處，峰值幅度為0.4，含水飽和度為40%；油的T2峰值在300ms處，峰值幅度0.6，含油飽和度60%；模擬采集間隔TE＝0.9ms，回波數(shù)N＝500；模擬得到的回波信號(hào)見圖2。

圖2 模擬信號(hào)處理結(jié)果對(duì)比

圖2中藍(lán)色實(shí)線為沒有加入噪聲的回波信號(hào)，即期望信號(hào)，黑色點(diǎn)線為帶噪聲的回波信號(hào)，紅色虛線為經(jīng)過SVD濾波處理反演前2部分疊加的回波信號(hào)。圖2（a）到圖2（e）分別為信噪比50、30、20、10和5的帶噪回波信號(hào)，相應(yīng)的經(jīng)過SVD濾波處理后的重構(gòu)信號(hào)與期望信號(hào)。從中可以看出濾波后的信號(hào)有效地恢復(fù)到了期望信號(hào)，但隨著信噪比下降過低時(shí)，SVD濾波處理仍舊保證了大弛豫成分信息的準(zhǔn)確性，小弛豫信息會(huì)造成少量的丟失。

同時(shí)，反演了相應(yīng)的弛豫譜。反演計(jì)算中，弛豫基分布區(qū)間0.3～3 000ms，布點(diǎn)數(shù)30（見圖3），藍(lán)色實(shí)線表示弛豫模型，黑色點(diǎn)實(shí)線表示帶噪回波信號(hào)直接反演得到的弛豫譜（Fn），紅色星實(shí)線表示經(jīng)本文濾波處理得到的弛豫譜（F）。

如圖3所示，應(yīng)用本文方法處理得到的T2譜更接近模型，反演譜峰的位置與模型近乎一致。隨著信噪比的降低，小弛豫成分峰值會(huì)有所丟失，但還是保證了弛豫位置的準(zhǔn)確性；而直接反演得到T2譜只在油峰位置與模型相近，隨著信噪比的降低，小弛豫位置會(huì)出現(xiàn)假峰。從對(duì)比的結(jié)果看，本文方法在低信噪比數(shù)據(jù)的處理上更接近真實(shí)情況。

圖3 模擬信號(hào)處理結(jié)果對(duì)比

為定量說明反演效果，令反演得到的T2譜用列向量f表示，T2譜模型用列向量m表示，定義相對(duì)誤差Er作為反演精度的度量，有

根據(jù)式（12），相應(yīng)反演譜的相對(duì)誤差見表1。由表1可知，隨著回波信號(hào)信噪比的降低，直接反演方法得到譜分布的相對(duì)誤差增加的很快，說明了當(dāng)信噪比低或過低時(shí)，直接反演結(jié)果的不可靠性；與直接反演方法相比，在較低信噪比情況下，采用本文濾波反演方法得到譜的相對(duì)誤差更小，誤差變化幅度也更小。

當(dāng)回波信號(hào)信噪比較高（SNR≥20）時(shí)，直接反演譜和經(jīng)濾波反演譜的流體弛豫位置都與模型相一致，但濾波反演譜的整體相對(duì)誤差要略大于直接反演結(jié)果，說明較高信噪比回波信號(hào)濾波處理時(shí)，由于信噪的分離難度加大，在設(shè)置濾除低強(qiáng)度噪聲的同時(shí)，部分有用信息也被濾除了。因此，較高信噪比回波信號(hào)處理時(shí)應(yīng)采用設(shè)置為全通能量的濾波反演模式或跳過濾波處理直接反演。

表1 反演譜相對(duì)誤差

3 油田實(shí)例分析

某區(qū)×井核磁共振測(cè)井?dāng)?shù)據(jù)，測(cè)量儀器為MRIL－Prime，采用雙TW測(cè)量模式識(shí)別可動(dòng)流體，采集參數(shù)為，A組采用長磁化參數(shù)TW＝12.7s，B組采用短磁化參數(shù)TW＝2.0s，2組采用相同的采集間隔TE＝0.9ms和回波數(shù)NE＝500。2組回波信號(hào)的統(tǒng)計(jì)信噪比為，A組大約在2～19，B組大約在1～19，屬于較低的信噪比數(shù)據(jù)。

圖4為本文方法處理結(jié)果與商業(yè)軟件DPP處理結(jié)果的對(duì)比圖。圖4中，第1道GR、SP；第2道為地層電阻率Rt、沖洗帶電阻率Rxo；第3道為深度道；第4、5道為長TW，L和短TW，S回波信號(hào)的 DPP反演結(jié)果；第6道為DPP軟件處理的長TW，L與短TW，S弛豫差譜；第7、8道為本文濾波反演方法處理的長TW，L和短TW，S結(jié)果；第9道為本文濾波反演方法處理的長TW，L與短TW，S弛豫差譜?？梢钥闯?，在××55～××73m深度段與DPP處理結(jié)果相比，在小弛豫位置弛豫譜沒有出現(xiàn)噪聲引起的跳變尖銳峰；在大弛豫位置壓制了噪聲引起的基線偏移造成的大弛豫成分；在A組和B組的弛豫譜上能夠看到該層段存在大量的束縛流體和可動(dòng)流體；從第6道和第9道的差譜信息上看，兩者處理的差譜信號(hào)都很弱，但文中方法處理結(jié)果仍能看出可動(dòng)流體中含有一定量的烴。

圖4 實(shí)測(cè)數(shù)據(jù)處理結(jié)果對(duì)比

圖5 巖性柱狀圖

圖5為該層段錄井巖性及含油級(jí)別測(cè)試結(jié)果。錄井巖性顯示，該層段巖性為深灰、灰黑色泥巖、灰色粉砂質(zhì)泥巖與淺灰、灰褐色細(xì)砂巖呈略等厚互層，含油測(cè)試結(jié)果為油跡、油斑級(jí)別。該層段的測(cè)井響應(yīng)資料顯示，自然伽馬曲線值為64.7～89.7API，均值為78.0API。三孔隙度曲線變化平緩，聲波時(shí)差、密度均值分別為244.8μs／m、2.41g／cm3，中子16.8%～27.9%，說明地層孔隙發(fā)育，物性好。該層電阻率值略有變化，深感應(yīng)電阻率值為5.0～15.0 Ω·m，均值為6.8Ω·m，陣列感應(yīng)電阻率曲線明顯分異，120in＊非法定計(jì)量單位，1ft＝12in＝0.304 8m，下同陣列感應(yīng)電阻率值為4.9～14.9 Ω·m，孔隙度為13.8%，滲透率為0.62×10－3μm2，表明地層含油性較好，該層解釋為油水同層。

4 結(jié)論與認(rèn)識(shí)

（1）受噪聲影響，核磁共振測(cè)井采集的回波信號(hào)表現(xiàn)為在大、中孔隙可動(dòng)流體成分造成基線偏移，在微、小孔隙束縛流體成分產(chǎn)生畸變；在反演的弛豫譜上表現(xiàn)為異常的長弛豫和小弛豫位置跳變的尖銳假峰。尤其當(dāng)信噪比低于5時(shí)，單從數(shù)學(xué)角度不能反演出真實(shí)可靠信息，必須對(duì)數(shù)據(jù)濾波或重新采樣。

（2）提出基于SVD濾波的BRD反演處理方法，在低信噪比（SNR≤20）保證了可動(dòng)流體和束縛流體信息在反演處理中的有效恢復(fù)，提高低信噪比核磁共振數(shù)據(jù)反演弛豫信息的準(zhǔn)確性，在數(shù)值模擬和油田實(shí)例中都取得了令人滿意的效果。從油田實(shí)例的應(yīng)用效果看，與DPP處理結(jié)果相比，在小弛豫位置弛豫譜沒有出現(xiàn)噪聲引起的跳變尖銳峰，在大弛豫位置壓制了噪聲引起的基線偏移造成的大弛豫成分，在弛豫譜形態(tài)上能更精細(xì)地識(shí)別束縛流體和可動(dòng)流體。

（3）SVD濾波處理利用Hankel矩陣與回波信號(hào)的相互構(gòu)建關(guān)系，根據(jù)噪聲強(qiáng)度選取適當(dāng)?shù)木仃嚹芰勘葔褐坪徒財(cái)嘁蛟肼曇鸬钠娈愔档淖兓哂泻軓?qiáng)的信噪分離能力。

（4）實(shí)現(xiàn)了對(duì)回波信號(hào)的分解和濾波。在處理過程中，回波信號(hào)分解采用低于最優(yōu)能量比的SVD濾波設(shè)置，保證信號(hào)分解的有效性和回波信號(hào)的反演精度，回波信號(hào)濾波采用高于最優(yōu)能量比的SVD濾波設(shè)置，盡可能保留有效信息使之能夠在反演中恢復(fù)。

［1］ Hansen P C.The Discrete Picard Condition for Discrete Ill－Posed Problems［J］.BIT Numerical Mathematics，1990，30（4）：658－672.

［2］ 肖立志，謝然紅，廖廣志.中國復(fù)雜油氣藏核磁共振測(cè)井理論與方法 ［M］.北京：科學(xué)出版社，2012.

［3］ 翁愛華，李舟波，莫修文，等.低信噪比核磁共振測(cè)井資料的處理技術(shù) ［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2003，33（2）：232－235.

［4］ MA S，KONG L，CHEN J.An Improved NMR Signal De－noising Algorithm Based on Wavelet Transform［J］.Computational Information Systems，2011，7（13）：4651－4659.

［5］ Ahmed O A W.Method for Removing Noise from Nuclear Magnetic Resonance Signals and Images：U.S.Patent 7，253，627［P］.2007－08－07.

［6］ Kalman D.A Singularly Valuable Decomposition：The SVD of a Matrix［C］∥College Math Journal，1996.

［7］ Aharon M，Elad M，Bruckstein A.SVD：An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation［C］∥Signal Processing，IEEE Transactions，2006，54（11）：4311－4322.

［8］ Di F，Changzhi L，Qinguang C.SVD Filter Based on Noise Singular Values Clustering ［C］∥Intelligent Computation Technology and Automation （ICICTA），2010International Conference，IEEE，2010，3：665－668.

［9］ Sanliturk K Y，Cakar O.Noise Elimination from Measured Frequency Response Functions ［J］.Mechanical Systems and Signal Processing，2005，19（3）：615－631.

［10］ 劉志鵬，趙偉，陳小宏，等.局部頻率域SVD壓制隨機(jī)噪聲方法 ［J］.石油地球物理勘探，2012，47（2）：202－206.

［11］ 鄧克俊，謝然紅.核磁共振測(cè)井理論及應(yīng)用 ［M］.東營：中國石油大學(xué)出版社，2010.

［12］ Butler J P，Reeds J A，Dawson S V.Estimating Solutions of the First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing［J］.SIAM Journal on Numerical Analysis，1981，18（3）：381－397.

［13］ Venkataramanan L，Song Y Q，Hurlimann M D.Solving Fredholm Integrals of the First Kind with Tensor Product Structure in 2and 2.5Dimensions［C］∥Signal Processing，IEEE Transactions，2002，50（5）：1017－1026.

［14］ LIAO G Z，XIAO L Z，XIE R H，et al.Influence Factors of Multi－exponential Inversion of NMR Relaxation Measurement in Porous Media ［J］.Chinese Journal of Geophysics，2007，50（3）：796－802.

［15］ Wang W，Li P，Ye C.Multi－exponential Inversions of Nuclear Magnetic Resonance Relaxation Signal［J］.Science in China Series A：Mathematics，2001，44（11）：1477－1484.