電阻抗成像探測地下溶洞的幾種方法*

吳澤宇，劉祚秋，汪利

中山大學航空航天學院，廣東 廣州 510006

電阻抗成像方法是可以用于探測地下溶洞的一種方法，其原理是根據物體內部不同物質的導電參數（如電阻率、電容率）的不同，通過對物體表面電流、電壓的施加和測量來獲知物體內部導電參數的分布，進而重建出反映物體內部結構的圖像，該方法的難點一是固有的病態性，即邊界電壓數據的微小擾動可能引起解的巨大變化；另一個困難是它的信息量小。如何實現系統的高精度、高分辨率和算法的快速收斂是其應用的主要難點［1］。

電阻抗成像的起源可以追溯到20 世紀20 年代，地球物理學研究者提出了線性電極陣列的電阻率成像（Resistivity Imaging）技術。19 世紀70年代，生物醫學研究者提出了圓形電極陣列的斷層電阻率測量技術（Tomographic Resistivity Measurement Technique）。1978 年，Henderson 和Webster 做出了第一幅電阻抗圖像。他們使用固定于胸部上的由一大電極和與之相對的若干小電極組成的電極系統。通過測量從各小電極流向大電極的電流所形成的等位差，獲得了可以清楚地顯示肺臟位置的阻抗圖像。但是這還不是斷層圖像，而是類似X 胸片的透視圖像。1982 年，英國Sheffield大學的Brown 和Barber 實現了第一個手臂的阻抗層析圖像，開辟了電阻抗層析成像（EIT）技術這一新的研究領域［2］。

本文采用電阻抗成像法探測了地下溶洞，并采用Tikhonov正則化和稀疏正則化方法以解決反問題固有的非適定性問題；在正則化的過程中運用了靈敏度分析法，又分別采用線性近似直接求解和迭代求解兩種方法進行求解，并進一步分析了幾種方法的識別精度與誤差敏感度，為實際工程中的地下溶洞識別提供了一定的參考。

1 問題的描述

1.1 正問題

考慮一個空間區域Ω 及其邊界?Ω［3-4］，電阻抗成像的正問題是在已知材料電阻抗率ρ(x)，x∈Ω（或其倒數，視電阻率c(x) = 1/ρ(x)）的情況下，求解區域內的電勢u(x)（或電壓）分布。考慮邊界?Ω上輸入輸出電流，正問題的控制方程［5］為

其中n為邊界上的外法向向量，SI為所有邊界電流點（或電極）的集合，考慮電流幾乎為單點輸入輸出，el表示的是包含電流點l的微小區域，Il是點l處的電流量。Il＞0 表示電流輸入，Il＜0 則為電流輸出。如果某個電極方案只有點1 輸入單位電流，而點2 輸出單位電流，則有I1= 1，I2= -1，Il=0，l∈SI{1，2}。此外，還可以測量一系列邊界點的電勢，這些測量的電勢點的集合記為Su；電勢點Su與電流點SI可以有部分重合。問題（1）～（2）是典型的Neumann 邊值問題，為了確保解的存在性，需要滿足所有輸入輸出電流總和為0，即

其中uk表示點k上的電勢。

基于以上描述，考慮電阻抗成像探測地下溶洞這一問題（如圖1 所示）。問題所屬區域簡化為圖1右側的矩形區域。電流點與電勢點不一致，但數目相等。

圖1 電阻抗成像探測地下溶洞示意圖Fig.1 Diagram of electrical impedance imaging for cavity detection

通?？梢允褂糜邢拊M行正問題求解，得到如下方程

其中Kj為單位電阻率時的單元剛度矩陣。

1.2 反問題

2 電阻抗成像的理論與方法

2.1 初值或平均視電阻率估計

此時，可以選擇初值為c0=c01。

2.2 靈敏度分析

接下來，將根據已有的初值c0和靈敏度矩陣進行參數識別?？紤]使用不同的正則化方法（Tikhonov 正則化［8］和稀疏正則化［9］）以及是否迭代的因素，主要可以分為4 種方法：Tikhonov 更新法、稀疏更新法、Tikhonov迭代法、稀疏迭代法。

2.3 Tikhonov更新法

這里假設初值c0與待識別結果c比較接近，此時可以對R(c)使用Taylor展開，得到

其中正則化稀疏λ＞0 可通過L-curve 方法［10-11］進行確定。

2.4 稀疏更新法

同樣考慮如公式（12）～（13）的線性近似的過程，但使用稀疏正則化處理非適定性。實際溶洞探測問題中，溶洞的存在相對于探測區域通常體積較小，在空間分布上具有一定稀疏性（或更新量(c-c0)是稀疏的），因此可以使用稀疏正則化進行處理。稀疏更新法的求解問題為

2.5 Tikhonov迭代法

2.6 稀疏迭代法

與Tikhonov迭代法類似，稀疏迭代法則表現為使用稀疏更新法重復、迭代求解的過程。迭代式為

3 數值算例

選取矩形區域長寬均為4 m×4 m，每邊均勻劃分為8 個單元，假設初始區域的視電阻率是均勻的，取值為c0= 10。考慮如圖2 中三種溶洞情況（溶洞處的視電阻率等效為3，相對于土體的10，有較大的折減），其中Case 1、Case 2 為單個溶洞的情況、Case 3為多個溶洞的情況，圖塊顏色與視電阻率的對應關系如圖右側色標所示。

圖2 溶洞探測的三種情況Fig.2 Three cases of cavity detection

考慮如圖1右側的矩形區域和電極方案，通過給定的電流與真實視電阻率分布得到電勢結果，分別利用Tikhonov 更新法、稀疏更新法、Tikhonov迭代法和稀疏迭代法識別以上三種溶洞情況。其中Tikhonov 法的正則化參數通過L-curve 方法進行確定，而稀疏正則化參數取10-8。

另外，考慮誤差存在時的情形，此時電勢情況u=u0(1+ fr·Randn)；其中，u0為真實電勢；Randn是與電勢矩陣同階的隨機數矩陣，其每個元素都是[-1，1]上的正態分布隨機變量；fr為變量的系數，這里取fr = 0.001。

Case 1在無誤差和有誤差情況下的識別結果分別如圖3～4 所示，Case 2 在無誤差和有誤差情況下的識別結果分別如圖5～6 所示，Case 3 在無誤差和有誤差情況下的識別結果分別如圖7～8所示。從結果中可以看出：

圖3 無誤差時Case 1的識別結果Fig.3 Results of Case 1 without error

1）Tikhonov 更新法在無誤差的情況下識別單個溶洞（圖3 和圖5）時，在溶洞位置可以識別到視電阻率的折減，但與真實情況相比，視電阻率在數值上相差較大，結果較為粗糙；而在有誤差的情況下（圖4和圖6），在溶洞位置也能看出視電阻率的差別，但與無誤差情況相比，結果的差別較大；而在識別多個溶洞的情況時（圖7 和圖8），無論是否有誤差，視電阻率的分布較為混亂，多處出現電阻率折減，電阻率在數值上與真實情況有較大差別，無法識別出正確的溶洞位置。

圖4 有誤差時Case 1的識別結果Fig.4 Results of Case 1 with error

圖5 無誤差時Case 2的識別結果Fig.5 Results of Case 2 without error

圖6 有誤差時Case 2的識別結果Fig.6 Results of Case 2 with error

2）稀疏更新法在無誤差的情況下識別單個溶洞（圖3 到圖6）時，可以明顯地看到溶洞位置有明顯的視電阻率折減，而在其他位置的視電阻率在數值上與真實情況十分接近；識別相鄰的多個溶洞的情況時（圖7和圖8），較大溶洞處有視電阻率的折減，但不能完全識別出溶洞形狀，數值上也與真實情況有較大差別；且較小溶洞處沒有視電阻率的折減，其他位置的視電阻率也與真實情況有一定差別，但各處視電阻率之間相差不大，且無誤差和有誤差結果相似。

圖7 無誤差時Case 3的識別結果Fig.7 Results of Case 3 without error

圖8 有誤差時Case 3的識別結果Fig.8 Results of Case 3 with error

3） Tikhonov 迭代法在無誤差的情況下識別單個溶洞（圖3 和圖5）時，在溶洞位置可以識別到視電阻率的折減，其他位置的視電阻率與真實情況在數值上存在較大差別，且視電阻率分布不均勻。在有誤差的情況下（圖4和圖6），視電阻率分布混亂，有多處位置存在明顯的視電阻率折減，無法識別出溶洞位置；而在無誤差的情況下識別多個溶洞時（圖7），能準確識別到溶洞位置處的視電阻率折減，其他位置的視電阻率也較為均勻。但在有誤差的情況下（圖8），整個識別區域的視電阻率分布混亂，無法識別到溶洞位置。

4） 稀疏迭代法溶洞探測效果非常好，在無誤差的情況下識別單個與多個溶洞（圖3、圖5 和圖7）時，結果與真實情況幾乎完全一致；而在有誤差的情況下（圖4、圖6 和圖8），結果也基本保持不變。

4 結 論

Tikhonov 更新法方法簡單，可以大體識別單個溶洞的位置，但結果粗糙，對誤差較為敏感，且無法識別相鄰的多個溶洞；稀疏更新法也比較簡單，識別單個溶洞時結果較好，對誤差敏感度也較低，但在識別相鄰的多個溶洞時存在遺漏；Tikhonov 迭代法比Tikhonov 更新法復雜，在無誤差時能識別單個或多個溶洞的位置，但對誤差的敏感程度非常大，在有誤差情況下完全無法識別出溶洞位置；稀疏迭代法計算量大，但識別結果非常好，對誤差的敏感度也很低，無論是否有誤差，無論是識別單個還是多個溶洞，都能得到非常理想的結果。

總而言之，從結果上看，幾種方法的排序（由好到壞）為：稀疏迭代法＞稀疏更新法＞Tikhonov更新法＞Tikhonov迭代法。