井下采空區時域不連續伽遼金法電磁正演

王樹奇，黃茸茸，吳 楠

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

0 引 言

煤炭是中國能源支柱產業，然而大范圍的煤炭開采在礦山中留下了大量的采空區，導致地下巖體原有的力學平衡被打破，扭曲的巖體隨時可能發生位移、巖爆等事故；更嚴重的是采空區會被瓦斯、地下水等充填，給礦區的開采作業埋下重大的安全隱患，采用電磁數值計算對礦井復雜地質結構進行探測分析具有重要的意義。

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

時域不連續伽遼金方法(discontinuous galerkin time domain,DGTD)[7-12]繼承了FETD 采用非結構網格對復雜幾何外形擬合好，便于使用高階基函數的優點；該算法對Maxwell 方程采用伽遼金加權法得到弱解形式，放松了單元之間的邊界條件，單元之間通過數值通量聯系交換數據；高階疊層矢量基函數的引入是一個突破，它滿足單元邊界上的 場切向連續性，消除了節點法中存在的“偽解”并且可顯著提高該方法的計算精度；在時間離散方面既支持顯式時間格式也支持隱式時間格式；在空間離散方面由于單元間的高度獨立性，對于每個單元可以選用不同階的基函數進行離散，需要高精度時采用高階基函數，對精度要求低時采用低階基函數。可以將大型多尺度問題被拆分成為一個個相對獨立的小問題，使得分析問題更為簡單方便，以上原因更加說明了研究DGTD算法的價值與意義。

文中首次將基于高階疊層矢量基函數的DGTD方法應用于井下采空區三維正演研究，建立了井下采空區三維空間模型，分析了疊層矢量基函數階數對正演結果的影響，試驗驗證了井下采空區模型正演的計算精度。

1 DGTD方法

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

(1)

式中ε為介電系數；μ為磁導系數；σe為電導率；σm為導磁率；E和H為電場強度與磁場強度。

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

(2)

φq′為加權函數，使用以下恒等式對公式(2)中的×E和×H項進行代數變換

(3)

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

對上述數學關系式 y=134.1x2-182.8x+385.3進行一階求導，確定出在定義域范圍內日耗電最低時的抽油機懸點載荷利用率為68.2%，稱之為第一抽油機懸點載荷利用率。按照類似的方法對所得的40組數據分別進行分析，得出40組在各自電動機負載率下的最優抽油機懸點載荷利用率，具體數據如表1。

(4)

圖1 數值通量Fig.1 Numerical flux

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

(5)

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

(6)

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

(7)

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

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

(8)

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

(9)

(10)

其中系數為

(11)

2 高階疊層矢量基函數

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

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

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

表1 疊層矢量基函數Table 1 Hierarchical vector basis functions

3 矩形諧振腔模擬計算

為了分析不同高階疊層矢量基函數的DGTD方法的計算精度，建立了尺寸為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個均勻網格尺寸的四面體，圖3所示計算域被劃分為2 292個網格大小不同的四面體單元。在點(0，0，0)m處激發一個中心頻率為310 MHz，極化方向為(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階以及混合階疊層矢量基函數時，DGTD方法在觀測點(0.7，0，0)m處得到的時域波形分別如圖4(a)～圖4(d)所示。

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

圖4 基于疊層基函數的DGTD方法接收點數值結果Fig.4 Numerical results of the receiving point of the DGTD method based on the hierarchical vector basis functions

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

4 井下采空區電磁正演

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

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

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

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

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

式中E為時域波形數值解；Eref為仿真軟件計算結果。

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

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

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

5 結 論

1)利用基于高階疊層矢量基函數的DGTD方法對井下采空區地質結構進行三維電磁正演，解決了傳統時域電磁方法在模擬曲線邊界時存在階梯效應，對復雜、曲面目標不能準確建模、不支持高階基函數、隱式離散以及“偽解”等問題。

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