王樹奇，黃茸茸，吳 楠

(1.西安科技大學(xué) 通信與信息工程學(xué)院，陜西 西安 710054;2.陜西學(xué)前師范學(xué)院 信息工程學(xué)院，陜西 西安 710100)

0 引 言

煤炭是中國能源支柱產(chǎn)業(yè)，然而大范圍的煤炭開采在礦山中留下了大量的采空區(qū)，導(dǎo)致地下巖體原有的力學(xué)平衡被打破，扭曲的巖體隨時(shí)可能發(fā)生位移、巖爆等事故；更嚴(yán)重的是采空區(qū)會被瓦斯、地下水等充填，給礦區(qū)的開采作業(yè)埋下重大的安全隱患，采用電磁數(shù)值計(jì)算對礦井復(fù)雜地質(zhì)結(jié)構(gòu)進(jìn)行探測分析具有重要的意義。

國內(nèi)外學(xué)者針對時(shí)域電磁正演的方法做了大量研究，但對于井下采空區(qū)的時(shí)域電磁正演的研究多采用時(shí)域有限差分(finite difference time domain，F(xiàn)DTD)方法和時(shí)域矢量有限元(finite element time domain，F(xiàn)ETD)[1-3]。于景邨和常江浩采用三維有限差分法對老空水的礦井瞬變電磁響應(yīng)進(jìn)行了研究[4]。楊道學(xué)提出了基于卷積完全匹配層的交替方向隱式的FDTD方法，并將其應(yīng)用于探地雷達(dá)正演[5]。FDTD方法具有理論難度低，計(jì)算效率高等優(yōu)點(diǎn)，但是該方法中使用的六面體剖分單元在復(fù)雜幾何結(jié)構(gòu)中具有局限性，且不支持高階基函數(shù)等問題將直接影響計(jì)算精度。張永超和李宏杰等研究了時(shí)域矢量有限元三維正演[6]，拓展了礦井瞬變電磁正演對復(fù)雜地質(zhì)模型的適用性，棱邊矢量基函數(shù)的使用很好地解決了節(jié)點(diǎn)方法中的“偽解”，但只是用了低階棱邊基函數(shù)且隱式離散格式，計(jì)算效率低。

時(shí)域不連續(xù)伽遼金方法(discontinuous galerkin time domain,DGTD)[7-12]繼承了FETD 采用非結(jié)構(gòu)網(wǎng)格對復(fù)雜幾何外形擬合好，便于使用高階基函數(shù)的優(yōu)點(diǎn)；該算法對Maxwell 方程采用伽遼金加權(quán)法得到弱解形式，放松了單元之間的邊界條件，單元之間通過數(shù)值通量聯(lián)系交換數(shù)據(jù)；高階疊層矢量基函數(shù)的引入是一個(gè)突破，它滿足單元邊界上的 場切向連續(xù)性，消除了節(jié)點(diǎn)法中存在的“偽解”并且可顯著提高該方法的計(jì)算精度；在時(shí)間離散方面既支持顯式時(shí)間格式也支持隱式時(shí)間格式；在空間離散方面由于單元間的高度獨(dú)立性，對于每個(gè)單元可以選用不同階的基函數(shù)進(jìn)行離散，需要高精度時(shí)采用高階基函數(shù)，對精度要求低時(shí)采用低階基函數(shù)?？梢詫⒋笮投喑叨葐栴}被拆分成為一個(gè)個(gè)相對獨(dú)立的小問題，使得分析問題更為簡單方便，以上原因更加說明了研究DGTD算法的價(jià)值與意義。

文中首次將基于高階疊層矢量基函數(shù)的DGTD方法應(yīng)用于井下采空區(qū)三維正演研究，建立了井下采空區(qū)三維空間模型，分析了疊層矢量基函數(shù)階數(shù)對正演結(jié)果的影響，試驗(yàn)驗(yàn)證了井下采空區(qū)模型正演的計(jì)算精度。

1 DGTD方法

三維非均勻各向同性介質(zhì)Maxwell旋度方程為

(1)

式中ε為介電系數(shù)；μ為磁導(dǎo)系數(shù)；σe為電導(dǎo)率；σm為導(dǎo)磁率；E和H為電場強(qiáng)度與磁場強(qiáng)度。

對公式(1)采用伽遼金加權(quán)法，在四面體單元內(nèi)進(jìn)行體積分，可以得到以下弱解形式

(2)

φq′為加權(quán)函數(shù)，使用以下恒等式對公式(2)中的×E和×H項(xiàng)進(jìn)行代數(shù)變換

(3)

在相鄰區(qū)域單元交界面上采用數(shù)值通量[13-14]的形式來保證單元之間切向場的連續(xù)性，數(shù)值通量格式如下

對上述數(shù)學(xué)關(guān)系式 y=134.1x2-182.8x+385.3進(jìn)行一階求導(dǎo)，確定出在定義域范圍內(nèi)日耗電最低時(shí)的抽油機(jī)懸點(diǎn)載荷利用率為68.2%，稱之為第一抽油機(jī)懸點(diǎn)載荷利用率。按照類似的方法對所得的40組數(shù)據(jù)分別進(jìn)行分析，得出40組在各自電動機(jī)負(fù)載率下的最優(yōu)抽油機(jī)懸點(diǎn)載荷利用率，具體數(shù)據(jù)如表1。

(4)

圖1 數(shù)值通量Fig.1 Numerical flux

將數(shù)值通量代入公式(3)中可得

(5)

公式(5)的空間部分采用基函數(shù)在每個(gè)單元中展開

(6)

整理簡化后得到基于四面體單元的DGTD顯式半離散方程組

(7)

在得到空間離散格式之后將進(jìn)行時(shí)間離散，最后根據(jù)時(shí)間離散格式完成時(shí)間迭代。可采取不同方式對得到的空間離散方程進(jìn)行時(shí)間偏導(dǎo)數(shù)的離散，若采用隱式時(shí)間離散方案則將降低計(jì)算效率，文中采用已經(jīng)廣泛應(yīng)用于時(shí)域算法中的蛙躍算法進(jìn)行時(shí)間離散[15]，蛙躍算法是一種顯式時(shí)域離散方法，電磁場相隔半個(gè)時(shí)間步交替進(jìn)行迭代計(jì)算。規(guī)定電場E在整數(shù)時(shí)間步tn采樣，磁場H在半整數(shù)步tn+1/2進(jìn)行采樣。

在公式(7)中對時(shí)間的導(dǎo)數(shù)采取中心差分的格式處理，由二階中心差分代替一階時(shí)間導(dǎo)數(shù)[16]

(8)

由于在迭代過程中會出現(xiàn)電場的半整數(shù)步、磁場的整數(shù)步，對其取平均近似[16]有以下形式

(9)

(10)

其中系數(shù)為

(11)

2 高階疊層矢量基函數(shù)

為了克服低階矢量基函數(shù)精度較低的缺點(diǎn)，NOTAROS等提出了高階插值型矢量基函數(shù)[17-19]。高階插值型矢量基函數(shù)具有線性獨(dú)立性好、物理解釋明確、編程實(shí)現(xiàn)方便等優(yōu)點(diǎn)。但在一個(gè)四面體內(nèi)插值型矢量基函數(shù)具有更多的未知量，為了滿足交界面處切向場連續(xù)，插值矢量基不允許不同階次基函數(shù)的混合，這導(dǎo)致了在原本低階單元即可足夠精確描述的區(qū)域引入了大量冗余自由度，從而增加了存儲量和計(jì)算時(shí)間，降低了算法分析的效率。

高階疊層型矢量基函數(shù)[20-21]的提出正是為了解決高階插值型矢量基函數(shù)的這些缺點(diǎn)。首先，高階疊層型基函數(shù)具有疊層嵌套的特點(diǎn)，高階疊層構(gòu)矢量基函數(shù)的構(gòu)建是基于低階基函數(shù)逐階遞增的，即高階的疊層基函數(shù)里面包含有階數(shù)較低的基函數(shù)，這在保證高精度的情況下降低了基函數(shù)構(gòu)造難度以及存儲量。其次，借助高階疊層矢量基函數(shù)，時(shí)域不連續(xù)伽遼金算法可以實(shí)現(xiàn)計(jì)算域內(nèi)高階和低階單元的混合計(jì)算，數(shù)值結(jié)果表明高階疊層矢量基函數(shù)的引用使得對目標(biāo)離散時(shí)采用較大尺寸單元得到的結(jié)果與采用較小尺寸單元剖分的低階基函數(shù)方案有相當(dāng)?shù)木?。疊層矢量基函數(shù)是一種基于棱邊、面和體積的基函數(shù)，它將自由度賦予棱邊、面和體積而不是單元節(jié)點(diǎn)。它隱含了散度為零的邊界條件，消除了偽解，非常適合用來表示矢量場。

綜上可知,高階疊層矢量基函數(shù)在處理大規(guī)模電磁問題和多尺度電磁問題上都具有非常大的優(yōu)勢，具有很高的實(shí)用性與靈活性。表1給出了文中所采用的基于四面體單元的0.5階到2.5階疊層矢量基函數(shù)形式。當(dāng)然，疊層型矢量基函數(shù)還有更高的階的形式，但隨著基函數(shù)階數(shù)的增加，自由度會以數(shù)倍增加，這將增加計(jì)算難度和時(shí)間成本。并且當(dāng)基函數(shù)階數(shù)為2.5階時(shí)，該算法就可以獲得較好的計(jì)算精度。因此文中就0.5階至2.5階基函數(shù)展開研究。

表1 疊層矢量基函數(shù)Table 1 Hierarchical vector basis functions

3 矩形諧振腔模擬計(jì)算

為了分析不同高階疊層矢量基函數(shù)的DGTD方法的計(jì)算精度，建立了尺寸為2.0 m×0.3 m×0.2 m的無損耗矩形PEC邊界諧振腔模型，如圖2、圖3所示。

圖2 矩形PEC邊界均勻單元尺寸諧振腔Fig.2 Uniformly meshed rectangular resonator with PEC boundary

圖2中腔體被劃分為2 261個(gè)均勻網(wǎng)格尺寸的四面體，圖3所示計(jì)算域被劃分為2 292個(gè)網(wǎng)格大小不同的四面體單元。在點(diǎn)(0，0，0)m處激發(fā)一個(gè)中心頻率為310 MHz，極化方向?yàn)?1，1，1)的布萊克曼-哈里斯(blackman harris window，BHW)電偶極子源。

圖3 矩形PEC邊界不同單元尺寸邊界諧振腔Fig.3 Non-uniformly meshed rectangular resonator with PEC boundary

分別選取0.5階、1.5階、2.5階以及混合階疊層矢量基函數(shù)時(shí)，DGTD方法在觀測點(diǎn)(0.7，0，0)m處得到的時(shí)域波形分別如圖4(a)～圖4(d)所示。

從圖4可以看出,當(dāng)疊層基函數(shù)為0.5階時(shí)，其結(jié)果的計(jì)算精度較差。隨著高階疊層矢量基函數(shù)的階數(shù)增加，數(shù)值計(jì)算結(jié)果的精度有了顯著提升，經(jīng)分析見表2所示結(jié)果。當(dāng)基函數(shù)階數(shù)分別為0.5階、1.5階、2.5階以及混合階疊層矢量基函數(shù)時(shí)，計(jì)算結(jié)果的相對誤差為43.94%、22.44%、2.29%和2.42%。從0.5階到2.5階，DGTD方法的計(jì)算相對誤差降低了41.65%，計(jì)算精度得到了明顯的改善；混合階基函數(shù)的計(jì)算中基函數(shù)的階數(shù)根據(jù)四面體單元的尺寸進(jìn)行選擇，其結(jié)果的相對誤差與2.5階基函數(shù)結(jié)果相近，該結(jié)果表明基于高階疊層矢量基函數(shù)的DGTD方法在保證計(jì)算精度的同時(shí)可實(shí)現(xiàn)計(jì)算域內(nèi)高階和低階單元的混合計(jì)算，對復(fù)雜目標(biāo)的電磁特性研究具有很高的實(shí)用性。

圖4 基于疊層基函數(shù)的DGTD方法接收點(diǎn)數(shù)值結(jié)果Fig.4 Numerical results of the receiving point of the DGTD method based on the hierarchical vector basis functions

表2 矩形諧振腔模擬計(jì)算誤差分析Table 2 Error analysis of simulation of rectangular resonator

4 井下采空區(qū)電磁正演

建立煤層區(qū)域?yàn)?0 m×10 m×10 m，井下采空區(qū)尺寸為4 m×4 m×2 m的三維正演地質(zhì)模型，如圖5，圖6所示。背景設(shè)置為均勻介質(zhì)，取相對介電常數(shù)為4.09褐煤進(jìn)行模擬驗(yàn)證。以空間模型的中心位置作為原點(diǎn)，掘進(jìn)工作面位于z方向上。在(0 m，0 m，4 m)處設(shè)置中心頻率為31 MHz的BHW電偶極子激勵源。接收點(diǎn)放置于掘進(jìn)工作面上，其位置坐標(biāo)為(0 m，1 m，5 m)。

圖5 井下采空區(qū)三維模型Fig.5 Three-dimensional model of underground goaf

圖6 井下采空區(qū)三維模型剖分Fig.6 Meshes of three-dimensional model of underground goaf

為了研究所選取的基函數(shù)對井下電磁響應(yīng)特性的影響，設(shè)計(jì)了地質(zhì)模型，分別將0.5階、1.5階、2.5階疊層矢量基函數(shù)應(yīng)用于電磁正演中，結(jié)果如圖7所示。對比不同階基函數(shù)對采空區(qū)的電磁響應(yīng)特征的影響，當(dāng)基函數(shù)為0.5階、1.5階、2.5階3種不同階數(shù)時(shí)均可分辨反射波信號。從圖7所示的基于DGTD方法井下采空區(qū)電場響應(yīng)結(jié)果可看出，隨著基函數(shù)階數(shù)的增加，直達(dá)波和反射波的幅值精確度有了顯著提高，且采空區(qū)反射波曲線拖延現(xiàn)象得到了顯著改善。對比DGTD方法中不同階疊層基函數(shù)正演波形，誤差分析結(jié)果見表3。計(jì)算中將接收點(diǎn)的場值記錄并與仿真軟件結(jié)果進(jìn)行比較，其相對誤差定義為

Error=20log(|E-Eref|/Eref max)(dB)

式中E為時(shí)域波形數(shù)值解；Eref為仿真軟件計(jì)算結(jié)果。

經(jīng)分析得：井下采空區(qū)正演的相對誤差從0.5階的-19.12 dB縮小至2.5階的-73.19 dB，計(jì)算精度提高了74%。由以上分析可得高階疊層矢量基函數(shù)的階數(shù)將直接影響DGTD方法的正演計(jì)算精度，采空區(qū)正演時(shí)域結(jié)果隨基函數(shù)階數(shù)的增加，幅值的精確度大幅提升。

圖7 基于DGTD方法煤層采空區(qū)正演結(jié)果Fig.7 Forward modeling results of coal seam goaf based on DGTD method

表3 井下采空區(qū)電磁正演誤差分析Table 3 Error analysis of electromagnetic forward modeling results in underground goaf

5 結(jié) 論

1)利用基于高階疊層矢量基函數(shù)的DGTD方法對井下采空區(qū)地質(zhì)結(jié)構(gòu)進(jìn)行三維電磁正演，解決了傳統(tǒng)時(shí)域電磁方法在模擬曲線邊界時(shí)存在階梯效應(yīng)，對復(fù)雜、曲面目標(biāo)不能準(zhǔn)確建模、不支持高階基函數(shù)、隱式離散以及“偽解”等問題。

2)高階疊層矢量基函數(shù)的引入提高了DGTD方法的計(jì)算精度，隨著高階疊層矢量基函數(shù)階數(shù)由0.5階、1.5階到2.5階增加，基于DGTD方法的井下采空區(qū)電磁正演計(jì)算結(jié)果的精確度分別提高了69%和74%；當(dāng)基函數(shù)為2.5階時(shí)，基于DGTD方法的井下采空區(qū)電磁正演計(jì)算結(jié)果可精準(zhǔn)得出直達(dá)波絕對幅值的大小、反射波絕對幅值的大小、反射波形的拖尾程度等表征，達(dá)到較好的電磁正演的效果。該方法為井下采空區(qū)的電磁正演提供了一種新的思路。