何可，郭明，胡章榮，易國財,王仕興

(1.西華師范大學 教育信息技術中心，四川 南充 637002； 2.成都理工大學 地球物理學院,四川 成都 610059)

0 引言

半航空瞬變電磁法(semi-airborne transient electromagnetic，SATEM)采用接地長導線源向地下供給階躍電流，使用無人機搭載接收線圈進行空中觀測，接收二次場[1]。半航空瞬變電磁法與航空瞬變電磁法相比，具有信噪比高、勘探深度大等特點[2]；與地面瞬變電磁比較，具有快速、大面積在復雜地形進行探測的能力，對良導體探測效果較好[3-5]。半航空瞬變電磁法數(shù)據(jù)量大，相較于航空瞬變電磁采用CDI快速成像，由于需考慮線源長度、偏移距等影響，實現(xiàn)難度大，而三維反演由于計算時間較長，目前也很難在實際中應用,所以實踐中數(shù)據(jù)處理仍以一維反演為主。地面瞬變電磁反演研究主要以“煙圈”理論和最優(yōu)化算法為主[6-7]。“煙圈”算法快速、近似，適合作現(xiàn)場解釋工作[8]，但由于不能提供層厚信息，所以只能作為定性研究。最優(yōu)化算法主要以阻尼最小二乘[9-11]和以Occam反演[12-13]引入模型正則化[14-16]的思想為主，阻尼最小二乘法通過改變阻尼因子大小使算法向高斯-牛頓方向或最速下降方向靠近，改進了最小二乘法收斂不穩(wěn)定的特性，該方法對初始模型仍然有較高要求，對復雜模型可能出現(xiàn)結(jié)果不收斂等現(xiàn)象；Occam反演屬于一種L2范數(shù)的正則化反演算法，最為穩(wěn)定，但也主要有兩個方面的不足，一是計算時間較長，每次迭代都需要通過多次正演搜索最佳的正則因子，二是反演結(jié)果過于光滑，對真實電性界面分辨不足。許多學者對正則化因子的選取作出改進并提出一些新方法，如:采用牛頓迭代法和二分法組合[17]降低正則因子搜索的正演次數(shù)；陳小斌等[18]提出MD和CMD兩種正則化因子自適應調(diào)整方案，應用較為廣泛。常規(guī)正則化反演算法正則項通常采用L2范數(shù)，利用L2范數(shù)對電阻率模型進行處理，其假定模型空間分布連續(xù)光滑，反演結(jié)果會弱化對突變電性界面的分辨能力。國外已有學者在重力、地震數(shù)據(jù)三維反演[19-20]中采用L1正則項，國內(nèi)在航空電磁數(shù)據(jù)反演[21]、三維大地電磁反演[22]中，有學者采用L1正則項來提高反演算法對層狀電性介質(zhì)層界面的分辨率。作為L1范數(shù)的一種特例，聚焦反演[23-25]已成熟應用于大地電磁[26]、重力、磁測等[27-28]研究方向，常規(guī)的聚焦反演采用NLCG(非線性共軛梯度)直接求解目標函數(shù)，然而正則項求解不可靠，容易陷入奇異系數(shù)解“陷阱”[22]。

本文正則項采用L1范數(shù)，分析了模型在解空間獲得稀疏解的原因，采用迭代重加權最小二乘法將原問題轉(zhuǎn)化為L2正則化子問題求解，解決L1范數(shù)存在不可導問題，并將兩種算法與Occam反演算法結(jié)果進行比較；采用OpenMP并行策略對雅克比矩陣進行并行運算，提高反演效率；對分段迭代法自適應正則化因子調(diào)整策略進行分析并改進，改進后的自適應正則化因子調(diào)整策略更適合半航空瞬變電磁法，通過與Occam反演對比，L1正則反演結(jié)果的模型分辨率要優(yōu)于Occam反演結(jié)果。

1 半航空瞬變電磁一維正演

長導線源半航空瞬變電磁法以有限長接地導線作為發(fā)射端，不能用電偶極子近似表示，需要沿導線進行積分，對X、Y、Z三個分量均可觀測，目前以垂直地面Z分量為主。假設線源中心為笛卡爾坐標系原點，x軸為長導線方向，z軸向上為正，水平層狀介質(zhì)頻率域Hz分量公式為[27]

(1)

2 L1正則反演理論

2.1 L1正則反演基本原理

常規(guī)L2正則化反演目標函數(shù)可表示為

(2)

式中：Φ(M)L2為總目標函數(shù)；右邊第一項為時間域數(shù)據(jù)擬合項，dobs為觀測數(shù)據(jù)向量，F(xiàn)(M)為理論正演響應，M為模型參數(shù)向量，有n個元素，每個元素均為電阻率自然對數(shù)；Wd為數(shù)據(jù)加權矩陣，一般為主對角矩陣，其元素為觀測數(shù)據(jù)噪聲的倒數(shù)；第二項為L2正則模型約束項,Mref為參考模型，可取CDI成像結(jié)果為參考模型或者以均勻半空間為參考模型，超參數(shù)λ為正則化因子，其控制反演更傾向于擬合數(shù)據(jù)或參考模型，Wm為數(shù)據(jù)加權矩陣，一般可取最小模型約束(單位算子)、最平緩模型約束(梯度算子)、最光滑模型約束(拉普拉斯算子)，本文采用模型約束項為：

L1正則反演目標函數(shù)為：

λ‖Wm[Mref-M]‖1。

(3)

對式(2)和式(3)，可以用如下形式表示：

(4)

(5)

L2范數(shù)是向量各元素平方和開方，具有連續(xù)光滑的特點[21]，這也使得代表層狀介質(zhì)模型各分層電性結(jié)構的最終結(jié)果呈連續(xù)變化，不能明顯刻畫各層界面。L1范數(shù)是指向量各元素mi絕對值之和，也稱為“稀疏規(guī)則算子”，其向量中部分元素與最終輸出沒有任何關系或不提供任何信息，可以允許電性介面發(fā)生陡變。最小化目標函數(shù)的時候考慮mi這些額外特征雖然可以獲得最小誤差，但是也將部分無用信息考慮進去，使得對模型增量的判斷有所偏差，稀疏算子可以去掉無用信息，將其權重置為零。

圖1 數(shù)據(jù)擬合項與模型約束項等值線Fig.1 Contour plots of data fitting items and model constraints

眾所周知，L1范數(shù)取絕對值在存在不可導情況下會使迭代不穩(wěn)定，L1范數(shù)正則項最速下降方向無法得到?？扇∫粋€小值ξ，將式(5)改寫為

隨著網(wǎng)絡時代的發(fā)展，出現(xiàn)了大量老年教育網(wǎng)站，老年人可以通過網(wǎng)絡學習、交友、娛樂等。通過對溫州老年人學習網(wǎng)、夕陽紅·江蘇老年學習網(wǎng)、山東老年大學遠程教育網(wǎng)、老年開放大學等18家國內(nèi)主流老年教育網(wǎng)站調(diào)研，發(fā)現(xiàn)存在課程資源缺乏老年特色，這些網(wǎng)站多數(shù)不能多終端共享，老年遠程教育網(wǎng)站缺乏真正的交互設計等問題。

(6)

針對正則項存在不可導問題，采用迭代重加權最小二乘法(IRLS),將正則項改寫為

(7)

式中：V為對角加權矩陣，主對角線上元素為：

(8)

對式(6)關于增量Δmk(下降方向)求導：

(9)

式(9)右邊第一項為數(shù)據(jù)擬合梯度項，第二項為正則梯度項，Jk為雅可比矩陣(與正演貝塞爾函數(shù)無關)。

2.2 自適應正則因子策略

關于正則化因子的調(diào)整，陳小斌等[18]提出了MD和CMD兩種方案。

MD方案:

(10)

CMD方案:

(11)

>5%，

(12)

則有λk=λk-1，否則λk=λk-1/ω。這里引入了另一個超參數(shù)ω[22]，ω需取適當?shù)闹挡拍苓_到比較好的效果。由于在沒有先驗信息情況下設置參考模型，Mref設為均勻半空間模型時正則因子初始值過大，模型增量較小，所以前幾次迭代模型幾乎不發(fā)生改變。隨著迭代時間和次數(shù)增加，對其作一定改進，使初始正則化因子處于一個較小的值，再按照適當比例逐漸衰減，過程如下。

設定反演結(jié)束條件為：①達到最大迭代次數(shù)；②相鄰兩次擬合差小于給定擬合誤差(Rms1-Rms2

(13)

2.3 OpenMP并行策略

OpenMP是一種共享內(nèi)存系統(tǒng)的多處理器多線程并行語言[36]，采用fork-join(分叉—合并)并行執(zhí)行模式。主線程作串行運算，當遇到并行模塊時，調(diào)用其他從線程構成線程組，同時訪問共享內(nèi)存區(qū)域，執(zhí)行命令，執(zhí)行完畢跳出并行區(qū)域，繼續(xù)執(zhí)行串行命令。OpenMP并行策略具有效率高、執(zhí)行快的特點，適合單機操作。

半航空瞬變電磁反演耗時主要在于雅克比矩陣求解時需要多次調(diào)用正演，采用二維數(shù)組存儲雅可比矩陣，形成二重循環(huán)，主要運算時間也集中在內(nèi)循環(huán)(每一行)參數(shù)的正演計算，參數(shù)越多耗時越長。OpenMP可以對嵌套循環(huán)體內(nèi)多個循環(huán)進行并行運算，本文采用如圖2所示方案，只對內(nèi)循環(huán)進行并行運算策略，對循環(huán)內(nèi)變量作無關處理，各線程對數(shù)據(jù)、函數(shù)的調(diào)用相對獨立。本文用3層模型進行試算，觀測值15道，初始模型為30層均勻半空間；計算機處理器為Intel(R) Core(MT)i5-8265U,主頻 1.6 GHz，4核8線程，分別進行5次、15次、30次、60次迭代，計算時間對比如表1所示，采用并行運算后，運算時間約縮短3/4,效率顯著提高。

圖2 OpenMP運算并行策略Fig.2 OpenMP computing parallel strategy

表1 串并行計算時間對比

3 理論模型反演分析

設定半航空瞬變電磁參數(shù)如下：線源長度1 km，電流20 A,線圈接收高度20 m，線圈面積為1，偏移距為250 m，即測點坐標為(0，250，20)；設初始模型30層，初始層厚為2 m，依次遞增，每層電阻率為50 Ω·m，對理論模型響應值均加入3%的高斯白噪聲，并與Occam反演結(jié)果相比較。

3.1 三層模型反演算例

3.1.1 H型模型反演

H型模型初始電阻率由淺至深分別為100、10、100 Ω·m，對應層厚分別為50、50 m，最后一層層厚設為無窮大。其不同方法反演結(jié)果見圖3。

由圖3a可見，Occam反演結(jié)果較為光滑，對分界面刻畫較弱，第一、二層電阻率值與真實模型對比有所偏離。圖3b中，L1正則反演結(jié)果對各層分界面刻畫清晰，電阻率值和層厚更接近真實模型。圖3c顯示了反演算法擬合差和正則化因子衰減過程，迭代較為平穩(wěn)。圖3d為觀測數(shù)據(jù)與反演擬合響應的對比，數(shù)據(jù)擬合較好，接近真實值。圖3e中可以看出第5次迭代時趨向于L2范數(shù)反演結(jié)果，隨著迭代次數(shù)增加，對角加權矩陣V變得更復雜，層界面刻畫逐漸清晰。第15次迭代時電阻率值接近真實模型，層界面刻畫較為準確，第35次迭代時達到擬合條件，與真實模型接近。

圖3 H形模型反演結(jié)果Fig.3 Inversion results of H-shaped model

3.1.2 K型模型反演

對K型模型，電阻率設為10、100、10 Ω·m，對應層厚分別為50、100 m，最后一層設為無窮大半空間。圖4a、圖4b顯示Occam和L1范數(shù)這兩種算法都能有效反映出高阻部分，但Occam反演在層厚較大時對層界面的刻畫更為光滑，第一層和第三層的電阻率值與真實模型比較出現(xiàn)較小偏差，而L1范數(shù)對層界面刻畫清晰， 電阻率值更接近真實模型。從圖4e可以看出，第20次迭代時基本能反映出層狀介質(zhì)的基本模型；通過和H型反演迭代次數(shù)比較，可以看出對高阻的擬合相對困難，迭代次數(shù)增加。

圖4 K形模型反演結(jié)果Fig.4 Inversion results of K-shaped model

3.2 四層HK型模型反演算例

四層層狀模型電阻率設為100、10、100、10 Ω·m,對應層厚分別為50、50、50 m，最后一層設為無窮。如圖5所示，總體結(jié)果與圖4結(jié)論相似，可以看出隨模型的復雜程度增加，迭代次數(shù)也隨之增加。本次擬合2種反演方法都迭代在60次左右，Occam反演低阻擬合較好，對第三層的高阻部分擬合較差；L1正則反演對層界面位置反應與真實模型吻合較好，同時電阻率值也接近真實模型，說明L1正則反演方法相對于Occam反演對高阻的反應相對靈敏。隨著模型的復雜度變高，反演迭代次數(shù)也隨之增加,正則化因子較小時約束減弱，擬合差也有波動現(xiàn)象。然而，實際工作過程中并不知道地下實際電性結(jié)構分布，但L1正則反演在一定迭代次數(shù)后可以反映出層狀介質(zhì)大致模型，所以可以根據(jù)實際地質(zhì)情況調(diào)整反演約束條件。

圖5 HK形模型反演結(jié)果Fig.5 Inversion results of HK-shaped model

4 總結(jié)

針對半航空瞬變電磁L2正則反演層界面刻畫不夠清晰的問題，提出了利用L1正則反演處理半航空瞬變電磁數(shù)據(jù)，對自適應正則化因子調(diào)整策略進行改進優(yōu)化，并用三層的H、K地電模型和四層的HK型地電模型對兩種反演算法進行了試算與比較，通過分析得到以下結(jié)論：

1)基于L1范數(shù)的正則化反演算法比L2更有利于得到稀疏解，允許模型電阻率在空間分布上存在尖銳；對層界面位置反映較好，相對復雜的層狀介質(zhì)模型也能得到較好的反演效果；對高阻的分辨率也較高。

2)L1范數(shù)正則項采用迭代重加權最小二乘法，解決了不可導問題，反演結(jié)果早期趨向于L2范數(shù)反演結(jié)果，隨著迭代次數(shù)增加，對層界面的分辨逐漸清晰；正則化因子分段迭代法中比例因子的選取會影響收斂速度，因子越大迭代越穩(wěn)定，但收斂較慢。

3)正則化因子采用分段迭帶法，和OpenMP并行策略使得反演較快，然而相較于Occam通過搜索方式得到的正則化因子不是最優(yōu)的，面對復雜模型的反演時可能會出現(xiàn)假異常。