當(dāng)量和爆高(埋深)反演的聲震分析方法

張亮永,張德志,肖衛(wèi)國,梁旭斌,王同東,郭權(quán)勢(shì),李 翱

(強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)全國重點(diǎn)實(shí)驗(yàn)室;西北核技術(shù)研究所：西安710024)

在國防和民用領(lǐng)域中,近地面爆源參數(shù)反演是非常重要的研究內(nèi)容,快速預(yù)測(cè)爆炸當(dāng)量和爆高(埋深)(height of burnst, HOB)等源參數(shù),對(duì)爆炸事故和恐怖襲擊預(yù)警[1-4]、地震和火山噴發(fā)事件預(yù)報(bào)[5-7]、核查監(jiān)測(cè)[8-10]及武器性能評(píng)估[11-12]等方面具有重要意義。采用聲學(xué)數(shù)據(jù)(近區(qū)為沖擊波)[13-21]或地震波數(shù)據(jù)(近區(qū)為地運(yùn)動(dòng))[2-3, 22-24]預(yù)測(cè)爆源參數(shù)是常見的兩種近地面爆源參數(shù)反演方法,但采用這兩種方法反演的爆炸當(dāng)量和爆高(埋深)分布在特定范圍內(nèi)存在多個(gè)極值點(diǎn)(稱為折中關(guān)系[22-27]),給源參數(shù)估計(jì)造成了極大的不確定性。為提高源參數(shù)預(yù)測(cè)的精度,研究人員提出了聲震分析方法[22-27],該方法基于數(shù)據(jù)融合算法綜合地震波(地運(yùn)動(dòng))數(shù)據(jù)和聲學(xué)(空氣超壓)數(shù)據(jù)反演爆源參數(shù),通過同時(shí)利用聲震波動(dòng)信息對(duì)源參數(shù)進(jìn)行多重約束,減小折中關(guān)系的影響。聲震分析方法分為聲學(xué)模型、地震波模型和數(shù)據(jù)融合算法3個(gè)部分,本文將對(duì)這3部分目前的進(jìn)展展開論述。

1 聲學(xué)模型

由于近地面爆炸能量會(huì)耦合到空氣中產(chǎn)生明顯的聲擾動(dòng),一般通過遠(yuǎn)場(chǎng)聲學(xué)數(shù)據(jù),利用超壓峰值、正向聲沖量及正向脈寬等波形特征量來預(yù)測(cè)地面爆炸當(dāng)量[22-27],研究人員相繼提出適合快速計(jì)算,以超壓、聲沖量和持續(xù)時(shí)間為特征量的半經(jīng)驗(yàn)聲學(xué)模型及描述波形細(xì)節(jié)的全波形聲學(xué)反演方法。半經(jīng)驗(yàn)聲學(xué)模型包括ANSI模型[28, 29]、KG85標(biāo)準(zhǔn)模型[30]、BOOM模型[15, 31]和IPM聲學(xué)模型[19]等自由場(chǎng)聲學(xué)模型及地面爆炸聲學(xué)模型[2]和近地面爆炸聲學(xué)模型[22-26]。半經(jīng)驗(yàn)自由場(chǎng)聲學(xué)模型未考慮爆炸能量耦合到地下的影響,對(duì)地面爆炸事件通常假設(shè)兩倍等效爆炸當(dāng)量進(jìn)行爆源參數(shù)估計(jì)[20, 28-30, 32];而對(duì)于近地面爆炸事件,Zhang等[33-34]通過考慮地面反射、耦合作用及大氣環(huán)境因素,建立了基于自由場(chǎng)聲學(xué)模型的近地面聲學(xué)模型(簡稱為自由場(chǎng)近地面聲學(xué)模型),進(jìn)行了爆炸當(dāng)量和爆高估計(jì)。表1為不同的自由場(chǎng)聲學(xué)模型。地面爆炸聲學(xué)模型[2]是基于地面爆炸聲學(xué)數(shù)據(jù)建立的模型,考慮了地面爆炸能量耦合情況,用于地面爆炸事件的當(dāng)量估計(jì),但不適用于偏離地面位置的爆炸事件當(dāng)量估計(jì)。而近地面爆炸聲學(xué)模型考慮了近地面聲震能量耦合特性,可用于近地面爆炸事件的當(dāng)量和爆高(埋深)估計(jì),但不適用比爆高為-2.4 ～-0.6 m·kg-1/3的淺埋情況。上述兩種模型需根據(jù)不同場(chǎng)地介質(zhì)提前確定多個(gè)待定系數(shù),適用的場(chǎng)地介質(zhì)有限。Zhang等[25]在近地面爆炸聲來源機(jī)制[35-36]的基礎(chǔ)上,通過考慮氣體噴出地表和地沖擊耦合影響對(duì)近地面聲學(xué)模型進(jìn)行改進(jìn),拓寬了比例爆高適用范圍,同時(shí)確定了花崗巖夾雜變質(zhì)巖場(chǎng)地的模型待定系數(shù)。圖1為近地面爆炸聲來源機(jī)制[35-36]和改進(jìn)前后近地面聲學(xué)模型[25]。

(a)Mechanism of near-surface explosion producing sound [35-36]

(b)Unimproved and improved near-surface acoustic models[25]

表1 自由場(chǎng)聲學(xué)模型

全波形聲學(xué)反演方法包括超壓波形模板匹配法[16]、經(jīng)驗(yàn)聲源模型法[20]和折合聲沖量法[14, 17, 24]等。超壓波形模板匹配法[16]以不同比距離下隨比時(shí)間變化的大量壓力波形數(shù)據(jù)為模板,將未知當(dāng)量的全壓力波形和模板壓力波形進(jìn)行匹配,對(duì)應(yīng)匹配偏差最小的當(dāng)量為估計(jì)當(dāng)量,但該方法需大量不同當(dāng)量不同比距離的壓力波形組成的樣本庫,但目前樣本數(shù)據(jù)庫還不夠充分;經(jīng)驗(yàn)聲源模型法[20]以G17HE為標(biāo)準(zhǔn)源模型,通過大量實(shí)測(cè)視壓力數(shù)據(jù)并假設(shè)標(biāo)準(zhǔn)高斯分布,獲取源模型的峰值壓力和持續(xù)時(shí)間平均值和方差(偏差分布),由此計(jì)算爆炸當(dāng)量,但該模型聲測(cè)點(diǎn)數(shù)據(jù)需選取20～100 m·kg-1/3的比距離范圍,在未知當(dāng)量時(shí)不好提前預(yù)設(shè)測(cè)點(diǎn)距離。圖2為超壓波形樣本庫[16]以及G17HE源模型和實(shí)測(cè)波形[20]。Kim等[14, 17]基于線性聲學(xué)理論反演爆源等效源時(shí)間函數(shù),以KG85模型折合聲沖量作為參考函數(shù),提出一種地面爆炸當(dāng)量反演思路(本文稱為折合聲沖量法),并將其應(yīng)用于剛性地面爆炸當(dāng)量反演,具有較高的當(dāng)量預(yù)測(cè)精度,但該方法只適用剛性地面爆炸當(dāng)量估計(jì),沒有考慮波傳播過程中地表結(jié)構(gòu)的影響。圖3為等效源時(shí)間函數(shù)及預(yù)測(cè)和實(shí)測(cè)波形。張亮永等[37]對(duì)該方法進(jìn)行了發(fā)展,通過考慮波傳播過程中地表結(jié)構(gòu)的吸收衰減影響及近地面爆炸聲震能量分配過程,建立了多孔彈性場(chǎng)地近地面折合聲沖量模型,解決了多孔彈性場(chǎng)地的當(dāng)量估計(jì)問題,給出當(dāng)量相對(duì)偏差為7%,爆高相對(duì)偏差不超過19%。圖4為多孔彈性場(chǎng)地聲傳播示意圖和爆源參數(shù)反演結(jié)果。

(a)Pressure wave template[16] producing sound [35-36]

(b)Source model of G17HE with measured waveforms [20] acoustic models[25]

圖3 等效源時(shí)間函數(shù)(頂端黑色曲線)及預(yù)測(cè)(紅色曲線)和實(shí)測(cè)(黑色曲線)波形[14]

(a)Sound propagation diagram in three-layer medium model

(b)Yield and height of burst predicted by seismoacoustic analysis

2 地震波模型

近地面爆炸耦合的地震波(近距離稱為地沖擊或地運(yùn)動(dòng))來源于爆炸直接耦合方式和空氣超壓(聲波)耦合方式,對(duì)近距離前者稱為直接地沖擊,后者稱為感生地沖擊[12, 38-39],而對(duì)遠(yuǎn)距離,Albert等[40]將前者耦合的地震波稱為前驅(qū)地震波,后者耦合的地震波稱為聲耦合地震波。圖5為近地面爆炸耦合的地震波來源方式示意圖。

圖5 近地面爆炸耦合的地震波來源方式示意圖[40]

目前,近地面地震波模型以前驅(qū)地震波模型為主,隨著DM/DB[2],HRI-III[16, 23, 41],SAY[22, 42],FSE和SPE[24]等系列化爆試驗(yàn)的開展,相繼建立了沖積土、沉積巖和花崗巖等場(chǎng)地介質(zhì)的近地面前驅(qū)地震波爆炸當(dāng)量估計(jì)模型[23-24, 26, 42]。其中,Templeton 等[42]基于HRI-III 和SAY 化爆試驗(yàn)地震波數(shù)據(jù)對(duì)地震波模型進(jìn)行了深入討論,對(duì)比和分析了沖積土和沉積巖兩種地質(zhì)成分下的近地面爆炸地震波模型,建立了地震波數(shù)據(jù)分析和模型建立的基本方法,對(duì)近地面化爆地震波建模具有指導(dǎo)意義。圖6為不同場(chǎng)地介質(zhì)的近地面地震波模型和實(shí)測(cè)數(shù)據(jù)。但上述化爆試驗(yàn)場(chǎng)地介質(zhì)類型主要是沖積土、沉積巖和花崗巖等,受試驗(yàn)地質(zhì)成分的限制,地震波模型局限在以上幾種地質(zhì)成分,數(shù)據(jù)還不夠充分,模型還不夠完善[22, 24, 42-43],且未考慮爆源附近地質(zhì)和其余位置地質(zhì)不一樣的情況。

(a)Alluvium

(b)Sedimentary rock

Zhang等[25, 34, 44]基于試驗(yàn)數(shù)據(jù)探討了花崗巖夾雜變質(zhì)巖場(chǎng)地和土石混合場(chǎng)地的近地面地震波實(shí)測(cè)數(shù)據(jù)和模型,分別如圖7和圖8所示。

(a)Seismic data

(b)Seismic model

(a)Model one

(b)Model two

(c)Model three

(d)Model four

(e)Model five

為解決爆源附近地質(zhì)和其余位置地質(zhì)不一樣的混合場(chǎng)地地震波模型建模問題,給出了模型的待定系數(shù),并通過引入放大系數(shù)。圖9為放大系數(shù)K不同取值對(duì)反演結(jié)果的影響。圖9中,綠色、紅色虛線分別為聲耦合地震波模型和地震波位移首峰值模型的折中曲線(最小偏差對(duì)應(yīng)的當(dāng)量和埋深),白色五角星為交叉和最小偏差位置。表2為土石混合場(chǎng)地不同模型的定義。

(a)K=2

(b)K=4

(c)K=6

(d)K=10

表2 土石混合場(chǎng)地模型的定義[44]

上述近地面前驅(qū)地震波模型,基于初至波首峰值建立當(dāng)量和爆高(埋深)關(guān)系,適用中近區(qū)爆炸當(dāng)量估計(jì)。對(duì)于遠(yuǎn)區(qū)情況,Pasyanos等[10, 24]基于近地面爆炸地震波能量耦合關(guān)系,采用地震波尾波包絡(luò)法,建立近地面遠(yuǎn)區(qū)地震波模型,進(jìn)一步拓展了近地面地震波模型的適用范圍。圖10為尾波包絡(luò)法示意圖[46]和典型反演結(jié)果[24]。

(a)Schematic diagram of coda-wave method

(b)Typical inversion results

聲耦合地震波模型可解決單一地震計(jì)的近地面爆源當(dāng)量估計(jì)問題,為爆源參數(shù)反演提供除前驅(qū)地震波外的新途徑,但相關(guān)研究較少。Zhang等[46-47]基于聲耦合系數(shù)建立聲耦合地震幅值和聲學(xué)幅值之間的關(guān)系,確定沙土場(chǎng)地的聲震耦合系數(shù),并聯(lián)合近地面聲學(xué)模型建立聲耦合地震波模型。沙場(chǎng)地聲耦合地震波典型波形及耦合系數(shù)衰減特性分別如圖11和圖12所示[46]。表3為聲耦合地震波特征值的定義[46]。

(a)East

(b)North

(c)Verticle

(a)Z1

(b)ENZ1 1&2

(c)ENZ1 0&1

(d)Z2

(e)Zhalf

(f)ENZ2 1&2

(g)ENZhalf1 1&2

(h)ENZ2 0&1

(i)ENZhalf1 0&1

表3 聲耦合地震波特征值的定義[46]

3 數(shù)據(jù)融合算法

近地面爆炸能量會(huì)耦合到大氣和地介質(zhì)中,產(chǎn)生大氣擾動(dòng)(遠(yuǎn)場(chǎng)為聲波)和地運(yùn)動(dòng)(遠(yuǎn)場(chǎng)為地震波)。圖13為近地面爆炸產(chǎn)生的典型聲震信號(hào)[23]。

(a)Seismo-acoustic energy distribution

(b)Seismic data of HRI and HRII

(c)Acoustic data of HRI and HRII

采用聲學(xué)數(shù)據(jù)或地震波數(shù)據(jù)對(duì)源參數(shù)進(jìn)行反演,當(dāng)量和爆高(埋深)之間通常會(huì)存在明顯的折中關(guān)系使求解不適定。而聲震分析方法通過融合地震波(地運(yùn)動(dòng))數(shù)據(jù)和聲學(xué)(空氣超壓)數(shù)據(jù)使得反演問題變得適定,從而得到當(dāng)量和爆高(埋深)估計(jì)值。圖14為典型試驗(yàn) (FSE-4) 的折中關(guān)系和聲震分析結(jié)果[24]。圖14中,五角星為真實(shí)值,空心三角形為地表處的最優(yōu)點(diǎn),實(shí)心三角形為整個(gè)區(qū)域的最優(yōu)點(diǎn)。

(a)Results by seismic data

(b)Results by acoustic data

(c)Results by seismic and acoustic data

聲震融合算法是聲震分析的關(guān)鍵組成,隨著聲震分析方法的不斷改進(jìn),聲震融合算法從相對(duì)偏差法[22]、對(duì)數(shù)偏差法[23]及折中曲線交叉法[23]到似然函數(shù)法[24]和貝葉斯MCMC方法[25-26, 48]再到機(jī)器學(xué)習(xí)方法[27],得到了不斷發(fā)展。表4為數(shù)據(jù)融合算法特點(diǎn)和計(jì)算公式。

表4 數(shù)據(jù)融合算法特點(diǎn)和計(jì)算公式

圖15為不同數(shù)據(jù)融合方法的典型反演結(jié)果。表4列出的前4種方法沒有考慮模型參數(shù)的先驗(yàn)信息,源參數(shù)缺少進(jìn)一步約束,且計(jì)算采用格點(diǎn)搜索法,效率較低。貝葉斯馬爾可夫鏈蒙特卡羅 (Markov chain Monte Carlo,MCMC)方法[25-26, 48]同時(shí)考慮了偏差信息和先驗(yàn)信息,采用統(tǒng)計(jì)的方法對(duì)爆源參數(shù)進(jìn)行估計(jì)可給出更合理的統(tǒng)計(jì)解釋,同時(shí)采用MCMC方法獲得較高的求解效率。Johannesson等[48]給出了聲震融合的貝葉斯MCMC反演框架,Zhang等[25]和Ford等[26]進(jìn)一步發(fā)展了該方法,基于MCMC方法求解源參數(shù),并應(yīng)用于近地面爆源參數(shù)估計(jì)。圖16為貝葉斯MCMC方法求解流程。以上數(shù)據(jù)融合方法基于已有的聲震模型反演爆源參數(shù)和不同場(chǎng)地巖土類型具有不同的聲震模型。數(shù)據(jù)融合方法的反演精度和聲震模型密切相關(guān)。研究表明[42, 45],聲震模型和場(chǎng)地條件的不匹配將導(dǎo)致較大的反演偏差。為減小因聲震模型匹配程度帶來的反演偏差,提高反演模型在不同場(chǎng)地的適用性,Stroujkova等[27]采用人工神經(jīng)網(wǎng)絡(luò)(artificial neural network,ANN)和貝葉斯正則化方法融合小樣本聲震數(shù)據(jù)來估計(jì)當(dāng)量和爆高(埋深),嘗試將機(jī)器學(xué)習(xí)算法應(yīng)用在爆源參數(shù)反演。圖17為基于深度學(xué)習(xí)的爆源參數(shù)反演流程。

(a)Method of relative deviation[22]

(b)Method of logarithmic deviation[23]

(c)Tradeoff curve crossing method[23]

(d)Likelihood function method[24]

(e)Bayesian MCMC method[25, 26, 48]

(f)Machine learning method[27]

圖16 貝葉斯MCMC方法求解流程[25]

圖17 基于深度學(xué)習(xí)的爆源參數(shù)反演流程[27]

4 總結(jié)與展望

經(jīng)過十幾年的研究,發(fā)展了多種近地面聲學(xué)模型和地震波模型,適用沖積土、沙土、沉積巖、花崗巖和土石混合體等多種地質(zhì)條件,聯(lián)合相對(duì)偏差法、折中曲線交叉法和貝葉斯MCMC等數(shù)據(jù)融合算法,或基于聲震數(shù)據(jù)采用機(jī)器學(xué)習(xí)方法,可有效解決多種場(chǎng)地近地面爆源的爆炸當(dāng)量和爆高(埋深)反演問題。聲震分析方法雖得到了快速發(fā)展,但還需重點(diǎn)解決以下問題：

(1)近區(qū)聲/震模型修正完善。受到場(chǎng)地條件限制,不同地質(zhì)、不同爆高和不同距離的聲震數(shù)據(jù)有限,還需更多的試驗(yàn)數(shù)據(jù)對(duì)近地面聲/震模型進(jìn)行補(bǔ)充、修正和完善,從而獲取適用不同地質(zhì)成分的通用模型。目前,聲/震數(shù)據(jù)主要通過外場(chǎng)試驗(yàn)獲取,需進(jìn)一步發(fā)展基于微藥量化爆模擬及超重力離心機(jī)等平臺(tái)的實(shí)驗(yàn)室模擬技術(shù)。

(2)遠(yuǎn)區(qū)聲/震模型快速建模。區(qū)域較大時(shí),受地質(zhì)地形和大氣環(huán)境等因素的影響,近地面爆炸耦合的地震波和聲波表現(xiàn)出復(fù)雜的傳播特性,給遠(yuǎn)區(qū)爆源參數(shù)快速準(zhǔn)確反演帶來了巨大挑戰(zhàn)。目前,聲震分析技術(shù)以近區(qū)為主,需進(jìn)一步發(fā)展大區(qū)域快速建模技術(shù),實(shí)現(xiàn)遠(yuǎn)區(qū)爆源參數(shù)的快速準(zhǔn)確反演。

(3)基于機(jī)器學(xué)習(xí)的快速反演方法。機(jī)器學(xué)習(xí)可有效解決數(shù)據(jù)融合和快速計(jì)算問題,但該方法在聲震分析應(yīng)用較少,需基于機(jī)器學(xué)習(xí)算法,進(jìn)一步發(fā)展小樣本聲震數(shù)據(jù)反演技術(shù)及基于物理機(jī)制的大區(qū)域爆源參數(shù)反演技術(shù)。

(4)未知爆源位置的聲震分析方法。目前聲震分析研究是在爆源位置已知條件下開展的,未考慮爆源位置估計(jì)偏差的影響,需進(jìn)一步開展基于聲震數(shù)據(jù)的爆源位置、當(dāng)量和爆高(埋深)的綜合反演方法研究?？赏ㄟ^貝葉斯MCMC或機(jī)器學(xué)習(xí)方法建立爆源多參數(shù)的綜合反演框架。