基于優(yōu)化MRVA的Φ-OTDR相干衰落抑制方法

丁文超 何海軍 鄭狄 鄒喜華

摘要：本文針對相位敏感光時(shí)域反射（Φ-OTDR）系統(tǒng)的一種相干衰落抑制方法旋轉(zhuǎn)矢量移動(dòng)平均（MRVA）進(jìn)行了優(yōu)化，獲得了性能更優(yōu)的相干衰落抑制. 該優(yōu)化方法從兩個(gè)方面展開，首先使用加權(quán)窗口代替MRVA方法中的矩形窗口，兼顧測量精度和空間分辨率；然后根據(jù)移動(dòng)平均與低通濾波的數(shù)理等效性，將時(shí)域移動(dòng)平均轉(zhuǎn)換為運(yùn)算速度更快的頻域低通濾波，提高信號處理效率.實(shí)驗(yàn)結(jié)果表明，相較于MRVA方法，該方法的解調(diào)相位精度提高了約23%，計(jì)算時(shí)間從0.957 s縮短至0.026 s，運(yùn)算速度提高約36倍.

關(guān)鍵詞：相位敏感光時(shí)域反射； 相干衰落抑制； 移動(dòng)平均； 頻域加速

中圖分類號：??TP212? 文獻(xiàn)標(biāo)識(shí)碼：A? DOI：10.19907/j.0490-6756.2023.053004

收稿日期： ?2023-08-27

基金項(xiàng)目： ??國家自然科學(xué)基金（62071395）； 四川省自然科學(xué)基金面上項(xiàng)目（2023NSFSC0493）； 中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金（2682023CG008）

作者簡介： ??丁文超（1999-）， 男， 內(nèi)蒙古赤峰人， 碩士研究生， 主要研究方向?yàn)楣饫w傳感. E-mail： Dinghan1005@163.com

通訊作者： ?鄒喜華. E-mail： zouxihua@126.com

A Φ-OTDR interference fading suppression method based on optimized MRVA

DING Wen-Chao， HE Hai-Jun， ZHENG Di， ZOU Xi-Hua

（Center of Information Photonics & Communications， School of Information Science and Technology，

Southwest Jiaotong University， Chengdu 611756， China ）

In this paper， an optimized rotating vector moving average （MRVA） method is proposed to suppress interference fading in phase-sensitive optical time-domain reflectometry （Φ-OTDR） system. The proposed optimization method is developed from two aspects. Firstly， the weighted window is used to replace the rectangular window in the MRVA method， which mitigates the tradeoff between measurement accuracy and spatial resolution. Then， the time domain moving average is converted to frequency-domain computation according to the mathematical equivalence of moving average and low-pass filtering， thus achieving faster processing speed. Experimental results show that， compared with MRVA method， the demodulation phase accuracy of this method is improved by about 23%， and the calculation time is reduced from 0.957 s to 0.026 s， achieving about 36-times improvement.

Phase-sensitive optical time domain reflectometry; Interference fading suppression; Moving average; Frequency domain acceleration

1 引 言

相位敏感光時(shí)域反射儀（Φ-OTDR）作為分布式光纖傳感技術(shù)的重要分支之一，已廣泛應(yīng)用于結(jié)構(gòu)健康監(jiān)測、水下聲學(xué)、地質(zhì)災(zāi)害監(jiān)測、周界安全防護(hù)等領(lǐng)域 ?［1-6］ . 通過提取Φ-OTDR信號的相位，不僅能準(zhǔn)確定位擾動(dòng)，還能準(zhǔn)確地重建外部振動(dòng)信號. 然而，在Φ-OTDR系統(tǒng)中使用高相干光源不可避免地會(huì)導(dǎo)致相干衰落. 當(dāng)在衰落區(qū)域中解調(diào)時(shí)，測量信號的低信噪比（SNR）不足以有效地解調(diào)相位，導(dǎo)致在某些時(shí)間點(diǎn)發(fā)生相位跳變. 相位跳變使得相位跡線上出現(xiàn)假警報(bào)峰值 ?［7］ ，導(dǎo)致無法提取干涉區(qū)域的相位信息，解調(diào)相位的精度和可靠性嚴(yán)重下降.

為了抑制相干衰落，提高解調(diào)相位信號的精度，相關(guān)研究人員提出了很多方法. Wang等人 ?［8］ 提出了一種相移雙脈沖技術(shù)，兩個(gè)脈沖的散射信號的衰落區(qū)域可以相互補(bǔ)償，消除了干擾衰落帶來的不利影響，但該方法犧牲了一半的響應(yīng)帶寬. 使用不同頻率的探測脈沖時(shí)，瑞利后向散射（RBS）跡線之間的相關(guān)性較弱，因此人們提出了多頻脈沖方法來降低發(fā)生相干衰落的概率 ?［9-15］ . 受多頻脈沖方法的啟發(fā)，空分復(fù)用方法抑制衰落噪聲的思路被提出. Zhang等人 ?［16］ 在Φ-OTDR中使用多芯光纖構(gòu)建空分復(fù)用，通過聯(lián)合解調(diào)沿多芯光纖不同芯傳輸?shù)亩鄺l獨(dú)立RBS跡線的相位有效抑制了衰落噪聲. 此外，通過弱反射器陣列 ?［17-19］ 來抑制干擾衰落的技術(shù)也受到了廣泛的關(guān)注. 然而，上述方案需要專用的探測信號或?qū)ｉT設(shè)計(jì)的光纖，這增加了系統(tǒng)的復(fù)雜性，昂貴的精密設(shè)備也增加了系統(tǒng)成本.

另一方面，數(shù)字信號處理是抑制相位解調(diào)Φ-OTDR中相干衰落的替代手段， 這種方法在數(shù)字域進(jìn)行信號處理，無需增加額外的器件和結(jié)構(gòu). 2019年，Wu等人 ?［20］ 采用了頻譜提取和混疊（SERM）方法來消除干擾衰落. SERM方法的優(yōu)點(diǎn)是在數(shù)字域?qū)π盘栠M(jìn)行處理. 然而，在混疊前需要用窄帶FIR濾波器分別提取拍頻信號的主瓣和兩個(gè)一階副瓣，這與多頻脈沖方法相比并沒有減少計(jì)算量，且為了不犧牲較多空間分辨率，僅僅將頻譜分成三部分，難以確保衰落抑制干凈. 2021年，He等人 ?［21］ 提出相移變換來抑制Φ-OTDR中的干擾衰落. 該方法需要輔助偏振分極探測，且為了完全抑制干擾衰落，需要進(jìn)行多次相移變換，所以計(jì)算時(shí)間過長，無法實(shí)現(xiàn)實(shí)時(shí)數(shù)據(jù)處理. 2022年，Qian等人 ?［22］ 提出了一種移動(dòng)旋轉(zhuǎn)矢量平均（MRVA）方法抑制干擾衰落. 在對RBS信號進(jìn)行旋轉(zhuǎn)矢量對齊后，選擇合適的窗口進(jìn)行移動(dòng)平均. 為了保持較高的空間分辨率，移動(dòng)平均的窗口受到了限制，對應(yīng)的相干衰落抑制效果和解調(diào)性能也受到了制約.

本文對MRVA方法進(jìn)行了優(yōu)化，將加權(quán)移動(dòng)窗口引入到相位解調(diào)Φ-OTDR系統(tǒng)中，在消除干擾衰落的同時(shí)，提高解調(diào)相位的精度. 同時(shí)，為了加快解調(diào)計(jì)算速度，將時(shí)域移動(dòng)平均轉(zhuǎn)換為運(yùn)算更快的頻域低通濾波. 實(shí)驗(yàn)結(jié)果表明，相較于傳統(tǒng)MRVA方法，該方法可以有效提高解調(diào)相位的精度和解調(diào)頻率的信噪比，同時(shí)顯著縮短數(shù)據(jù)處理時(shí)間.

2 原 理

2.1 加權(quán)移動(dòng)窗口

傳統(tǒng)MRVA方法采用矩形窗口進(jìn)行移動(dòng)平均，窗口越長，解調(diào)精度越高，但是空間分辨率劣化也更嚴(yán)重. 基于矩形窗的移動(dòng)平均，窗口內(nèi)所有數(shù)據(jù)對運(yùn)算結(jié)果給予相同的貢獻(xiàn)，因此空間分辨率劣化嚴(yán)重. 加權(quán)平均通常能夠有效地兼顧空間分辨率和解調(diào)精度 ?［23］ .基于此，本文對傳統(tǒng)移動(dòng)平均窗口進(jìn)行優(yōu)化，設(shè)計(jì)窗口中心的權(quán)重大于窗口邊緣，使得距離窗口中心數(shù)據(jù)對計(jì)算結(jié)果的貢獻(xiàn)大于窗口邊緣部分 ?［23］ ，如此便可以提高有效空間分辨率. 常規(guī)的加權(quán)窗口有高斯窗、三角窗、布萊克曼窗和漢寧窗等. 本文選擇了其中3種典型的加權(quán)窗口，如圖1所示.

圖1a為高斯窗，窗函數(shù)為

w ??Gaussian ??n =e ?- 1 2 ??α n （L-1）/2 ???2 ??（1）

其中 L 為窗口長度； n 為窗口內(nèi)各點(diǎn)的權(quán)重； α 為寬度因子，本文設(shè)α=1；圖1b為三角窗，窗函數(shù)為

w ??Triangular ??n = ??2n-1 ?L ， ?1≤n≤ L 2 ?2- ?2n-1 ?L ， ??L 2 +1≤n≤L ??（2）

其中 L 是窗口長度. 圖1c是布萊克曼窗，窗函數(shù)為

w ??Blackman ??n =0.42-0.5 ?cos ??2πn L-1 ?+ ??0.08 ?cos ??4πn L-1 ???（3）

由于加權(quán)值的分布特性，與矩形窗相比，加權(quán)窗口能夠獲得與目標(biāo)中心點(diǎn)更接近的均值，獲得良好平滑效果的同時(shí)兼顧空間分辨率，整體性能更加均衡 ?［24］ . 盡管加權(quán)窗口獲得與矩形窗相同空間分辨率的窗長更長，但可以獲得更高的解調(diào)精度，同時(shí)有效空間分辨率高于或等于矩形窗.

2.2 頻域加速

MRVA方法中移動(dòng)平均部分的計(jì)算量與窗口長度相關(guān)，移動(dòng)窗口越長，所需的計(jì)算量越大，耗費(fèi)的時(shí)間也越長. 如果從頻域的角度分析移動(dòng)平均方法，以矩形窗口的移動(dòng)平均為例，對于一個(gè)長度為 n 的序列，移動(dòng)平均的計(jì)算結(jié)果為

y n = 1 L ∑ ?L-1 ??k=0 ?x n-k ??（4）

式中， L 為正整數(shù). 輸入信號為沖擊函數(shù)

δ（n）= 1， n=0 0， n≠0 ??（5）

則移動(dòng)平均的輸出為

y n = 1 L ∑ ?L-1 ??k=0 ?δ n-k ??（6）

輸出結(jié)果對應(yīng)的頻域表達(dá)式為 ?［25］

H e ??jω ?= 1 L ∑ ?L-1 ??k=0 ?e ??-jωk = ??sin ?ωL/2 ?L sin ?ω/2 ??e ??-jω（L-1）/2 ??（7）

根據(jù)式（4）～式（7）可知，移動(dòng)平均本質(zhì)上就是低通濾波. 并且窗長越長，對應(yīng)的濾波帶寬越窄. 對圖1中的矩形窗、高斯窗、三角窗和布萊克曼窗函數(shù)分別做傅里葉變換，得到其幅頻響應(yīng)如圖2所示. 考慮到移動(dòng)平均方法時(shí)域運(yùn)算耗時(shí)較長，根據(jù)滑動(dòng)窗口的時(shí)頻對應(yīng)關(guān)系，可以將耗時(shí)較長的運(yùn)算轉(zhuǎn)換到頻域執(zhí)行，進(jìn)而實(shí)現(xiàn)運(yùn)算速度的提高.

綜合以上分析，本文從兩個(gè)方面對傳統(tǒng)MRVA方法進(jìn)行優(yōu)化. 一方面，使用不同類型的加權(quán)移動(dòng)窗口來替代原有的矩形窗，突破矩形窗解調(diào)精度和空間分辨率之間的制約；另一方面，以運(yùn)算更快速的頻域?yàn)V波代替耗時(shí)較長的時(shí)域滑動(dòng)平均，提高解調(diào)速度. 具體操作如下： （1） 構(gòu)造距離-時(shí)間信號. 將采集到的瑞利后向散射信號沿距離和時(shí)間重組為二維信號. （2） 旋轉(zhuǎn)矢量對齊. 沿距離軸進(jìn)行希爾伯特變換，對得到的復(fù)數(shù)向量進(jìn)行旋轉(zhuǎn)矢量對齊. （3） 設(shè)計(jì)移動(dòng)窗口. 根據(jù)探測脈沖寬度選擇合適長度的加權(quán)移動(dòng)窗口. （4） 頻域加速計(jì)算. 通過傅里葉變換得到該移動(dòng)窗口的幅頻特性曲線，與旋轉(zhuǎn)矢量對齊的信號在頻域相乘，再進(jìn)行傅里葉反變換得到頻域處理后的信號，相干衰落引起的相位假峰被消除. （5） 解調(diào)信號. 通過對處理后的數(shù)據(jù)進(jìn)行解調(diào)，可以準(zhǔn)確地得到擾動(dòng)的位置和擾動(dòng)波形.

3 實(shí)驗(yàn)與分析

本文搭建的相干探測Φ-OTDR實(shí)驗(yàn)系統(tǒng)如圖3所示. 系統(tǒng)采用線寬為100 Hz的窄線寬激光器（NLL）作為光源，其產(chǎn)生的光波經(jīng)過一個(gè)10∶90光耦合器（OC1）分成上下兩部分.上路的光（90%）的光經(jīng)由任意函數(shù)發(fā)生器（AFG）控制的聲光調(diào)制器（AOM，頻移為200 MHz）調(diào)制成為脈沖寬度 20 ns ，周期100 μs的光脈沖，該脈沖經(jīng)摻鉺光纖放大器（EDFA）放大后，通過光環(huán)行器（CIR）進(jìn)入9.6 km長的待測光纖（FUT）. 下路的光（10%）經(jīng)由偏振控制器（PC）調(diào)節(jié)偏振狀態(tài)后作為本振光與背向瑞利散射光進(jìn)行相干接收. 本振光和后向瑞利散射光通過50∶50的光耦合器（OC2）相互拍頻，然后通過平衡探測器（BPD）轉(zhuǎn)為電信號，并使用采樣率為1 GSa/s的采樣率進(jìn)行采集. 兩個(gè)壓電陶瓷振蕩器（PZT）作為振動(dòng)驅(qū)動(dòng)器對光纖施加振動(dòng)，第一個(gè)PZT在光纖7.8 km處對約1.2 m的裸纖施加500 Hz的正弦振動(dòng)，第二個(gè)PZT在光纖8.6 km處對約1.5 m的裸纖施加800 Hz的正弦振動(dòng).

首先使用MRVA方法對采集的數(shù)據(jù)進(jìn)行處理，差分相位的標(biāo)準(zhǔn)差（SD）可以準(zhǔn)確地反映相位波動(dòng)的分布，因此它可以用來評價(jià)解調(diào)相位的精度，而SD曲線的上升沿的寬度可以反映解調(diào)相位的空間分辨率. 將SD曲線上升部分以0.1 m的距離進(jìn)行差分，差分值大于0.1的區(qū)域定義為上升沿. 傳統(tǒng)MRVA方法中矩形移動(dòng)窗口的寬度取值應(yīng)為脈沖寬度對應(yīng)距離的2～3倍 ?［22］ ，在本實(shí)驗(yàn)中取2倍，即4.0 m. 以此寬度下矩形移動(dòng)窗口的SD曲線上升沿寬度為基準(zhǔn)，調(diào)節(jié)各加權(quán)移動(dòng)窗口的寬度，使得它們的SD曲線上升沿寬度相等，保證各種窗口的解調(diào)空間分辨率一致，從而可以比較不同類型窗口的移動(dòng)平均方法解調(diào)相位的精度，如圖4b和4c所示. 通過此方法得到各移動(dòng)窗口的寬度，將MRVA方法的矩形移動(dòng)窗口（窗長為4.0 m）分別替換為寬度因子 α =1的高斯型窗（窗長為7.5 m）、三角窗（窗長為6.4 m）和布萊克曼窗（窗長為8.1 m）后，對數(shù)據(jù)進(jìn)行移動(dòng)平均處理，解調(diào)的相位以4.0 m的差分距離前后做差，計(jì)算差分相位的標(biāo)準(zhǔn)差如圖4a所示. 經(jīng)計(jì)算，矩形窗、高斯窗、三角窗和布萊克曼窗移動(dòng)平均作用下的解調(diào)信號相位SD分別為0.0542、0.0417、0.0468和0.0494. 在同樣的空間分辨率下，高斯型窗口可以得到最高的解調(diào)精度.

圖5比較了四種移動(dòng)窗口處理下的解調(diào)信號的時(shí)域特性和頻域特性. 矩形窗、布萊克曼窗、三角窗和高斯窗移動(dòng)平均作用下的解調(diào)頻率信噪比分別為18 dB、25 dB、28 dB和30 dB，三種加權(quán)移動(dòng)窗口的信噪比均高于矩形移動(dòng)窗口.

移動(dòng)平均方法的計(jì)算時(shí)間取決于窗口長度，采用加權(quán)移動(dòng)窗口提高解調(diào)精度，增大了窗口長度，所以也增加了計(jì)算時(shí)間，這不利于信號的實(shí)時(shí)處理. 經(jīng)過第二節(jié)的分析可知，可以在頻域處理相同的數(shù)據(jù)，加快計(jì)算速度. 表1比較了四種移動(dòng)窗口的性能，使用的計(jì)算機(jī)處理器為11th Gen Intel（R） Core（TM） i7-11800H. 表中的所有結(jié)果都是在達(dá)到相同的空間分辨率的條件下取得的，可以看到，加權(quán)移動(dòng)窗口的解調(diào)相位精度和解調(diào)頻率的信噪比都優(yōu)于矩形窗口，但是時(shí)間消耗更長；采用頻域方法進(jìn)行加速過后，性能保持不變的同時(shí)顯著提高了計(jì)算速度.

4 結(jié) 論

本文對抑制Φ-OTDR系統(tǒng)中干擾衰落的傳統(tǒng)MRVA方法進(jìn)行了優(yōu)化，提出了使用加權(quán)移動(dòng)窗口代替MRVA中矩形移動(dòng)窗口的方法，并根據(jù)移動(dòng)窗口的時(shí)頻對應(yīng)關(guān)系將計(jì)算轉(zhuǎn)換到頻域執(zhí)行，大大提高了計(jì)算速度. 通過對比分析，該方法可以在保持空間分辨率不變的基礎(chǔ)上提高解調(diào)性能，有效抑制干擾衰落，準(zhǔn)確定位和恢復(fù)干擾信號.為了驗(yàn)證該方法的可行性，對兩點(diǎn)擾動(dòng)進(jìn)行了測試. 實(shí)驗(yàn)結(jié)果表明，該方法能夠正確地定位和恢復(fù)干擾信號.與MRVA方法相比，該方法在解調(diào)相位精度、解調(diào)頻率信噪比和處理時(shí)間等方面都得到了提高.本文為數(shù)字域抑制Φ-OTDR系統(tǒng)的相干衰落提供了一種新的方法.

參考文獻(xiàn)：

［1］ ??Juarez J C， Maier E W， Choi K N， ?et al . Distributed fiber-optic intrusion sensor system ［J］. J Lightwave Technol， 2005， 23： 2081.

［2］ ?Pan W X， Zhao K， Xie C， ?et al . Distributed online monitoring method and application of cable partial discharge based on φ-OTDR ［J］. IEEE Access， 2019， 7： 1.

［3］ ?Bao X Y， Zhou D P， Baker C， ?et al . Recent development in the distributed fiber optic acoustic and ultrasonic detection ［J］. J Lightwave Technol， 2017， 35： 3256.

［4］ ?Wu H J， Liu X R， Xiao Y， ?et al . A dynamic time sequence recognition and knowledge mining method based on the Hidden Markov Models （HMMs） for pipeline safety monitoring with Φ-OTDR ［J］. J Lightwave Technol， 2019， 37： 4991.

［5］ ?Fernández-Ruiz M R， Soto M A， Williams E F， ?et al . Distributed acoustic sensing for seismic activity monitoring ［J］. APL Photonics， 2020， 5： 030901.

［6］ ?Lu B， Wu B Y， Gu J F， ?et al . Distributed optical fiber hydrophone based on Φ-OTDR and its field test ［J］. Opt Express， 2021， 29： 3417.

［7］ ?Pang F F， He M， Liu H， ?et al . A fading-discrimination method for distributed vibration sensor using coherent detection of phi-OTDR ［J］. IEEE Photonic Tech L， 2016， 23： 28.

［8］ ?Wang X， Lu B， Wang Z Y， ?et al . Interference-fading-free Φ-OTDR based on differential phase shift pulsing technology ［J］. IEEE Photonic Tech L， 2018， 1： 99.

［9］ ?Zhou J， Pan Z Q， Ye Q， ?et al . Characteristics and explanations of interference fading of a φ-OTDR with a multi-frequency source ［J］. J Lightwave Technol， 2013， 31： 2947.

［10］ ?Zhou J， Pan Z Q， Ye Q， ?et al . Phase demodulation technology using a multi-frequency source for discrimination of interference-fading induced false alarms in a φ-OTDR system ［J］. J Lightwave Technol， 2013， 9： 119.

［11］ Chen D， Liu Q， He Z. Phase-detection distributed fiber-optic vibration sensor without fading-noise based on time-gated digital OFDR ［J］. Opt Express， 2017， 25： 8315.

［12］ Hartog A H， Liokumovich L B， Ushakov N A， ?et al . The use of multi frequency acquisition to significantly improve the quality of fiber-optic-distributed vibration sensing ［J］. Geophys Prospect， 2018， 66： 192.

［13］ Zabihi M， Chen Y， Zhou T， ?et al . Continuous fading suppression method for Φ-OTDR systems using optimum tracking over multiple probe frequencies ［J］. J Lightwave Technol， 2019， 37： 3602.

［14］ Zhang J， Wu H， Zheng H， ?et al . 80-km fading free phase-sensitive reflectometry based on multi-carrier NLFM pulse without distributed amplification ［J］. J Lightwave Technol， 2019， 37： 4748.

［15］ Wakisaka Y， Iida D， Oshida H， ?et al . Fading suppression of Φ-OTDR with the new signal processing methodology of complex vectors across time and frequency domains ［J］. J Lightwave Technol， 2021， 39： 4279.

［16］ Zhang Y X， Liu J X， Xiong F， ?et al . A space-division multiplexing method for fading noise suppression in the Φ-OTDR system ［J］. Sensors， 2021， 21： 1694.

［17］ Wu M S， Fan X Y， Liu Q W， ?et al . Highly sensitive quasi-distributed fiber-optic acoustic sensing system by interrogating a weak reflector array ［J］. Opt Lett， 2018， 43： 3594.

［18］ Zhu F， Zhang Y X， Xia L， ?et al . Improved Φ-OTDR sensing system for high-precision dynamic strain measurement based on ultra-weak fiber Bragg grating array ［J］. J Lightwave Technol， 2015， 33： 4775.

［19］ Li H， Sun Q Z， Liu T， ?et al . Ultra-high sensitive quasi-distributed acoustic sensor based on coherent OTDR and cylindrical transducer ［J］. J Lightwave Technol，2020， 38： 929.

［20］ Wu Y， Wang Z N， Xiong J， ?et al . Interference fading elimination with single rectangular pulse in Φ-OTDR ［J］. J Lightwave Technol， 2019， 37： 3381.

［21］ He H J， Yan L S， Qian H， ?et al . Suppression of the Interference fading in phase-sensitive OTDR with phase-shift transform［J］. J Lightwave Technol， ?2021， 39： 295.

［22］ Qian H， Luo B， He H J， ?et al . Fading-free Φ-OTDR evaluation based on the statistical analysis of phase hopping ［J］. Appl Optics， 2022， 61： 6729.

［23］ Han L， Wu J. The study of noise filtering algorithm experiment on spatial domain and frequency domain ofhyperspectral image ［C］∥International Symposium on Spatial Analysis， Spatial-Temporal Data Modeling， and Data Mining. Belingham： SPIE， 2009： 1313.

［24］ Hou J H， Tian J W， Liu J. Wavelet and spatial domain method for image filtering ［J］. Infrared and Laser Engineering， 2006， 35： 122.

［25］ 胡松， 江小煒， 楊光， 等. 移動(dòng)平均濾波在微弱脈沖信號檢測中的應(yīng)用［J］. 計(jì)算機(jī)與數(shù)字工程， 2007， 35： 3.