張必明， 蔣奇云， 莫丹， 肖龍英

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長沙 410083 2 湖南涉外經(jīng)濟(jì)學(xué)院信息科學(xué)與工程學(xué)院, 長沙 410205

?

電磁勘探數(shù)據(jù)粗大誤差處理的一種新方法

張必明1,2， 蔣奇云1*， 莫丹1， 肖龍英1

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長沙 410083 2 湖南涉外經(jīng)濟(jì)學(xué)院信息科學(xué)與工程學(xué)院, 長沙 410205

由于存在各種地電干擾，電磁法勘探采集到的原始電場數(shù)據(jù)中往往包含粗大誤差.電磁法勘探中信號量的測量與傳統(tǒng)的精密測量在誤差來源與特點、測量值分布等方面均存在較大差異.經(jīng)試驗，對電磁勘探采集到的原始電場數(shù)據(jù)采用傳統(tǒng)的萊伊達(dá)、格拉布斯、狄克遜等準(zhǔn)則進(jìn)行粗大誤差的自動判別和剔除，處理效果不好；采用Robust估計和中值濾波方法，也不能達(dá)到滿意的效果；采用手工方式挑揀剔除粗大誤差，處理效率太低，均不能滿足電磁勘探數(shù)據(jù)預(yù)處理的要求.作者提出了一種自適應(yīng)雙向均方差閾值法實現(xiàn)電磁勘探數(shù)據(jù)粗大誤差的自動判別和剔除，此方法對采集到的原始電場數(shù)據(jù)樣本進(jìn)行排序后，采用迭代或遞歸的方式，每次均以中點為界分別計算前后兩部分?jǐn)?shù)據(jù)的均方差，將較大的一個與預(yù)先設(shè)置的均方差閾值進(jìn)行比較，若其大于閾值，則判斷粗大誤差存在于相應(yīng)的一端，進(jìn)而剔除相應(yīng)端端點位置的數(shù)據(jù)點；若前后均方差值都小于閾值或樣本數(shù)量小于3個時算法結(jié)束.此方法具有自適應(yīng)優(yōu)化、閾值參數(shù)化控制、適應(yīng)小樣本數(shù)據(jù)以及計算簡單效率高等特點.大量實驗結(jié)果表明：在選取均方差閾值在30至90范圍內(nèi)時(經(jīng)驗值)，能夠有效地剔除電磁勘探原始電場數(shù)據(jù)中的粗大誤差，保留最可信數(shù)據(jù).目前已在多個實際勘探生產(chǎn)項目中應(yīng)用此方法處理粗大誤差，取得了令人滿意的處理效果.

電磁勘探； 粗大誤差處理； 均方差閾值

1 引言

在電磁法勘探中，接收機接收的電場信號中除了有效信號之外，同時還包含大量干擾信號，有效信號幅值往往在微伏級或毫伏級，而干擾信號特別是工業(yè)用電所造成的干擾信號幅度卻很大，往往在毫伏級甚至達(dá)到伏特級，干擾信號的頻率范圍也非常寬，從直流到幾十赫茲甚至上百兆赫茲范圍內(nèi)都存在(蔣奇云，2010)，因此，采集到的電磁勘探原始電場數(shù)據(jù)中往往混有一些偏離均值較大的粗大誤差數(shù)據(jù).一般采用在每個測點的各個頻率點都進(jìn)行冗余采集，再將剔除粗大誤差后的均值作為測量結(jié)果的方法來提高采集數(shù)據(jù)的可靠性和可信度.與傳統(tǒng)的精密測量(如機械工件尺寸測量等)不同，電磁勘探采集到的原始數(shù)據(jù)值的分布范圍要比精密測量大得多；尤其在采集低頻信號時(0.1 Hz以下)，由于信號周期長，野外勘探采集到的原始數(shù)據(jù)樣本集規(guī)模(即樣本數(shù)量)較??；在某些干擾源較多、干擾頻率范圍較寬以及干擾信號較強的工區(qū)勘探時，數(shù)據(jù)樣本中常常混有較大比例且偏離度較大的粗大誤差數(shù)據(jù)點.

目前，常見的粗大誤差處理方法從原理上大致可分為三類：基于統(tǒng)計學(xué)的方法；基于灰度理論的方法；基于其他理論的方法(如小波理論等).其中以基于統(tǒng)計學(xué)的方法較多，有經(jīng)典的粗大誤差判別準(zhǔn)則(萊伊達(dá)準(zhǔn)則、格拉布斯準(zhǔn)則、狄克遜準(zhǔn)則等)的方法(馬宏和王金波，2009)，也有基于RUBOST估計的穩(wěn)健統(tǒng)計方法(趙黎明，1992)，以及基于其他統(tǒng)計量的判別方法.如劉宗寶等(2012)采用結(jié)合萊伊達(dá)準(zhǔn)則和滑動窗口方法實現(xiàn)粗大誤差剔除和數(shù)據(jù)平滑處理，熊艷艷和吳先球(2010)對比了四種常用粗大誤差判別準(zhǔn)則，耿素軍和余劍(2005)應(yīng)用格拉布斯準(zhǔn)則在智能系統(tǒng)中處理粗大誤差，劉志平(2009)應(yīng)用格拉布斯和狄克遜準(zhǔn)則在檢測儀表中處理粗大誤差，梁暄(2006)和余劍(2003)采用統(tǒng)計估值理論處理粗大誤差等.在小樣本數(shù)據(jù)中，王廣林和李云峰(2005)、高寧等(2012)采用灰度系統(tǒng)理論進(jìn)行粗大誤差的處理.其他理論方面，高用竹和蔡冬曼(2011)、李登輝和李智(2009)采用基于信息熵的方法判別粗大誤差，吳維勇和王英惠(2005)采用小波變換方法去除粗大誤差，劉振宇等(2013)采用基于數(shù)據(jù)拓?fù)渚嚯x的方法進(jìn)行粗大誤差處理.在國外，Carvalho 和Bretas(2013)在電力系統(tǒng)控制中采用最大標(biāo)準(zhǔn)化殘差健壯性分析來處理加權(quán)最小二乘法狀態(tài)估計中的粗大誤差，Bretas等(2011)提出采用一個“創(chuàng)新指數(shù)”指標(biāo)來處理電力系統(tǒng)狀態(tài)估計中的測量粗大誤差，Manescu(2006)在電力系統(tǒng)狀態(tài)估計中采用修改的松弛方法檢測粗大誤差，Sproule等(2006)在澳大利亞陸地重力數(shù)據(jù)庫中采用DEM-9S和經(jīng)局部插值精化處理的布格重力異常進(jìn)行交叉驗證的方法檢測粗大誤差，Nyrnes等(2005)在油井產(chǎn)量可靠性分析中，采用了一種基于統(tǒng)計測試的方法檢測井眼定向測量中的粗大誤差，Badria等(2012)采用統(tǒng)計質(zhì)量控制方法來檢測和處理測量數(shù)據(jù)中的粗大誤差.以上方法中，傳統(tǒng)的基于統(tǒng)計學(xué)的粗大誤差處理方法主要適用于測量樣本量較大，基本呈正態(tài)分布，測量樣本值分布范圍集中，粗大誤差比例小等情況，如傳統(tǒng)的靜態(tài)精密測量，包括機件尺寸等靜態(tài)物理量狀態(tài)測量等.由于電磁勘探數(shù)據(jù)具有數(shù)值分布范圍較大，測量樣本量可能較小(低頻段)，數(shù)據(jù)分布不一定呈正態(tài)分布等特點.經(jīng)試驗，在采用傳統(tǒng)統(tǒng)計方法處理時，數(shù)據(jù)點被剔除的比例過大，處理結(jié)果不理想.

中值濾波作為一種具有穩(wěn)健特性的非線性濾波方法，在信號處理領(lǐng)域一直以來有著廣泛的應(yīng)用，尤其是在數(shù)字圖像處理領(lǐng)域有著非常多的基于中值濾波的方法研究和應(yīng)用(劉麗梅等，2004；索俊祺，2010；王歡和王修暉，2011).在地球物理勘探數(shù)據(jù)預(yù)處理階段，也經(jīng)常采用一維中值濾波進(jìn)行信號去噪(地震法中應(yīng)用非常廣泛).自20世紀(jì)90年代起，李承楚(1986)在垂直地震剖面法(VSP) 的資料處理中, 為了分離上行波和下行波, 采用中值濾波方法，并對其理論進(jìn)行了完整的分析和闡述；劉道安(1989)研究了矢量中值濾波方法，以提高傳統(tǒng)中值濾波的處理效率，并在零偏移距VSP資料處理中進(jìn)行了應(yīng)用并取得良好效果；在地震資料疊前處理剔除野值和提高信噪比中，徐常練(1992)利用了中值濾波的特點，使計算效率大為提高，數(shù)據(jù)試驗表明, 疊前中值濾波在自動剔除野值和不正常道、提高疊前資料的信噪比以及消除面波干擾三個方面都有較好的效果.劉財?shù)?2005)將二維多級中值濾波技術(shù)應(yīng)用于地震數(shù)據(jù)隨機噪聲消除，不僅能有效地壓制隨機噪聲，而且在一定程度上削弱了淺層折射干擾.王典(2006)在其博士學(xué)位論文中詳細(xì)研究了中值濾波在地震勘探數(shù)字信號處理中的應(yīng)用.

上述文獻(xiàn)中的各種方法一般也是針對其自身領(lǐng)域數(shù)據(jù)特點所提出，在一定范圍內(nèi)能夠獲得較好的效果.作者在綜合方法原理和電磁勘探數(shù)據(jù)特點的基礎(chǔ)上，提出一種自適應(yīng)的雙向均方差閾值方法實現(xiàn)電磁勘探數(shù)據(jù)粗大誤差處理，經(jīng)實驗驗證，此方法能取得良好的處理效果.

2 自適應(yīng)雙向均方差閾值法

依據(jù)測量理論，測量結(jié)果中可能存在系統(tǒng)誤差、隨機誤差和粗大誤差(馬宏和王金波，2009).對于電磁勘探來說，數(shù)據(jù)采集設(shè)備在使用前需通過一致性標(biāo)定，因而系統(tǒng)誤差基本可控制在能接受的范圍內(nèi)；勘探生產(chǎn)中存在的測量隨機誤差，一般采用冗余采集并計算均值方式可使其降低到能接受的范圍內(nèi)；而由各種不可預(yù)測的地電強干擾等因素造成的粗大誤差數(shù)據(jù)(簡稱粗差)，則必須采取有效的方法進(jìn)行判別和剔除.

2.1 方法原理

勘探采集到的原始數(shù)據(jù)樣本，在干擾較少的情況下，樣本分布一般接近正態(tài)分布(理論上應(yīng)服從正態(tài)分布)，而在干擾源較多，干擾較強的情況下，分布形態(tài)往往不規(guī)則，不呈近似正態(tài)分布.

在將測點原始數(shù)據(jù)樣本排序后，一般靠近兩端點位置的數(shù)據(jù)點偏離均值較大，而位于中點附近的數(shù)據(jù)點則更接近均值(中位值原理)，當(dāng)數(shù)據(jù)中存在粗大誤差時，一般位于兩端位置(極大或極小).

由于干擾因素不同，從幾何意義上分析排序后樣本從中點分成前后兩段的數(shù)據(jù)分布形態(tài)，一般有三種情況(如圖1)：中心點對稱形態(tài)；右側(cè)陡峭形態(tài)；左側(cè)陡峭形態(tài).從統(tǒng)計指標(biāo)量來看，形態(tài)較平緩則均方差較小，形態(tài)較陡峭則均方差較大；反之，均方差較小則形態(tài)較平緩，均方差較大則形態(tài)較陡峭.因此，若要采用迭代的方式判別并剔除樣本中的粗大誤差，則每次均應(yīng)從形態(tài)較陡峭一側(cè)進(jìn)行判別，若符合剔除條件，則刪除端點處數(shù)據(jù)點.遵循上述原則，一般情況下每次迭代剔除的是整個樣本中偏離度最大的數(shù)據(jù)點.在不斷迭代處理的過程中，中間點的位置也不斷變化，樣本分布逐漸向中心點對稱形態(tài)收斂(自適應(yīng)優(yōu)化).進(jìn)一步地，在樣本數(shù)量較少的情況下，此方法在統(tǒng)計學(xué)原理上仍成立，因此適應(yīng)于小樣本數(shù)據(jù)處理.使用本方法時，應(yīng)根據(jù)可信樣本值分布范圍選取一個合適的均方差作為粗大誤差的判別閾值(對于電磁勘探數(shù)據(jù)，閾值范圍可在1到90之間，一般情況下取經(jīng)驗值30)，閾值越大判別條件越寬松(允許更大的樣本離散度)，反之則越嚴(yán)格.

圖1 數(shù)據(jù)樣本典型分布形態(tài)(a)中心點對稱形態(tài)；(b)右側(cè)陡峭形態(tài)；(c) 左側(cè)陡峭形態(tài).Fig.1 Typical distribution shape of sample data(a) Center point symmetry shape; (b) Right steep shape; (c) Left steep shape.

2.2 算法描述

本文提出的自適應(yīng)雙向均方差閾值法先對原始數(shù)據(jù)樣本排序，然后采用迭代(或遞歸)方式，每次均以中點為界分別計算前后兩段數(shù)據(jù)的均方差，判斷并剔除均方差大于閾值的較大一端的端點數(shù)據(jù)點，在前后兩段的均方差值都小于閾值或樣本數(shù)量小于3時算法結(jié)束.算法描述如下：

(1) 將含N個數(shù)據(jù)點的樣本數(shù)據(jù)由小到大排序.

(2) 計算樣本數(shù)據(jù)前段末元素位置m1和后段首元素位置m2(公式(1)).為使前后兩段具有相關(guān)性，需在中點位置相互重疊(N為奇數(shù)時重疊一個數(shù)據(jù)點，為偶數(shù)時重疊兩個數(shù)據(jù)點)，即m1≥m2.

(3) 采用貝塞爾公式(公式(2))分別計算前段[1,m1]和后段[m2,N]的均方差σfront和σrear.

(4) 比較得到σlarger=Max(σfront,σrear)，將σlarger與均方差閾值σthreshold比較.

(5) 若σlarger≤σthreshold，則算法結(jié)束；否則轉(zhuǎn)到步驟(6).

(6) 若σlarger=σfront則刪除樣本的第一個數(shù)據(jù)點P1，否則刪除最后一個數(shù)據(jù)點PN；樣本數(shù)量N=N-1.

(7)若樣本數(shù)量N<3，則算法結(jié)束；否則轉(zhuǎn)到步驟(2).

中間點位置計算公式如下：

(1)

式中INT()為取整函數(shù)，K為奇偶因子(當(dāng)N為偶數(shù)時K=0，奇數(shù)時K=1).

均方差采用貝塞爾公式計算，公式如下：

(2)

(3)

2.3 算法特點分析

根據(jù)前述分析，本算法具有以下特點.

(1) 自適應(yīng)優(yōu)化

算法在處理過程中每次迭代都只刪除位于端點處的一個數(shù)據(jù)點，刪除后樣本數(shù)量減1，在接下來的一次迭代中，中點位置計算將自動調(diào)整前后兩段分界點(公式(1)).在循環(huán)迭代處理中，算法不斷趨向樣本數(shù)據(jù)中分布最集中的一部分(即最可信的部分)，樣本數(shù)據(jù)的分布形態(tài)也不斷趨向中心點對稱形態(tài)，使得算法具有自適應(yīng)優(yōu)化的特點.

(2) 閾值參數(shù)化控制

算法在處理過程中粗大誤差的判斷依據(jù)是均方差閾值，此閾值是算法從外部接收的控制參數(shù)，用于控制判斷粗大誤差的寬嚴(yán)程度.在實際應(yīng)用時，根據(jù)樣本數(shù)據(jù)值的分布范圍選取寬嚴(yán)條件不同的閾值，以獲得較好的處理效果.

(3) 適用于小樣本數(shù)據(jù)

本算法基于統(tǒng)計學(xué)原理，因此樣本數(shù)量越大，處理結(jié)果越接近真值.同時，根據(jù)均方差計算公式和算法中前后段的分段方式，在僅有3個樣本數(shù)據(jù)點的情況下，仍可分成前后兩段計算均方差并進(jìn)行比較和剔除，最后保留兩個數(shù)據(jù)樣本，理論上適用于小樣本數(shù)據(jù).

(4) 適用于誤差比例較大樣本

根據(jù)算法原理和實驗驗證，本方法在樣本中的誤差比率較大的情況下，仍能有效地判別和剔除粗大誤差數(shù)據(jù)，而相比較的其他方法均已失效或效果較差.

(5) 算法實現(xiàn)簡單，處理效率高

算法中采用貝塞爾公式(公式(2))計算均方差，每次迭代時只計算兩個均方差值，計算代價較小、效率較高；同時，算法的過程較為簡單，容易實現(xiàn)，適用于各種計算能力一般的計算設(shè)備.

2.4 閾值選取方法

均方差作為一個統(tǒng)計指標(biāo)量，其大小表示的是樣本數(shù)據(jù)集的離散程度，值越大表明樣本數(shù)據(jù)越分散，越小則越集中，對于測量數(shù)據(jù)來說，數(shù)據(jù)分布越集中其質(zhì)量相對越好.本文算法中使用的均方差閾值作為一個參數(shù)傳遞給算法，其大小直接影響最終的處理結(jié)果，是一個關(guān)鍵參數(shù)，在使用時須給定一個具體數(shù)值，如何確定這個給定的閾值，目前沒有通用的計算公式，而是根據(jù)應(yīng)用中樣本數(shù)據(jù)的分布范圍由試驗得出.

電磁勘探采集到的有效信號幅值往往在微伏級或毫伏級(蔣奇云，2010)，反映在數(shù)值上為幾十到幾千之間(實際情況一般為10至3000之間)，在此范圍之外的數(shù)據(jù)值一般為干擾信號，實際上，一個勘探測點的某一個頻率點的測量樣本數(shù)據(jù)值分布一般會比上述范圍更小.經(jīng)過反復(fù)試驗，采用按90，75，60，45，30，15，10，5逐步收緊(數(shù)值減小)的方式可選取出適合電磁勘探數(shù)據(jù)預(yù)處理的均方差閾值，大多數(shù)情況下在30至90范圍內(nèi)總可以選取到適合的閾值(經(jīng)驗值).在實際應(yīng)用時，一般選取閾值為30進(jìn)行處理，若處理結(jié)果不理想，導(dǎo)致多數(shù)頻率點的樣本數(shù)據(jù)點剔除比例過大，則可按相反順序逐步放寬閾值(數(shù)值加大)進(jìn)行試驗，反之則可逐步收緊.

2.5 數(shù)據(jù)質(zhì)量評價

一般情況下，電磁勘探采集到的原始數(shù)據(jù)均由儀器設(shè)備采集，從信號感應(yīng)到采樣量化等過程無人工干預(yù).為保證原始數(shù)據(jù)的真實性，在進(jìn)行數(shù)據(jù)預(yù)處理時，只是剔除粗大誤差數(shù)據(jù)點，而不對樣本數(shù)據(jù)值進(jìn)行任何修改.

數(shù)據(jù)質(zhì)量的好壞，可從兩個方面綜合評價.一是按傳統(tǒng)的統(tǒng)計學(xué)方法，根據(jù)各頻點數(shù)據(jù)樣本的相對均方誤差大小評價.相對均方誤差為均方差與算術(shù)均值之比的百分比值(公式(4))，一般而言在10%及以下可以接受，5%及以下為良好.二是在對數(shù)坐標(biāo)系中觀察測點的頻率數(shù)值曲線的形態(tài)，一般頻率曲線形態(tài)變化平緩(即曲線較為平滑)表明數(shù)據(jù)質(zhì)量相對較好，如圖2中0.03125 Hz及以上頻率，反之，若曲線在某些頻點位置有明顯波折、上下尖峰或呈鋸齒狀，則相應(yīng)頻點的數(shù)據(jù)質(zhì)量不良，如圖2中0.03125 Hz以下頻率.

為直觀了解數(shù)據(jù)相對均誤差，本文中大部分圖中在頻率曲線的頻率點位置一般繪有一條垂直的淺灰色線條，此線條代表相應(yīng)頻點樣本數(shù)據(jù)的相對均方誤差，通常稱為誤差棒(見圖2)，誤差棒越長表示相對均方誤差越大，數(shù)據(jù)質(zhì)量相對越差，反之則越好.

相對均方誤差計算公式如下：

(4)

3 實驗結(jié)果與分析

本算法的研究目的，主要為解決電磁勘探生產(chǎn)項目中的數(shù)據(jù)處理的實際問題，限于篇幅，選取兩個項目的兩個測點數(shù)據(jù)進(jìn)行實驗與分析.為簡化文字，在本文的余下部分均以“本文方法”指代本文提出的自適應(yīng)雙向均方差域值法.

3.1 本文方法與經(jīng)典粗大誤差處理方法實例處理對比與分析

以某頁巖氣電磁勘探項目中某個測點數(shù)據(jù)為例進(jìn)行本文方法與經(jīng)典粗大誤差處理方法的處理結(jié)果分析對比研究.

圖2 原始數(shù)據(jù)頻率曲線圖Fig.2 Frequency curve of raw data

圖2為本實例原始數(shù)據(jù)的頻率曲線圖.如圖所示，此測點0.03125 Hz以上頻段頻率曲線較為平滑，且相對均方誤差小(多數(shù)5%以下，部分10%以下)，表明此頻段數(shù)據(jù)質(zhì)量較好；而0.03125 Hz及以下低頻部分曲線呈現(xiàn)明顯尖峰形態(tài)(尖峰出現(xiàn)在0.023438 Hz頻率點處)，且各頻點相對均方誤差較大(見誤差棒)，數(shù)據(jù)質(zhì)量較差.對0.03125 Hz及以下各頻率點(共4個)原始樣本數(shù)據(jù)分布情況進(jìn)行對比分析，如圖3所示，各頻率點樣本數(shù)量范圍為8～24個，分布形態(tài)為非近似正態(tài)分布.

在實驗中，對此測點的所有頻點樣本數(shù)據(jù)分別用本文方法和三種經(jīng)典的基于統(tǒng)計準(zhǔn)則的粗差處理方法，即萊伊達(dá)、格拉布斯以及狄克遜準(zhǔn)則方法，進(jìn)行處理結(jié)果的對比分析.對于三種經(jīng)典準(zhǔn)則方法，萊伊達(dá)準(zhǔn)則是以測量樣本充分大為前提(N?10時適用)，在樣本較小時可靠性不高；當(dāng)樣本較小時格拉布斯準(zhǔn)則可靠性最高(當(dāng)20

采用本文方法(閾值=30)和三種經(jīng)典準(zhǔn)則方法處理結(jié)果對比如圖4所示，四種方法對0.03125 Hz及以上數(shù)據(jù)質(zhì)量較好部分的頻率曲線形態(tài)基本沒有影響，而在0.03125 Hz以下的低頻部分，本文方法明顯消除了頻率曲線的尖峰形態(tài)，使頻率曲線呈平緩變化形態(tài)，同時明顯降低了各頻點的相對均方誤差；萊伊達(dá)和狄克遜準(zhǔn)則方法基本上對尖峰形態(tài)無效，且對相關(guān)頻率點相對均方誤差無改善；從頻率曲線形態(tài)觀察，格拉布斯準(zhǔn)則有效地消除了尖峰，同時所有頻點(全部40個頻率點)的相對均方誤差都變得非常小，看似處理效果很好，但從剔除的樣本點數(shù)量和比例來分析，見表1，格拉布斯準(zhǔn)則處理后每個頻率點均只保留了兩個數(shù)據(jù)樣本(表1 未列出的其余各個頻率點也都僅只保留了兩個數(shù)據(jù)樣本)，其余均認(rèn)為是粗大誤差而進(jìn)行了剔除，而實際情況并非如此，如圖3所示，此例中0.03125 Hz及以下頻率范圍的每個頻率點的有效樣本均不止兩個，因此格拉布斯準(zhǔn)則的條件對電磁勘探數(shù)據(jù)來說過于嚴(yán)格，導(dǎo)致了樣本數(shù)據(jù)的過度剔除，其結(jié)果將影響后續(xù)的勘探數(shù)據(jù)處理和解釋.此例中三種經(jīng)典準(zhǔn)則方法的適用前提條件、判別準(zhǔn)則條件等是其不能適用于電磁勘探數(shù)據(jù)處理的根本原因.

圖片說明(本圖及以下所有數(shù)據(jù)樣本分布圖均參照此說明):1) 各坐標(biāo)圖中左上部均直方圖為數(shù)據(jù)樣本分布直方圖(樣本數(shù)量小于5個時無分布直方圖)；2) 由于各頻率點樣本數(shù)據(jù)數(shù)量和樣本數(shù)據(jù)值范圍可能不同，為了完整顯示樣本數(shù)據(jù)，各頻率樣本坐標(biāo)系中的X和Y坐標(biāo)刻度單位和范圍不一定相同.圖3 0.03125 Hz及以下頻率點原始樣本分布圖(a) 0.03125 Hz；(b) 0.023438 Hz；(c) 0.015625 Hz；(d) 0.0117 Hz.Fig.3 Raw sample data distribution and histogram of frequencies below 0.03125 Hz

圖4 4種方法處理結(jié)果對比(a)本文方法處理結(jié)果(閾值=30)；(b)萊伊達(dá)準(zhǔn)則粗大誤差處理結(jié)果；(c)格拉布斯準(zhǔn)則粗大誤差處理結(jié)果；(d)狄克遜準(zhǔn)則粗大誤差處理結(jié)果.Fig.4 Results comparison of the 4 methods(a) Our method result (threshold=30)；(b) Pauta criterion result； (c) Grubbs criterion result； (d) Dixon criterion result.

進(jìn)一步分析本文方法處理后0.03125 Hz及以下各頻點的樣本數(shù)據(jù)，如圖5所示，處理后原始數(shù)據(jù)樣本中離群較大的粗大誤差數(shù)據(jù)點均已有效剔除，結(jié)果數(shù)據(jù)樣本分布形態(tài)已近似正態(tài)分布.

為了能夠在平滑頻率曲線和降低各頻率點相對均方誤差的情況下盡可能多地保留可信的樣本數(shù)據(jù)，可選取不同寬嚴(yán)程度的均方差閾值進(jìn)行進(jìn)一步處理，并對結(jié)果進(jìn)行比較和分析.下面是選取90、60、30和15四個不同的閾值時本文方法處理結(jié)果的對比，見圖6，隨著閾值的逐步減小，各個頻率點的相對均方誤差相應(yīng)變小(相對均方誤差與均方差直接相關(guān))，頻率曲線的整體平滑形態(tài)變化不大，在此例中選取閾值為30時即可取得滿意的效果，并盡可能多地保留了各頻率點的樣本數(shù)據(jù).

3.2 本文方法與Robust估計方法實例處理對比與分析

以某油氣資源電磁法勘探項目中某個測點數(shù)據(jù)為例，對比原始數(shù)據(jù)、Robust估計結(jié)果以及本文方法的處理結(jié)果并進(jìn)行分析.

圖7a為本測點原始數(shù)據(jù)采用算術(shù)均值計算的頻率曲線，由圖可見，此測點采集到的數(shù)據(jù)質(zhì)量不佳，反映在頻率曲線上，除48 Hz的工頻干擾形成一個尖峰突起外，0.03125 Hz以下頻率由于干擾較大(相對均方誤差很大)、樣本數(shù)量小，頻率曲線形態(tài)也呈現(xiàn)出明顯的尖峰突起.圖7b為原始數(shù)據(jù)采用Robust估計后的曲線形態(tài)，與原始數(shù)據(jù)對比，首先，消除了48 Hz頻率點處由粗大誤差造成的尖峰，其次，0.03125 Hz以下頻率估計后的曲線形態(tài)較原始數(shù)據(jù)也明顯改善，但粗大誤差的影響仍未完全消除，仍有較明顯的小尖峰存在，處理結(jié)果的相對均方誤差也明顯降低.圖7c為本文方法處理結(jié)果數(shù)據(jù)(采用算術(shù)均值計算)的頻率曲線形態(tài)，除了完全消除48 Hz頻率點處的尖峰外，0.03125 Hz以下頻段曲線形態(tài)也較Robust估計結(jié)果好，完全消除了粗大誤差造成的影響，曲線形態(tài)更加平滑，能夠正確地反映實際信號變化趨勢，同時相對均方誤差也明顯降低.表2中列出了部分頻點的處理結(jié)果對比數(shù)據(jù)，并在圖8和圖9中分別顯示了相應(yīng)頻點樣本處理前后的數(shù)據(jù)分布圖.

表1 0.03125 Hz及以下各頻率點樣本數(shù)量及剔除比例對比分析表

注: 除原始樣本列外，每種方法標(biāo)題下面的兩列分別是處理后結(jié)果中的樣本數(shù)據(jù)數(shù)量和判別為粗大誤差并被剔除的樣本比例.

圖5 0.03125 Hz及以下頻率點樣本本文方法處理后數(shù)據(jù)分布圖(a) 0.03125 Hz；(b) 0.023438 Hz；(c) 0.015625 Hz；(d) 0.0117 Hz.Fig.5 Result distribution and histogram comparison of frequencies below 0.03125 Hz handled by our method

圖6 4種不同閾值處理結(jié)果曲線形態(tài)對比圖(a) 閾值=90；(b) 閾值=60；(c) 閾值=30；(d) 閾值=15.Fig.6 Results curve shapes comparison of 4 different threshold values(a) Threshold=90；(b) Threshold=60；(c) Threshold=30；(d) Threshold=15.

頻率值原始樣本大小原始樣本算術(shù)均值Robust估計均值本文方法結(jié)果算術(shù)均值剔除數(shù)剔除比48Hz120個466．2581174．6878176．405316個13．33%0．023438Hz12個982．8742108．949472．16061個8．33%0．015625Hz10個369．795387．902177．74841個10%0．0117Hz10個117．159176．981374．61521個10%

總體而言，在本例中Robust估計和本文方法均較好地改善了原始數(shù)據(jù)中數(shù)據(jù)質(zhì)量較差的大部分頻點，只在0.03125 Hz以下頻率范圍內(nèi)兩種方法的效果有所區(qū)別.根據(jù)圖8、圖9和表2分析，首先，全部頻率點的數(shù)據(jù)樣本數(shù)量分為兩類：以48 Hz為代表的頻率點，其樣本數(shù)量較大，一般都有100個以上；以0.03125 Hz以下頻率為代表的頻率點，樣本數(shù)量較小，不超過20個.其次，觀察處理前后的樣本分布直方圖(見圖8、圖9)，處理后的樣本分布形態(tài)較原始形態(tài)更趨向正態(tài)分布.

根據(jù)表2分析，在樣本數(shù)量較大的頻率區(qū)間，兩種方法處理結(jié)果比較接近，處理效果都比較好.而在0.03125 Hz以下頻率范圍，在樣本數(shù)量較小，粗大誤差數(shù)值較大的情況下本文方法的處理結(jié)果明顯比Robust估計效果更好.

圖10是3.1節(jié)實驗用例采用本文方法和Robust估計結(jié)果的對比圖.由圖可見，0.023438 Hz頻率處的尖峰突起在Robust估計中無法有效消除，據(jù)圖3、圖5和表1分析，0.023438 Hz頻率點的樣本小(8個)，粗大誤差比例高(37.5%)且粗大誤差數(shù)值較大，是Robust估計方法無效的原因.

圖7 本文方法與Robust估計處理結(jié)果對比圖(a) 原始數(shù)據(jù)算術(shù)均值頻率曲線；(b) 原始數(shù)據(jù)Robust估計均值頻率曲線；(c) 本文方法處理結(jié)果算術(shù)均值曲線；(d) 曲線對比(淺灰:原始，深灰:Robust估計，黑:本文方法).Fig.7 Results comparison of our method and Robust Estimation(a) Raw data frequency curve of arithmetic mean；(b) Raw data frequency curve of Robust Estimation mean；(c) Result data frequency curve of arithmetic mean handled by our method；(d) Curves comparison (light gray: raw, dark gray: Robust Estimation, black: our method).

圖8 部分頻點樣本處理前數(shù)據(jù)分布圖(a) 48 Hz；(b) 0.023438 Hz；(c) 0.015625 Hz；(d) 0.0117 Hz.Fig.8 Pictures of sample data distribution and histogram of some frequencies of raw data

圖9 部分頻點樣本處理后數(shù)據(jù)分布圖(a) 48 Hz； (b) 0.023438 Hz； (c) 0.015625 Hz； (d) 0.0117 Hz.Fig.9 Pictures of sample data distribution and histogram comparison of some frequencies handled by our method

圖10 本文方法與Robust估計處理結(jié)果對比(a) 本文方法處理結(jié)果曲線；(b) Robust估計結(jié)果曲線.Fig.10 Results comparison of our method and Robust Estimation(a) Our method result curve；(b) Robust Estimation result curve.

據(jù)實驗結(jié)果和分析，本文方法在數(shù)據(jù)樣本較大，粗大誤差較小的情況下與Robust估計的結(jié)果比較一致，均有較好的處理效果；在樣本數(shù)量較小，粗大誤差比例較高或誤差數(shù)值很大等情況下，Robust估計往往無法消除粗大誤差的影響，而本文方法仍能取得較好的處理效果，較Robust估計有更好的適應(yīng)性.

3.3 本文方法與中值濾波方法實例處理對比與分析

中值濾波作為一種非線性濾波方法，廣泛應(yīng)用于信號處理領(lǐng)域，是一種有效的經(jīng)典濾波抑噪方法.為使本文方法能夠更全面地與各種傳統(tǒng)經(jīng)典方法進(jìn)行對比，本節(jié)采用標(biāo)準(zhǔn)中值濾波方法(SMF，Standard Median Filtering)對3.1和3.2節(jié)的實例數(shù)據(jù)進(jìn)行處理，并將結(jié)果與本文方法進(jìn)行對比分析.

圖11是3.1節(jié)實驗用例采用標(biāo)準(zhǔn)中值濾波方法設(shè)置窗口大小為3、5、7和9的處理結(jié)果.

如圖所示，4種不同大小的濾波窗口不能消除0.03125 Hz以下頻率段曲線的尖峰(位于0.023438 Hz頻率點位置)，對0.03125 Hz以上頻率的正常曲線形態(tài)基本無影響.實際實驗時進(jìn)行多次迭代濾波或是選擇更大的濾波窗口效果也沒有改善(未附圖).進(jìn)一步查看不同濾波窗口大小處理后0.03125 Hz及以下各個頻率點的樣本數(shù)據(jù)分布情況，如圖12(WS表示濾濾窗口大小，即Window Size的縮寫).

圖12與原始數(shù)據(jù)樣本分布圖(圖3)對比.首先，位于樣本中非端點位置出現(xiàn)的異常離群值，SMF方法能夠在不同大小的窗口條件下有效抑制，如圖12中的子圖(e)、(g)、(i)、(k)、(m)、(o)，在窗口大小為5、7、9時，均能有效抑制0.03125 Hz和0.015625 Hz頻率點中的異常離群值.其次，對于0.023438 Hz頻率點，其樣本分布中有多個連續(xù)異常值(3個異常大值)位于樣本末端，此時受到SMF方法自身特性約束(劉道安，1987)，即使改變?yōu)V波窗口大小也無法有效消除異常，如圖12中的子圖(b)、(f)、(j)、(n)；同樣地，對于0.0117 Hz頻率點，樣本中有一個異常值(異常小值)位于樣本末端，在不同大小濾波窗口條件下也無法有效消除，如圖12中的子圖(d)、(h)、(l)、(p).圖12與本文方法處理結(jié)果樣本數(shù)據(jù)分布圖(圖5)對比，本文方法能夠在不修改樣本數(shù)據(jù)值的情況下，有效判別并剔除離群異常值，保留分布最為集中的可信的樣本結(jié)果.

圖13是3.2節(jié)實驗用例采用標(biāo)準(zhǔn)中值濾波方法設(shè)置窗口大小為3、5、7和9的處理結(jié)果.

與圖7a對比.首先，位于48 Hz處的曲線尖峰形態(tài)已明顯抑制(尚未完全消除)，且隨濾波窗口的加大而更趨平滑；其次，位于0.03125 Hz以下頻率段的尖峰形態(tài)也明顯抑制，其中，0.023438 Hz頻率點隨著濾波窗口的變大，其值越來越小，而0.0117 Hz頻率點對窗口大小變化無反應(yīng)，進(jìn)而造成0.023438 Hz頻率點處的曲線形態(tài)隨濾波窗口變大而逐漸向下突出.進(jìn)一步查看不同濾波窗口大小處理后上述頻率點樣本數(shù)據(jù)分布情況，如圖14.

圖14與原始數(shù)據(jù)樣本分布圖(圖8)對比，同圖12類似，位于樣本中非端點位置出現(xiàn)的異常離群值，SMF方法能夠在不同大小的窗口條件下有效抑制，且隨濾波窗口變大，可濾除的脈沖干擾寬度相應(yīng)增大，但對于樣本中有異常值位于樣本兩端時，無法有效消除，如圖14中的子圖(d)、(h)、(l)、(p).圖14與本文方法處理結(jié)果樣本數(shù)據(jù)分布圖(圖9)對比，本文方法能夠在不修改樣本數(shù)據(jù)值的情況下，有效判別并剔除離群異常值，保留分布最為集中的可信的樣本結(jié)果.

圖11 采用4種不同濾波窗口大小結(jié)果對比(a) 濾波窗口大小=3；(b) 濾波窗口大小=5；(c) 濾波窗口大小=7；(d) 濾波窗口大小=9.Fig.11 Results comparison with 4 different filter window size(a) Filter window size =3；(b) Filter window size =5；(c) Filter window size =7；(d) Filter window size =9.

圖12 0.03125 Hz及以下頻率點樣本標(biāo)準(zhǔn)中值濾波處理后數(shù)據(jù)分布圖(a) 0.03125 Hz, WS=3； (b) 0.023438 Hz, WS=3； (c) 0.015625 Hz, WS=3； (d) 0.0117 Hz, WS=3； (e) 0.03125 Hz, WS=5； (f) 0.023438 Hz, WS=5； (g) 0.015625 Hz, WS=5； (h) 0.0117 Hz, WS=5； (i) 0.03125 Hz, WS=7； (j) 0.023438 Hz, WS=7； (k) 0.015625 Hz, WS=7； (l) 0.0117 Hz, WS=7； (m) 0.03125 Hz, WS=9； (n) 0.023438 Hz, WS=9； (o) 0.015625 Hz, WS=9； (p) 0.0117 Hz, WS=9.Fig.12 Result distribution and histogram comparison of frequencies below 0.03125 Hz handled by SMF

圖13 采用4種不同濾波窗口大小結(jié)果對比(a) 濾波窗口大小=3；(b) 濾波窗口大小=5；(c) 濾波窗口大小=7；(d) 濾波窗口大小=9.Fig.13 Results comparison with 4 different filter window size(a) Filter window size =3；(b) Filter window size =5；(c) Filter window size =7；(d) Filter window size =9.

圖14 部分頻點樣本處理后數(shù)據(jù)分布圖(a) 48 Hz, WS=3； (b) 0.023438 Hz, WS=3； (c) 0.015625 Hz, WS=3； (d) 0.0117 Hz, WS=3； (e) 48 Hz, WS=5； (f)0.023438 Hz, WS=5； (g)0.015625 Hz, WS=5； (h)0.0117 Hz, WS=7； (i)48 Hz, WS=7； (j)0.023438 Hz, WS=7； (k)0.015625 Hz, WS=7； (l) 0.0117 Hz, WS=7；(m) 48 Hz, WS=9；(n) 0.023438 Hz, WS=9； (o) 0.015625 Hz, WS=9； (p) 0.0117 Hz, WS=9.Fig.14 Sample data distribution and histogram comparison of some frequencies handled by SMF

圖15 兩個樣例處理結(jié)果頻率曲線對比(黑:本文方法，灰:SMF中值濾波)(a) 樣例1處理結(jié)果對比；(b) 樣例2處理結(jié)果對比.Fig.15 Results comparison of the 2 examples (Black: our method, Gray: SMF method)(a) Comparison of example 1 results; (b) Comparison of example 2 results.

圖15為3.1和3.2節(jié)的SMF濾波處理結(jié)果與本文方法處理結(jié)果的頻率曲線形態(tài)對比.如圖所示，本文方法處理結(jié)果的曲線形態(tài)在此兩例中要比SMF濾波結(jié)果更好.

從方法理論基礎(chǔ)和上述實驗分析處理結(jié)果來看，本文方法與中值濾波方法的主要區(qū)別在于：

(1)方法原理和應(yīng)用領(lǐng)域不同：中值濾波方法是一種非線性低通濾波方法，在處理過程中，除能夠抑制和消除隨機噪聲、脈沖干擾外，對正常信號的值也會產(chǎn)生影響，即除修改異常值外，也會修改一些正常的樣本值，常常用于時間域(一維信號濾波)或空間域(二維空間信號濾波，如數(shù)字圖像處理)連續(xù)采樣信號的濾波處理；本文方法是一種基于統(tǒng)計原理，對多次觀測得到的測量值樣本進(jìn)行粗大誤差的判別與剔除方法，一般用于隨機變量樣本中判別并剔除離群異常值，而對正常樣本點不做任何修改，一般不用于時間域或空間域的連續(xù)采樣信號處理(剔除異常點后會造成連續(xù)采樣信號間斷缺失).

(2)算法局限性：標(biāo)準(zhǔn)中值濾波算法(SMF)自身存在一個局限，即對于首末任意一端出現(xiàn)異常值(一個或多個)時無法對其進(jìn)行有效處理(由算法定義決定)；本文方法不存在此局限性，即方法與樣本異常數(shù)據(jù)點的位置無關(guān).

(3)算法參數(shù)確定：中值濾波窗口大小的選擇要根據(jù)信號中連續(xù)異常值的寬度(即干擾脈沖寬度)確定，一般需要根據(jù)人工觀察數(shù)據(jù)來確定選用的窗口大小以獲得滿意的處理結(jié)果，可濾除的脈沖寬度應(yīng)小于L/2(L為濾波窗口大小)(劉道安，1989)；本文方法由均方差閾值決定對離群值的容忍度，可由試驗了解應(yīng)用領(lǐng)域的測量數(shù)值分布范圍，從而在相同應(yīng)用領(lǐng)域中可確定較為通用和一致的閾值.

(4)迭代處理：中值濾波處理在不改變窗口大小的情況下可多次迭代處理；本文方法在不縮小閾值的情況下迭代操作無效，即一次處理的結(jié)果和N次迭代處理的結(jié)果一樣.

4 結(jié)論

(1)本文以統(tǒng)計學(xué)理論為基礎(chǔ)，提出了一種粗大誤差處理的新方法，其具有自適應(yīng)優(yōu)化、參數(shù)化控制、適應(yīng)小樣本及大誤差比例數(shù)據(jù)以及算法實現(xiàn)簡單和處理效率高等特點.

(2)根據(jù)電磁勘探采集到的電場數(shù)據(jù)的特點，提出均方差閾值選取方法，并由試驗得出：在30至90范圍內(nèi)選取均方差閾值(經(jīng)驗值)，在作者所處理的電磁勘探頻率范圍為0.01 Hz至8192 Hz的電場數(shù)據(jù)時取得較為滿意的效果.

(3)經(jīng)理論分析和實驗驗證，在電磁勘探數(shù)據(jù)的處理應(yīng)用中，在經(jīng)典的粗誤差判別準(zhǔn)則方法，包括萊伊達(dá)、格拉布斯和狄克遜準(zhǔn)則，不適用的情況下，本文方法具有更好的適應(yīng)性，并能取得滿意的處理效果.

(4)經(jīng)實例數(shù)據(jù)對比實驗分析，在數(shù)據(jù)樣本數(shù)量較大(大于100)時，本文方法與Robust估計方法的結(jié)果非常一致，但在樣本數(shù)量較小(小于20)，粗差比例較大或誤差數(shù)值很大的情況下，本文方法的適應(yīng)性要比Robust估計更強一些，能夠取得更好的處理結(jié)果.

(5)經(jīng)實例數(shù)據(jù)對比實驗分析，在應(yīng)用領(lǐng)域的合理性、算法自身局限性、算法參數(shù)通用性等方面，本文方法比中值濾波更顯合理，得到的處理結(jié)果也更令人滿意.

5 建議與展望

均方差是衡量數(shù)據(jù)樣本集分散程度的一個常用統(tǒng)計指標(biāo).實際上，數(shù)據(jù)樣本集的分散程度指標(biāo)還可通過其他方式計算和表達(dá)，如絕對殘差和、線性擬合直線斜率等.本方法選取均方差作為計算和判別指標(biāo)主要是考慮到此指標(biāo)的常用性和通用性，在將來的應(yīng)用中，也可考慮選擇其他的統(tǒng)計指標(biāo)替代均方差，以獲得更好的計算效率或適應(yīng)更多的應(yīng)用場合.

此方法的原理具有一般性，理論上具備在其他領(lǐng)域推廣應(yīng)用的基礎(chǔ).若在其他領(lǐng)域應(yīng)用時，需根據(jù)測量數(shù)據(jù)可信值的分布范圍，經(jīng)試驗后選取與之匹配的閾值，可期待取得良好的處理效果.此方法目前主要應(yīng)用于數(shù)據(jù)采集后的預(yù)處理階段，若要在采集過程中應(yīng)用此方法實時檢測和剔除粗大誤差，也可通過軟件編程方式進(jìn)行實驗和實現(xiàn).此方法處理時首先需要對原始數(shù)據(jù)集進(jìn)行排序，處理完成后結(jié)果數(shù)據(jù)已非原序，若需要在處理完成后恢復(fù)為原序，作者已通過程序數(shù)據(jù)結(jié)構(gòu)設(shè)計的方法予以了解決(另有文章發(fā)表).

致謝 本文在撰寫過程中，得到了導(dǎo)師何繼善院士的指導(dǎo)，何老師以杖朝之年在百忙之中仔細(xì)審閱初稿，并指出了文章中一些不夠嚴(yán)謹(jǐn)?shù)牡胤?，?jīng)作者多次修改最終完成了本文的撰寫工作，在此非常感謝何老師對學(xué)生的指導(dǎo).同時，李帝銓老師在論文的修改過程給予了幫助，在此表示感謝.

Bretas N G, Bretas A S, Piereti S A. 2011. Innovation concept for measurement gross error detection and identification in power system state estimation.IETGeneration,Transmission&Distribution, 5(6): 603-608.

Carvalho B, Bretas N. 2013. Analysis of the largest normalized residual test robustness for measurements gross errors processing in the WLS state estimator.Systemics,CyberneticsandInformatics, 11(7): 1-6.

Elgazooli B A G, Ibrahim A M. 2012. Multiple gross errors detection in surveying measurements using statistical quality control.JournalofScienceandTechnology, 13: 36-47.

Gao N, Cui X M, Gao C Y. 2012. Grey envelope detection method for multidimensional gross error of deformation data.JournalofGeodesyandGeodynamics(in Chinese), 32(6): 99-102.

Gao Y Z, Cai D M. 2011. Non-statistics judging method for gross error in measurement.Science&TechnologyInformation(in Chinese), (30): 59, doi: 10.3969/j.issn.1672-3791.2011.30.047.

Geng S J, Yu J. 2005. Gross error processing in intelligent measuring system.JournalofElectrical&ElectronicEducation(in Chinese), 27(3): 37-39, doi: 10.3969/j.issn.1008-0686.2005.03.010. Jiang Q Y. 2010. Study on the key technology of wide field electromagnetic sounding instrument [Ph. D. thesis] (in Chinese). Changsha: Central South University. Li C C. 1986. The basic theory of median filtering.OilGeophysicalProspecting(in Chinese), 21(4): 372-379.

Li D H, Li Z. 2009. Processing of gross error in small samples based on measurement information theory.MechanicalEngineering&Automation(in Chinese), (6): 115-117. Liang X. 2006. A processing method of gross error in intelligent electronic measuring system.JournalofShanxiNormalUniversity(NaturalScienceEdition) (in Chinese), 20(4): 46-49, doi: 10.3969/j.issn.1009-4490.2006.04.012.

Liu C, Wang D, Liu Y, et al. 2005. Preliminary application study on of using 2D multi-level median filtering technology on eliminating random noise.OilGeophysicalProspecting(in Chinese), 40(2): 163-167.

Liu D A. 1989. Vector median filtering and its application.OilGeophysicalProspecting(in Chinese), 24(4): 460-468.

Liu L M, Sun Y R, Li L. 2004. Research of median filter technology.JournalofYunnanNormalUniversity(in Chinese), 24(1): 23-27. Liu Z B, Gao S Q, Du J Q. 2012. An algorithm for measurement data eliminating gross error and smoothing (in Chinese).∥ Proceedings of the 31st Chinese Control Conference. 7582-7587. Liu Z P. 2009. Treatment analysis of rough error about measurement instrument.ElectronicMeasurementTechnology(in Chinese), 32(11): 55-57, doi: 10.3969/j.issn.1002-7300.2009.11.016.

Liu Z Y, Zhou S H, Tan J R, et al. 2013. Multiple parameters correlation analysis and prediction method of product performance based on multi-criteria modification.JournalofMechanicalEngineering(in Chinese), 49(15): 105-114.

Ma H, Wang J B. 2009. Instrument Precision Theory (in Chinese). Beijing: Beihang University Press, 80-87.

Manescu L G. 2006. Detection of gross error by modified relaxation state estimation in power system.AnnalsoftheUniversityofCraiova,ElectricalEngineeringSeries, 30: 275-281.

Nyrnes E, Torkildsen T, Nahavandchi H. 2005. Detection of gross errors in wellbore directional surveying for petroleum production with emphasis on reliability analyses.KartOgPlan, 65(2): 1-19. Sproule D M, Featherstone W E, Kirby J F. 2006. Localised gross-error detection in the Australian land gravity database.ExplorationGeophysics, 37(2): 175-179.

Suo J Q. 2010. A new improved filtering algorithm based on median filter algorithm [Master′s thesis] (in Chinese). Beijing: Beijing University of Posts and Telecommunications.

Wang D. 2006. Several new digital techniques and its application in seismic prospecting [Ph. D. thesis] (in Chinese). Changchun: Jilin University.

Wang G L, Li Y F. 2005. An effective method for distinguishing the gross error.JournalofComputationandMeasurementTechnology(in Chinese), 25(6): 15-17, doi: 10.3969/j.issn.1674-5795.2005.06.005. Wang H, Wang X H. 2011. An improved adaptive median filter based on noise detection.JournalofChinaUniversityofMetrology(in Chinese), 22(4): 382-385.

Wu W Y, Wang Y H. 2005. Method for eliminating the measured bulky errors based on wavelet transformation.JournalofMachineDesign(in Chinese), 22(6): 61-62, doi: 10.3969/j.issn.1001-2354.2005.06.023.

Xiong Y Y, Wu X Q. 2010. The generalizing application of four judging criterions for gross errors.PhysicalExperimentofCollege(in Chinese), 23(1): 66-68, doi: 10.3969/j.issn.1007-2934.2010.01.022.

Xu C L. 1992. Prestack mid-value filtering.OilGeophysicalProspecting(in Chinese), 27(3): 387-396.

Yu J. 2003. Gross error processing technology in intelligence measuring system with high precision.JournalofTestandMeasurementTechnology(in Chinese), 17(3): 258-261, doi: 10.3969/j.issn.1671-7449.2003.03.017.

Zhao L M. 1992. Robust statistics method applied in sample mixed with outliers.EnvironmentSurveillanceManagementandTechnology(in Chinese), 4(4): 58-59.

附中文參考文獻(xiàn)

高寧, 崔希民, 高彩云. 2012. 變形數(shù)據(jù)多維粗差的灰包絡(luò)線探測方法. 大地測量與地球動力學(xué), 32(6): 99-102.

高用竹, 蔡冬曼. 2011. 測量中粗大誤差的非統(tǒng)計判定方法. 科技資訊, (30): 59, doi: 10.3969/j.issn.1672-3791.2011.30.047. 蔣奇云. 2010. 廣域電磁測深儀關(guān)鍵技術(shù)研究[博士論文]. 長沙: 中南大學(xué).

李承楚. 1986. 關(guān)于中值濾波的理論基礎(chǔ). 石油地球物理勘探, 21(4): 372-379.

李登輝, 李智. 2009. 基于測量信息論的小樣本粗大誤差處理研究. 機械工程與自動化, (6): 115-117.

梁暄. 2006. 一種智能電子測量系統(tǒng)中粗大誤差的處理方法. 山西師范大學(xué)學(xué)報(自然科學(xué)版), 20(4): 46-49, doi: 10.3969/j.issn.1009-4490.2006.04.012.

劉財, 王典, 劉洋等. 2005. 二維多級中值濾波技術(shù)在隨機噪聲消除中的應(yīng)用初探. 石油地球物理勘探, 40(2): 163-167.

劉道安. 1989. 矢量中值濾波及其應(yīng)用. 石油地球物理勘探, 24(4): 460-468.

劉麗梅, 孫玉榮, 李莉. 2004. 中值濾波技術(shù)發(fā)展研究. 云南師范大學(xué)學(xué)報, 24(1): 23-27.

劉宗寶, 高世橋, 杜井慶. 2012. 測量數(shù)據(jù)剔除粗大誤差與平滑處理的一種算法. ∥ 第三十一屆中國控制會議論文集. 7582-7587.

劉志平. 2009. 檢測儀表中粗大誤差的剔除分析. 電子測量技術(shù), 32(11): 55-57, doi: 10.3969/j.issn.1002-7300.2009.11.016.

劉振宇, 周思杭, 譚建榮等. 2013. 基于多準(zhǔn)則修正的產(chǎn)品性能多參數(shù)關(guān)聯(lián)分析與預(yù)測方法. 機械工程學(xué)報, 49(15): 105-114.

馬宏, 王金波. 2009. 儀器精度理論. 北京: 北京航空航天大學(xué)出版社, 80-87. 索俊祺. 2010. 一種新的基于中值濾波的優(yōu)化濾波算法[碩士論文]. 北京: 北京郵電大學(xué). 王典. 2006. 地震勘探幾種數(shù)字新技術(shù)及其應(yīng)用[博士論文]. 長春: 吉林大學(xué).

王廣林, 李云峰. 2005. 一種有效的粗大誤差判別方法. 計測技術(shù), 25(6): 15-17, doi: 10.3969/j.issn.1674-5795.2005.06.005.

王歡, 王修暉. 2011. 基于閾值判斷的自適應(yīng)中值濾波算法. 中國計量學(xué)院學(xué)報, 22(4): 382-385.

吳維勇, 王英惠. 2005. 基于小波變換的粗大誤差去除算法. 機械 設(shè)計, 22(6): 61-62, doi: 10.3969/j.issn.1001-2354.2005.06.023. 熊艷艷, 吳先球. 2010. 粗大誤差四種判別準(zhǔn)則的比較和應(yīng)用. 大學(xué)物理實驗, 23(1): 66-68, doi: 10.3969/j.issn.1007-2934.2010.01.022.耿素軍, 余劍. 2005. 智能測量系統(tǒng)中粗大誤差的處理. 電氣電子教學(xué)學(xué)報, 27(3): 37-39, doi: 10.3969/j.issn.1008-0686.2005.03.010. 徐常練. 1992. 疊前中值濾波. 石油地球物理勘探, 27(3): 387-396. 余劍. 2003. 高精度智能測量系統(tǒng)中粗大誤差的處理技術(shù). 測試技術(shù)學(xué)報, 17(3): 258-261, doi: 10.3969/j.issn.1671-7449.2003.03.017. 趙黎明. 1992. 穩(wěn)健統(tǒng)計法在含離群值樣本中的應(yīng)用. 環(huán)境監(jiān)測管理與技術(shù), 4(4): 58-59.

(本文編輯 何燕)

A novel method for handling gross errors in electromagnetic prospecting data

ZHANG Bi-Ming1,2, JIANG Qi-Yun1*, MO Dan1, XIAO Long-Ying1

1CollegeofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China2CollegeofInformationScienceandEngineering,HunanInternationalEconomicsUniversity,Changsha410205,China

Because of various interference factors, there are more or less gross errors mixed in raw electric field data from electromagnetic (EM) prospecting. The sources, features of such errors, as well as distribution and range of measurements of EM differ much from traditional precise surveys. We have tried several traditional methods in experiments such as Pauta, Grubbs and Dixon criteria to detect and eliminate gross errors in EM prospecting data, but the result is not satisfactory. Even using the Robust Estimation and Median Filtering methods can not obtain satisfactory results either. Meanwhile, the efficiency is very poor by manual operation for this work. So, we need a suitable method/algorithm to solve this problem.We proposed a new adaptive bidirectional MSD (Mean Square Deviation) threshold method to detect and eliminate gross errors in raw data of EM prospecting. This method sorts raw sample data set at first, and then adopts an iterative or recursive way to separately calculate two MSD values repeatedly for the front and rear part of the sample data set which is separated at the array midpoint. Afterwards, it picks the larger one of the two calculated MSDs to compare with the MSD threshold value. If the picked one is larger than the threshold, then we delete the very one data point at endpoint in corresponding part, otherwise the algorithm finished. In addition, if the size of data sample is less than 3 after deletion, the algorithm finished as well. To test and verify this method, we used lots of real raw data of prospecting sites from two industrial projects; one from a shale gas prospecting project, the other from an oil-gas prospecting project. We try various values as the MSD threshold to obtain satisfactory handled results. Consequently the value between 30 and 90 shows pretty good. Based on experimental data, we compare our method with some traditional methods such as Pauta, Grubbs and Dixon criteria, and Robust Estimation and Median Filtering methods.In the first experiment case, bad data quality causes the tail part (below 0.03125 Hz) of full 40 frequencies curve to exhibit a sharp peak shape. After handling gross errors by our new method and Pauta, Grubbs and Dixon criteria, the results show that the new method smoothes the sharp peak shape segment well and keeps the rest of the curve almost unchanged, while the Pauta and Dixon criterion methods have no effect to the sharp peak shape segment, actually producing no change for the whole curve. While the Grubbs criterion method smoothes the sharp peak shape of the curve, but also deletes the most of samples (include most of credible samples) while only 2 samples left for every frequency point after processing. So this result is not good for subsequent procedures such as further data processing and interpretation. To go deep into details, compared the raw sample data set of some frequencies below 0.03125 Hz with handled results of our method, we find that the apparent gross error sample points are eliminated correctly and the most credible samples are preserved. In the same case, we tried 4 typical MSD threshold values as 90, 60, 30 and 15, respectively. The results show that the smaller threshold value cannot make the result obviously better below 30, but more suspected sample points are removed. It means the value 30 is pretty good. In the second experiment case, we compare our method with the classical Robust Estimation method. The result shows that the two methods are quite consistent above 0.03125 Hz, both smooth the sharp peak at 48 Hz and keep the rest curve shape almost no change, but apparent difference appears at the tail part below 0.03125 Hz, in which our method is better than the Robust Estimation. Furthermore, if the sample size is small (less than 20) and the proportion of gross error samples is not small or the value of gross error samples is very big, especially those frequencies below 0.03125 Hz, the Robust Estimation method cannot obtain satisfactory results. At the same time, both the Robust Estimation and our method work pretty well at 48 Hz. In the third experiment case, we adopt the two same experiment data samples used at the first and the second cases and use the SMF (Standard Median Filtering) method to handle these data samples and compared the results with our method. The results show that in first example the SMF has no effect on the peak shape at 0.023438 Hz. The reason is that there are three big outliers just at the position of rear endpoint of the sample and SMF can not handle outliers at this point (either head or tail) as an inherent limitation. The second example is better than the first, but not better than our method. Finally, we analyze differences between Median Filtering and our method at the end of section in this case. Besides these three experiment cases, lots of experiments illustrate that our method can eliminate reliably gross errors mixed in EM data, and keep the most credible data as much as possible. By choosing an appropriate MSD threshold value between 30 and 90, we can achieve satisfactory results.This work suggested a novel method based on the statistical theory. It has some good features such as adaptive optimization, threshold parameterization adaptive for small sample data size and large proportion of gross error sample, being easy to program implementation and low computation cost, high efficiency. The MSD value between 30 and 90 is determined by mass of practice data handling, which be proved appropriate for handling EM prospecting data. Theoretical analysis and experiments show this method has better adaptability and handling results than the classical methods as Pauta, Grubbs and Dixon criteria. Furthermore, in conditions of small sample data size (less than 20), large proportion of gross error and huge values of gross errors, our method is better in method adaptability and can obtain more satisfactory results than Robust Estimation. In other conditions, the handling results of our method and Robust Estimation are quite consistent. Based on experiments and theoretical analysis, in the suitability of the application field, inherent limitation of algorithm and the consistency of parameters of algorithm, our method is also superior to the Median Filtering method.

Electromagnetic prospecting；Gross error handling；MSD threshold

10.6038/cjg20150623.

國家自然科學(xué)基金重大科研儀器設(shè)備研制專項(41227803)資助.

張必明，男，1977年生，湖南懷化人，中南大學(xué)在讀博士，講師，地球探測與信息技術(shù)專業(yè)，研究方向：電磁法勘探數(shù)據(jù)處理、數(shù)據(jù)可視化及并行計算. E-mail:zbm319@163.com

*通訊作者 蔣奇云，男，1974年生，講師，地球探測與信息技術(shù)專業(yè)，研究方向：地球物理儀器與信號處理.E-mail:jqy7412@163.com

10.6038/cjg20150623

P631

2014-02-27，2015-04-22收修定稿

張必明， 蔣奇云， 莫丹等. 2015. 電磁勘探數(shù)據(jù)粗大誤差處理的一種新方法.地球物理學(xué)報，58(6):2087-2102,

Zhang B M, Jiang Q Y, Mo D, et al. 2015. A novel method for handling gross errors in electromagnetic prospecting data.ChineseJ.Geophys. (in Chinese),58(6):2087-2102,doi:10.6038/cjg20150623.