基于FCM 聚類模型約束的二維初至旅行時(shí)反演

劉佳成，張志勇*，周欽淵，李曼，李紅立

（1. 東華理工大學(xué)地球物理與測(cè)控技術(shù)學(xué)院，江西南昌 330013；2. 湖南省交通規(guī)劃勘察設(shè)計(jì)研究院有限公司，湖南長(zhǎng)沙 410008；3. 廣西科技大學(xué)，廣西柳州 545006）

0 引言

地震旅行時(shí)成像是以地震波在不同介質(zhì)中傳播規(guī)律為基礎(chǔ)，在地表或井中觀測(cè)地震信號(hào)，通過旅行時(shí)反演重構(gòu)地質(zhì)體內(nèi)部速度分布的一種地球物理勘探方法。該方法廣泛應(yīng)用于工程檢測(cè)、油氣勘察、地球內(nèi)部結(jié)構(gòu)研究等［1-3］。地震初至波旅行時(shí)成像能識(shí)別、追蹤目標(biāo)地質(zhì)體，穩(wěn)定性較好，對(duì)淺層成像精度較高。

地震旅行時(shí)成像的過程主要包括正演和反演兩部分。早期的地震旅行時(shí)正演采用經(jīng)典射線追蹤方法，如打靶法［4］、高斯射線束法［5］、彎曲射線法［6］，但可能存在多條射線路徑，反演極易出現(xiàn)多解。為此，后續(xù)發(fā)展了基于程函方程的地震旅行時(shí)算法，如有限差分算法［7］、快速推進(jìn)算法（Fast Marching Method，F(xiàn)MM）［8］、最短路徑算法［9］等。其中FMM 是一種網(wǎng)格數(shù)值算法，通過有限差分或有限元方法求解程函方程追蹤地震波前界面，具有較高的精度和效率，因而被廣泛應(yīng)用［10］。由于地震資料觀測(cè)的局限性、觀測(cè)數(shù)據(jù)中存在誤差等原因，地震旅行時(shí)反演是一個(gè)不適定問題，反演算法主要有阻尼最小二乘法［11］、基于正交分解的最小二乘法（LSQR）［12］、子空間反演法［13］、代數(shù)重建技術(shù)（ART）［14］等，這些算法在不同程度上能降低反演問題的不適定性。此外，正則化技術(shù)［15］通過引入地球物理先驗(yàn)信息減少了反問題的多解性、提高了解的穩(wěn)定性，因此被廣泛用于地球物理反演［16］。最小結(jié)構(gòu)模型約束是將最小模型穩(wěn)定泛函與最大光滑穩(wěn)定泛函加入到正則化項(xiàng)中，可有效降低反演問題的不適定性，但地質(zhì)體邊界刻畫模糊，尤其是對(duì)射線分布稀疏的地質(zhì)體分辨率較低。

近年來，模糊C 均值（Fuzzy C-means，F(xiàn)CM） 聚類算法已用于地球物理反演［17-20］。該算法通過反復(fù)修改聚類中心及每個(gè)網(wǎng)格單元對(duì)聚類中心的隸屬度，實(shí)現(xiàn)對(duì)反演結(jié)果的自動(dòng)分類［17］。此外，通過將先驗(yàn)信息作為參考聚類中心的方式，將FCM 聚類算法從無監(jiān)督學(xué)習(xí)變成監(jiān)督學(xué)習(xí)，引導(dǎo)每個(gè)網(wǎng)格單元向先驗(yàn)信息的方向修正，可有效提高反演結(jié)果的精度［19-20］。

本文在地震初至波旅行時(shí)成像中，采用非結(jié)構(gòu)三角形進(jìn)行網(wǎng)格剖分，采用Sethian 等［8］提出的FMM 算法進(jìn)行正演計(jì)算，在地震旅行時(shí)正演過程中直接計(jì)算靈敏度但不存儲(chǔ)射線路徑［16］。采用以模型靈敏度信息為依據(jù)的漸進(jìn)反演網(wǎng)格優(yōu)化技術(shù)［21-22］，在最小結(jié)構(gòu)模型約束正則化反演基礎(chǔ)上，引入FCM 聚類模型約束，實(shí)現(xiàn)了基于FCM 聚類模型約束的二維地震初至波旅行時(shí)漸進(jìn)反演；對(duì)比了無監(jiān)督學(xué)習(xí)與先驗(yàn)信息監(jiān)督學(xué)習(xí)反演結(jié)果，并應(yīng)用簡(jiǎn)單模型討論了FCM 聚類模型約束的權(quán)重、先驗(yàn)信息引導(dǎo)項(xiàng)權(quán)重等參數(shù)選取策略；最后，通過實(shí)測(cè)數(shù)據(jù)的反演驗(yàn)證了方法的有效性和實(shí)用性。

1 方法原理

1.1 地震初至旅行時(shí)反演

在最小結(jié)構(gòu)模型約束正則化技術(shù)［23］的基礎(chǔ)上，引入FCM 聚類模型約束項(xiàng)，構(gòu)建地震初至波旅行時(shí)反演的目標(biāo)函數(shù)

式中：m為模型速度向量；φd（m）為初至波旅行時(shí)觀測(cè)數(shù)據(jù)dobs與正演響應(yīng)A（m）間的數(shù)據(jù)擬合項(xiàng)；φm（m）為模型擬合項(xiàng)；φFCM（m）為FCM 模型約束項(xiàng)；λ為正則化因子；β為FCM 聚類模型約束項(xiàng)的權(quán)重；Wm為模型誤差矩陣；mref為參考模型向量；Wd為數(shù)據(jù)方差矩陣，形式為

式中：σi為的數(shù)據(jù)方差；N為觀測(cè)的初至波旅行時(shí)數(shù)據(jù)個(gè)數(shù)；ζ為極小的正有理數(shù)，保證分母不為0。

本文采用最小結(jié)構(gòu)穩(wěn)定因子［24］作為模型擬合項(xiàng)

式中：αs為最小模型穩(wěn)定泛函的比例系數(shù)；αx、αz分別為最大光滑穩(wěn)定泛函在x和z方向的比例系數(shù)；?m/?x和?m/?z分別表示模型向量在x和z方向的偏導(dǎo)數(shù)；積分是對(duì)任意單元的面積分。

1.2 FCM 聚類模型約束

FCM 聚類算法的思想是在迭代過程中，反復(fù)修改聚類中心和聚類中心的隸屬度，是一種無監(jiān)督學(xué)習(xí)。FCM 聚類的優(yōu)化目標(biāo)函數(shù)為

式中：C為聚類中心個(gè)數(shù)；M為模型單元總數(shù)；mj為第j個(gè)模型單元速度；隸屬度uj，l表示第j個(gè)模型單元屬于第l個(gè)參考簇的程度；模糊化參數(shù)q控制隸屬度uj，l模糊化程度，本文取q=2；vl為第l個(gè)簇中心的速度值。uj，l為非負(fù)數(shù)，滿足

為了改善反演效果，將先驗(yàn)信息引入FCM 聚類目標(biāo)函數(shù)φFCM（m）［17］，將聚類算法從無監(jiān)督學(xué)習(xí)升級(jí)為監(jiān)督學(xué)習(xí)，即

式中：tl為第l個(gè)參考聚類中心；κl為第l個(gè)參考聚類中心的權(quán)重因子。

φFCM（m）的第一項(xiàng)可表示為

式中

1.3 目標(biāo)函數(shù)優(yōu)化求解

本文反演采用高斯—牛頓法求解模型參數(shù)向量m。假設(shè)第n+1 次迭代，對(duì)正演算子A（m）做一階泰勒展開，并將第n+1 次迭代目標(biāo)函數(shù)φ（n+1）（m）對(duì)Δm（n）進(jìn)行求導(dǎo)并取0，可得高斯—牛頓方程

式（9）可寫為

采用穩(wěn)定雙共軛梯度法（BICGSTAB）［21］求解Δm(n)

再根據(jù)線性搜索步長(zhǎng)更新模型參數(shù)

式中η為搜索步長(zhǎng)，可表示為

式中ξ為極小的正實(shí)數(shù)，確保分母不為0。

1.4 基于模型靈敏度的自適應(yīng)反演

本文采用變網(wǎng)格技術(shù)實(shí)現(xiàn)地震初至波旅行時(shí)反演。在反演的初期采用較粗的網(wǎng)格，通過減少單元數(shù)降低反演的多解性，再利用模型靈敏度信息進(jìn)行自適應(yīng)網(wǎng)格加密［22］。以上一輪反演結(jié)果作為網(wǎng)格細(xì)化后反演的初始和參考模型，確保反演的穩(wěn)定性，同時(shí)改善反演效果。定義模型靈敏度為

式中Jij為靈敏度矩陣的元素。模型靈敏度Sj表示單元的模型參數(shù)mj改變時(shí)對(duì)全部觀測(cè)數(shù)據(jù)集的影響。對(duì)于測(cè)點(diǎn)處單元的模型靈敏度較大，網(wǎng)格的細(xì)化會(huì)從測(cè)點(diǎn)逐漸向反演區(qū)域的深部進(jìn)行，因此具有對(duì)反演區(qū)域整體優(yōu)化的特點(diǎn)?；谀Ｐ挽`敏度的自適應(yīng)反演過程如下。

（1）讀入觀測(cè)數(shù)據(jù)，設(shè)置每重網(wǎng)格最大的迭代次數(shù)、最小單元面積、網(wǎng)格優(yōu)化比例（為反演單元總數(shù)的2%）。

（2）從數(shù)據(jù)信息生成正演模型的幾何描述，設(shè)置反演范圍，并做非結(jié)構(gòu)三角形網(wǎng)格剖分。

（3）對(duì)當(dāng)前網(wǎng)格下使用BICGSTAB 最優(yōu)化算法求解高斯—牛頓優(yōu)化的反演目標(biāo)函數(shù)，直至滿足BICGSTAB 的迭代次數(shù)，并求取模型改變量。

（4）通過模型靈敏度網(wǎng)格優(yōu)化方案計(jì)算單元模型靈敏度信息，按照網(wǎng)格優(yōu)化比例選取需要細(xì)化的網(wǎng)格單元，進(jìn)行網(wǎng)格的自適應(yīng)加密。

（5）將上一輪反演結(jié)果作為細(xì)化后反演的初始和參考模型，并判斷是否達(dá)到反演終止的條件。滿足則終止反演，否則返回第三步開始網(wǎng)格細(xì)化后的反演。

1.5 參數(shù)選取

地震初至波旅行時(shí)反演的目標(biāo)函數(shù)φ（m）中的正則化因子λ是平衡數(shù)據(jù)擬合項(xiàng)與模型擬合項(xiàng)的關(guān)鍵參數(shù)。本文采用李曼等［25］提出的改進(jìn)“L 曲線”算法求取λ。該算法無需利用導(dǎo)數(shù)信息，可以自動(dòng)舍棄錯(cuò)誤的曲率點(diǎn)，能提高正則化因子的求取精度。

引入FCM 聚類模型約束的同時(shí)，也引入了聚類中心的個(gè)數(shù)C、參考聚類中心tl、參考聚類中心的權(quán)重因子κl、FCM 聚類模型約束項(xiàng)權(quán)重β等眾多參數(shù)。這些參數(shù)勢(shì)必會(huì)影響反演的穩(wěn)定性及效果。因此如何選取合適參數(shù)是引入FCM 聚類模型約束反演的一個(gè)關(guān)鍵問題。其中，聚類中心的個(gè)數(shù)C及參考聚類中心tl可通過獲取有效且精確的先驗(yàn)信息確定。參考聚類中心的權(quán)重因子κl的作用是引導(dǎo)聚類中心向真值靠近，F(xiàn)CM 聚類模型約束項(xiàng)權(quán)重β的作用是控制反演中FCM 聚類模型約束項(xiàng)的比重。為此，對(duì)κl和β參數(shù)不同組合的反演效果進(jìn)行了對(duì)比分析，討論了兩個(gè)參數(shù)在反演過程中的作用；結(jié)合漸進(jìn)反演網(wǎng)格優(yōu)化技術(shù)的特點(diǎn)，提出了一種FCM 聚類模型約束參數(shù)的自動(dòng)選取方案。該方案在選取一個(gè)適當(dāng)?shù)木垲愔行臋?quán)重因子κ=100 后，通過減少反演前期FCM 聚類模型約束比重，讓正則化項(xiàng)占主導(dǎo)，后期增大FCM 聚類模型約束比重實(shí)現(xiàn)FCM 聚類模型約束。將每重網(wǎng)格的FCM 聚類模型約束項(xiàng)權(quán)重修改為

式中K為每重網(wǎng)格的迭代次數(shù)。

2 算例分析

2.1 模型算例1

為了驗(yàn)證算法的有效性，設(shè)計(jì)了圖1 左所示的理論模型。在背景速度為2000 m/s的均勻半空間中，設(shè)計(jì)了兩個(gè)尺寸為30 m×30 m 的高速異常體，速度均為3000 m/s，頂面埋深分別為10和90 m。

圖1 異常體模型（左）及最小結(jié)構(gòu)模型約束反演結(jié)果（右）

采用雙井激發(fā)、井中和地面同時(shí)接收的觀測(cè)系統(tǒng)。井中深度2.5 m 處開始設(shè)置炮點(diǎn)和檢波器各60個(gè)，地表均勻布設(shè)60 個(gè)檢波器，炮點(diǎn)距和檢波點(diǎn)距都為2.5 m，井間距為120 m?；谏鲜鰠?shù)進(jìn)行地震初至旅行時(shí)正演，將理論旅行時(shí)加入5%的隨機(jī)噪聲作為觀測(cè)數(shù)據(jù)（共14400個(gè)）。

反演的初始速度模型是速度為2000 m/s 的均勻半空間，初始正則化因子λ設(shè)置為10000。最小結(jié)構(gòu)模型約束反演結(jié)果如圖1 右所示，可以大致判斷兩個(gè)高速異常體的位置，淺部的異常體的規(guī)模和速度與實(shí)際值接近，但深部異常體規(guī)模較大，速度與實(shí)際值相差亦較大。

為了分析基于模型靈敏度漸進(jìn)反演過程，對(duì)于圖1右所示反演模型，設(shè)計(jì)了四重網(wǎng)格，第一重網(wǎng)格的迭代次數(shù)為5，其他的最大迭代次數(shù)為10，反演過程的網(wǎng)格變化如圖2 所示，從測(cè)點(diǎn)處逐漸向反演區(qū)域內(nèi)部細(xì)化。表1 給出了每重網(wǎng)格的反演單元個(gè)數(shù)、迭代次數(shù)及反演開始與結(jié)束的均方根誤差（RMS）。單元數(shù)由1961 逐漸增加到3883；初始網(wǎng)格和第一次細(xì)化網(wǎng)格反演的RMS 下降很快，后兩次的RMS 下降量減少并趨于穩(wěn)定。

表1 漸進(jìn)反演次數(shù)、RMS 統(tǒng)計(jì)

圖2 多重網(wǎng)格細(xì)化示意圖

首先根據(jù)精確的先驗(yàn)信息，設(shè)置了FCM 聚類模型約束反演的兩個(gè)參考聚類中心：v=2000 和3000 m/s。其次，為分析參考聚類中心的權(quán)重因子κ（本文算例中不同參考聚類中心的權(quán)重因子設(shè)為相同）和FCM 聚類模型約束項(xiàng)的權(quán)重因子β對(duì)反演結(jié)果的影響，設(shè)置κ和β分別為1、10、100、1000、10000，對(duì)比不同參數(shù)組合的反演結(jié)果（圖3）。

由圖3可見：

（1）當(dāng)聚類中心權(quán)重因子κ=1或10時(shí)，隨β增大，反演結(jié)果的聚類特征更明顯，但速度值較真值差異變大（圖3的第1、第2行）。

（2）當(dāng)FCM 聚類模型約束項(xiàng)權(quán)重β=1、10、100時(shí)，隨聚類中心權(quán)重因子κ的增大，異常體的特征改善不明顯，也沒有明顯的聚類特征（圖3 的第1～第3 列）。原因是聚類中心權(quán)重因子κ較大時(shí)引導(dǎo)每個(gè)單元的速度向真實(shí)值靠近，但控制聚類程度的FCM 聚類模型約束項(xiàng)權(quán)重β較小，聚類效果不明顯。

（3）適當(dāng)?shù)膮?shù)組合可以得到理想的反演結(jié)果，如κ=100、β=1000 或10000 時(shí)的反演結(jié)果異常體尺寸與真實(shí)模型接近（圖3 第3 行的第4 和第5 列），而β=10000 時(shí)反演的異常體速度與真實(shí)更接近（圖3第3行的第5列）。

（4）越大的κ和β組合的反演結(jié)果的聚類特征越明顯（圖3第5列的第4和第5行），但異常體的尺寸偏小。

由以上分析可知，適當(dāng)?shù)木垲愔行臋?quán)重因子κ能有效地引導(dǎo)聚類中心向真值靠近，而FCM 聚類模型約束項(xiàng)權(quán)重β可以控制聚類程度。因此，本文采用前文提出的FCM 聚類模型約束參數(shù)的選取方案，即通過減小反演前期FCM 聚類模型約束項(xiàng)比重β，讓模型正則化項(xiàng)占主導(dǎo)，后期增大FCM 聚類模型約束項(xiàng)的權(quán)重實(shí)現(xiàn)FCM 模型約束。新參數(shù)選取方案的反演結(jié)果如圖4所示，與未加入FCM 聚類模型約束的反演結(jié)果（圖1 右）相比，背景更干凈，兩個(gè)高速異常體尺寸、位置及速度值與實(shí)際更吻合，邊界刻畫清晰。與圖3第3行第5列的反演結(jié)果接近。

圖4 新參數(shù)選取方案的FCM 聚類模型約束反演結(jié)果

為進(jìn)一步分析FCM 聚類模型約束的地震初至波旅行時(shí)反演效果，對(duì)最小結(jié)構(gòu)模型約束和加入FCM模型約束反演結(jié)果的速度頻數(shù)進(jìn)行統(tǒng)計(jì)（圖5）?？梢?，F(xiàn)CM 聚類模型約束通過引入先驗(yàn)信息，反演結(jié)果中出現(xiàn)兩個(gè)峰值與實(shí)際給定的兩個(gè)聚類中心接近，形成了有效的聚類效果。

（4）該水源地遠(yuǎn)期應(yīng)對(duì)措施應(yīng)嚴(yán)格執(zhí)行水功能區(qū)劃的目標(biāo)要求，通過流域協(xié)同治理，削減流域污染負(fù)荷；同時(shí)構(gòu)建區(qū)域協(xié)同管理機(jī)制，實(shí)現(xiàn)信息共享。

圖5 FCM 聚類模型約束引入前（左）、后（右）反演結(jié)果的速度頻數(shù)統(tǒng)計(jì)直方圖對(duì)比

2.2 模型算例2

為進(jìn)一步驗(yàn)證基于FCM 聚類模型約束的地震初至波旅行時(shí)反演能提高成像分辨率，設(shè)計(jì)了如圖6 左所示的層狀模型，區(qū)域1、區(qū)域2 和區(qū)域3 是速度分別為1000、2000、3000 m/s 的均勻介質(zhì)。采用與模型算例1相同的觀測(cè)方法和參數(shù)進(jìn)行正演計(jì)算。

圖6 層狀模型（左）與加入FCM 聚類模型約束前（中）、后（右）反演結(jié)果對(duì)比

反演的初始模型是速度為2000 m/s 的均勻半空間，初始正則化因子λ為10000，設(shè)置了v=1000、2000、3000 m/s 三個(gè)參考聚類中心；采用模型算例1中FCM 聚類模型約束參數(shù)的選取方案。

圖6 中為最小結(jié)構(gòu)模型約束的反演結(jié)果，可以大致判斷層狀速度帶，但各層邊界模糊，界面附近的速度值與實(shí)際相差較大；圖6右為加入FCM 聚類模型約束的反演結(jié)果，與圖6 中相比，能明顯區(qū)分3 個(gè)區(qū)域的層狀速度帶，界面位置較準(zhǔn)確，刻畫清晰，速度值與真值更吻合。

圖6 兩種方法反演結(jié)果的速度頻數(shù)統(tǒng)計(jì)如圖7所示，加入FCM 聚類模型約束后，三個(gè)峰值的速度與引入的三個(gè)先驗(yàn)信息吻合，形成了有效的聚類特征。

圖7 層狀模型加入FCM 聚類模型約束前（左）、后（右）反演結(jié)果的速度頻數(shù)統(tǒng)計(jì)直方圖對(duì)比

2.3 實(shí)測(cè)數(shù)據(jù)反演

使用強(qiáng)風(fēng)化層中的孤石（微風(fēng)化或未風(fēng)化的花崗巖）探測(cè)的井中地震初至波旅行時(shí)數(shù)據(jù)驗(yàn)證本文方法的有效性和實(shí)用性。實(shí)測(cè)數(shù)據(jù)由單井激發(fā)、單井接收觀測(cè)系統(tǒng)采集。在孤石兩側(cè)各布置一個(gè)深度為20 m的鉆孔，鉆孔間距為9.4 m，共布置17 個(gè)炮點(diǎn)和36 個(gè)檢波器，左側(cè)從3.0 m 深度開始布設(shè)炮點(diǎn)，炮點(diǎn)距為1 m；右側(cè)從深度2.5 m 開始布置檢波器，道間距為0.5 m。共采集有效旅行時(shí)觀測(cè)數(shù)據(jù)612個(gè)。

在測(cè)線5 m 處的驗(yàn)證鉆孔揭示：0～2.6 m 為人工堆積黏性土；2.6～7.0 m 為坡積而成的粉質(zhì)黏土；7.0～9.8 m 為花崗巖風(fēng)化殘積而成的砂質(zhì)黏土；9.8～12.2 m為花崗巖微風(fēng)化體（孤石）；12.2～16.5 m為花崗巖風(fēng)化殘積而成的砂質(zhì)黏性土；16.5～33.6 m為全風(fēng)化層，局部夾強(qiáng)風(fēng)化花崗巖；穩(wěn)定的潛水面深度為2.5 m。

圖8 左為未加入FCM 聚類模型約束的反演速度剖面，可以大致判斷孤石（高速異常體）位置，與鉆孔信息吻合，但尺寸較小，且速度與背景差異不夠明顯。圖8右為加入FCM 聚類模型約束后的反演速度剖面，為了凸顯高速異常體的明顯特征，設(shè)置了v=2000 和3000 m/s兩個(gè)參考聚類中心，采用模型算例1中FCM聚類模型約束參數(shù)的選取方案。由圖8右能清晰地判斷孤石的位置及邊界，異常體的尺寸和速度值都與實(shí)際鉆孔信息吻合。

3 結(jié)束語(yǔ)

為了提高地震初至波旅行時(shí)對(duì)淺層結(jié)構(gòu)的成像精度，本文提出了基于FCM 聚類模型約束的地震初至波旅行時(shí)反演方法。通過理論研究、模型數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)試算，取得以下幾點(diǎn)認(rèn)識(shí)：

（1）FCM 聚類模型約束易與傳統(tǒng)正則化反演結(jié)合，可有效提高反演結(jié)果的分辨率；

（2）FCM 聚類模型約束通過引入精確的速度先驗(yàn)信息，可有效提高速度反演精度；

（3）FCM 聚類模型約束權(quán)重、先驗(yàn)信息引導(dǎo)項(xiàng)權(quán)重直接影響反演結(jié)果，本文的自動(dòng)選擇策略有效改善了反演效果，降低了參數(shù)選擇難度。