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

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

（1. 東華理工大學地球物理與測控技術學院，江西南昌 330013；2. 湖南省交通規劃勘察設計研究院有限公司，湖南長沙 410008；3. 廣西科技大學，廣西柳州 545006）

0 引言

地震旅行時成像是以地震波在不同介質中傳播規律為基礎，在地表或井中觀測地震信號，通過旅行時反演重構地質體內部速度分布的一種地球物理勘探方法。該方法廣泛應用于工程檢測、油氣勘察、地球內部結構研究等［1-3］。地震初至波旅行時成像能識別、追蹤目標地質體，穩定性較好，對淺層成像精度較高。

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

近年來，模糊C 均值（Fuzzy C-means，FCM） 聚類算法已用于地球物理反演［17-20］。該算法通過反復修改聚類中心及每個網格單元對聚類中心的隸屬度，實現對反演結果的自動分類［17］。此外，通過將先驗信息作為參考聚類中心的方式，將FCM 聚類算法從無監督學習變成監督學習，引導每個網格單元向先驗信息的方向修正，可有效提高反演結果的精度［19-20］。

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

1 方法原理

1.1 地震初至旅行時反演

在最小結構模型約束正則化技術［23］的基礎上，引入FCM 聚類模型約束項，構建地震初至波旅行時反演的目標函數

式中：m為模型速度向量；φd（m）為初至波旅行時觀測數據dobs與正演響應A（m）間的數據擬合項；φm（m）為模型擬合項；φFCM（m）為FCM 模型約束項；λ為正則化因子；β為FCM 聚類模型約束項的權重；Wm為模型誤差矩陣；mref為參考模型向量；Wd為數據方差矩陣，形式為

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

本文采用最小結構穩定因子［24］作為模型擬合項

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

1.2 FCM 聚類模型約束

FCM 聚類算法的思想是在迭代過程中，反復修改聚類中心和聚類中心的隸屬度，是一種無監督學習。FCM 聚類的優化目標函數為

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

為了改善反演效果，將先驗信息引入FCM 聚類目標函數φFCM（m）［17］，將聚類算法從無監督學習升級為監督學習，即

式中：tl為第l個參考聚類中心；κl為第l個參考聚類中心的權重因子。

φFCM（m）的第一項可表示為

式中

1.3 目標函數優化求解

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

式（9）可寫為

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

再根據線性搜索步長更新模型參數

式中η為搜索步長，可表示為

式中ξ為極小的正實數，確保分母不為0。

1.4 基于模型靈敏度的自適應反演

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

式中Jij為靈敏度矩陣的元素。模型靈敏度Sj表示單元的模型參數mj改變時對全部觀測數據集的影響。對于測點處單元的模型靈敏度較大，網格的細化會從測點逐漸向反演區域的深部進行，因此具有對反演區域整體優化的特點?；谀Ｐ挽`敏度的自適應反演過程如下。

（1）讀入觀測數據，設置每重網格最大的迭代次數、最小單元面積、網格優化比例（為反演單元總數的2%）。

（2）從數據信息生成正演模型的幾何描述，設置反演范圍，并做非結構三角形網格剖分。

（3）對當前網格下使用BICGSTAB 最優化算法求解高斯—牛頓優化的反演目標函數，直至滿足BICGSTAB 的迭代次數，并求取模型改變量。

（4）通過模型靈敏度網格優化方案計算單元模型靈敏度信息，按照網格優化比例選取需要細化的網格單元，進行網格的自適應加密。

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

1.5 參數選取

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

引入FCM 聚類模型約束的同時，也引入了聚類中心的個數C、參考聚類中心tl、參考聚類中心的權重因子κl、FCM 聚類模型約束項權重β等眾多參數。這些參數勢必會影響反演的穩定性及效果。因此如何選取合適參數是引入FCM 聚類模型約束反演的一個關鍵問題。其中，聚類中心的個數C及參考聚類中心tl可通過獲取有效且精確的先驗信息確定。參考聚類中心的權重因子κl的作用是引導聚類中心向真值靠近，FCM 聚類模型約束項權重β的作用是控制反演中FCM 聚類模型約束項的比重。為此，對κl和β參數不同組合的反演效果進行了對比分析，討論了兩個參數在反演過程中的作用；結合漸進反演網格優化技術的特點，提出了一種FCM 聚類模型約束參數的自動選取方案。該方案在選取一個適當的聚類中心權重因子κ=100 后，通過減少反演前期FCM 聚類模型約束比重，讓正則化項占主導，后期增大FCM 聚類模型約束比重實現FCM 聚類模型約束。將每重網格的FCM 聚類模型約束項權重修改為

式中K為每重網格的迭代次數。

2 算例分析

2.1 模型算例1

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

圖1 異常體模型（左）及最小結構模型約束反演結果（右）

采用雙井激發、井中和地面同時接收的觀測系統。井中深度2.5 m 處開始設置炮點和檢波器各60個，地表均勻布設60 個檢波器，炮點距和檢波點距都為2.5 m，井間距為120 m?；谏鲜鰠颠M行地震初至旅行時正演，將理論旅行時加入5%的隨機噪聲作為觀測數據（共14400個）。

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

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

表1 漸進反演次數、RMS 統計

圖2 多重網格細化示意圖

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

由圖3可見：

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

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

（3）適當的參數組合可以得到理想的反演結果，如κ=100、β=1000 或10000 時的反演結果異常體尺寸與真實模型接近（圖3 第3 行的第4 和第5 列），而β=10000 時反演的異常體速度與真實更接近（圖3第3行的第5列）。

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

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

圖4 新參數選取方案的FCM 聚類模型約束反演結果

為進一步分析FCM 聚類模型約束的地震初至波旅行時反演效果，對最小結構模型約束和加入FCM模型約束反演結果的速度頻數進行統計（圖5）。可見，FCM 聚類模型約束通過引入先驗信息，反演結果中出現兩個峰值與實際給定的兩個聚類中心接近，形成了有效的聚類效果。

（4）該水源地遠期應對措施應嚴格執行水功能區劃的目標要求，通過流域協同治理，削減流域污染負荷；同時構建區域協同管理機制，實現信息共享。

圖5 FCM 聚類模型約束引入前（左）、后（右）反演結果的速度頻數統計直方圖對比

2.2 模型算例2

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

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

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

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

圖6 兩種方法反演結果的速度頻數統計如圖7所示，加入FCM 聚類模型約束后，三個峰值的速度與引入的三個先驗信息吻合，形成了有效的聚類特征。

圖7 層狀模型加入FCM 聚類模型約束前（左）、后（右）反演結果的速度頻數統計直方圖對比

2.3 實測數據反演

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

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

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

3 結束語

為了提高地震初至波旅行時對淺層結構的成像精度，本文提出了基于FCM 聚類模型約束的地震初至波旅行時反演方法。通過理論研究、模型數據和實測數據試算，取得以下幾點認識：

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

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

（3）FCM 聚類模型約束權重、先驗信息引導項權重直接影響反演結果，本文的自動選擇策略有效改善了反演效果，降低了參數選擇難度。