基于動態空腔膨脹模型的土壤Mohr-Coulomb參數反演方法

鄒 歐，薛 健，李 娜，王 鑫，楊 榮，李玉瓊*

（1.中國科學院大學 工程科學學院，北京 100049；2.中國科學院 力學研究所 流固耦合系統力學重點實驗室，北京 100190；3.中國科學院 力學研究所 非線性力學國家重點實驗室，北京 100190）

1 引言

土壤強度參數解譯在定量評估軍事裝備機動性和指導裝備輕量化結構設計等方面具有重要作用，在地質勘測、裝備全域機動、地外行星科考等重大需求領域中有廣泛的應用背景［1-5］。傳統“取樣-室內實驗”是典型的非原位測試方法，取樣造成的土體擾動——特別是原位應力的解除——使得所采集試樣的力學性質不能真實反映土壤原本的力學特征，從而影響對土壤承載性能及裝備車輛通行性的可靠評估［6］。因此，基于現場原位貫入觸探測試并結合相關理論實現土壤原位力學參數的精準、快速反演的需求日益緊迫。目前，土壤原位貫入測試方法主要有靜力觸探（Cone Penetration Test，CPT）和自由落體貫入觸探。靜力觸探能獲得土壤強度參數隨深度的變化情況，一般需要外載荷持續作用以保持勻速貫入（2.5 cm/s），設備安裝耗時長，操作費用高［7-8］，難以滿足復雜地質環境和大范圍區域內土壤力學參數的測量，限制了其在裝備機動性評估方面的應用。自由落體貫入觸探技術是在靜力觸探基礎上發展起來的一種新型原位測試技術，它依靠自由落體貫入儀（Free Fall Penetrometer，FFP）以一定初始速度貫入土壤，通過傳感器實時記錄貫入過程中貫入儀的加速度和受到的阻力，來反演土壤強度參數。FFP 結構輕便、操作簡單且較為經濟，能夠實現復雜環境和較大區域內土壤強度參數的快速測量［7，9］。

目前，在海洋工程中，已有相關研究利用FFP 獲取土壤強度參數。Aubeny 等人［10］提供了基于數值模擬建立的不排水抗剪強度計算方法，給出了承載力系數的經驗公式；Stoll 等人［11］通過對貫入實驗所得數據進行最小二乘擬合，得到了貫入阻力和最大減速度之間的關系。研究表明貫入過程中的阻力可表示為一個與貫入速度有關的阻力函數［12-14］。通過確定土壤強度參數與阻力函數之間的關系，結合貫入儀上傳感器記錄的加速度數據可反演土壤強度參數。True［15］給出了一個描述FFP 貫入海床土過程的半經驗阻力函數，其中考慮了端部阻力、側壁摩擦阻力、土壤拖曳阻力以及土體浮力的作用。Chow 等人［9］利用離心機實驗中傳感器測量的加速度數據，反演了飽和軟黏土的不排水抗剪強度，并考慮了率效應對軟黏土強度的影響。劉君等人［16］采用數值仿真的方法模擬了貫入儀在軟黏土中的運動，研究其所受阻力與貫入速度、土壤強度和界面摩擦系數之間的關系，提出了一種反演軟黏土不排水抗剪強度的方法。以上反演方法均基于True 提出的半經驗阻力函數，其中包括承載力系數，界面摩擦系數，拖曳阻力系數等許多經驗參數，當缺少場地土壤條件前期資料時，經驗參數取值會對反演方法的準確性產生重要影響。此外，以上方法僅能針對海洋工程中常見飽和軟黏土的不排水抗剪強度進行反演。這表明，由于實際地質環境的復雜性，基于自由落體貫入技術的沖擊貫入機理還有待更深入的研究，參數反演方法還需完善，其適用性還需拓展。

在貫入阻力的理論研究方面，現有方法主要包括承載力理論［17］，應變路徑理論［18］和空腔膨脹理論［19-26］。承載力理論僅考慮土體破壞區域的整體平衡，而破壞區域內各點并不處于平衡狀態。Janbu 等人［27］通過黏性土的靜力觸探實驗表明，利用承載力理論求得的貫入阻力相比于實測值低約40%。應變路徑法只能近似求解飽和黏性土的情況［18］，對于砂土，由于貫入時的流場難以構造，該方法并不適用。

空腔膨脹理論是研究侵徹問題的重要模型。由于FFP 的剛度遠大于土壤，在貫入過程中變形較小，因此FFP 貫入土壤的過程可視作剛性彈侵徹靶體的過程。若將貫入儀看作由一系列與其內切的球形空腔組成，則貫入儀與土體之間的相互作用就可看作這些球形空腔在無限大介質中均勻膨脹時與土體介質間的相互擠壓［28］。通過求解介質中空腔膨脹的徑向應力分布，就可得到彈體貫入過程中受到的貫入阻力。Bishop［19］首先給出了準靜態柱/球形空腔膨脹理論的控制方程，并利用它求解了尖錐體對金屬靶的壓入阻力問題。Chadwick［20］將這個理論應用到理想土壤材料，結合Mohr-Coulomb 屈服準則，得到了土壤徑向位移的表達式。Cao 等人［21］假定土的本構為修正劍橋模型，通過求解準靜態空腔膨脹理論控制方程，得到了土體內的超孔隙水壓力和總應力。Hill［22］則首先給出動態空腔膨脹理論控制方程，得到了不可壓縮理想彈塑性靶體對彈體的阻力。基于動態空腔膨脹理論，Forrestal 等［23-25］開展了針對于不同靶體材料（如巖石、混凝土及土壤等）的侵徹研究，大大拓寬了這一理論的應用范圍。Shi 等人［26］針對干砂的沖擊貫入，建立了考慮壓縮性的動態空腔膨脹理論。動態空腔膨脹模型準確地刻畫了穩定貫入過程中貫入阻力和貫入速度平方間的線性關系，展現了清晰的物理圖像，而土壤強度參數則反映在阻力系數當中，這為土壤強度參數反演提供了理論基礎。

本文基于沖擊貫入觸探技術和動態空腔膨脹模型，開展了土壤Mohr-Coulomb 參數（內聚力和內摩擦角）反演方法的研究。利用土壤沖擊貫入實驗中測得的阻力-速度曲線、加速度時程曲線、速度時程曲線以及貫入深度-初速度關系，對貫入阻力系數進行了擬合，進而對Mohr-Coulomb 參數進行了反演。最后，對參數可解域和參數敏感性進行了討論。

2 動態球形空腔膨脹模型

一般來講，在對土壤的沖擊貫入過程中，貫入儀可視作剛體。此時，動態球形空腔膨脹模型為計算彈-靶相互作用提供了一個簡單的物理圖像。這個模型將貫入儀看作是由一系列與其內切的球形空腔組成，空腔球心在貫入儀的對稱軸上［28］。貫入儀與土體的相互作用可看作是這些球形空腔在無限大介質中均勻膨脹時與土體介質間的相互擠壓。這樣一來，貫入儀側表面上每一點的法向應力可近似為球形空腔膨脹時內表面的徑向應力，它與空腔膨脹的速度有關，而空腔膨脹速度為貫入儀側表面上該點的法向速度分量［25，28］。由于球對稱性的引入，經過以上簡化后，貫入儀側表面上的壓力分布就可轉化為一個一維問題進行求解。貫入儀圖示及球形空腔膨脹模型如圖1。

球形空腔的膨脹將導致土體介質出現相應的動態響應，形成彈塑性分區［25］（如圖1（b））。記空腔膨脹速度為V、彈塑性邊界擴張速度為c、初始時刻空腔的半徑為零，則在t時刻塑性區的范圍為（Vt，ct），彈性區的范圍為（ct，+∞）。在彈性區內，土壤的本構關系采用Hooke 定律來描述，并且假定彈性區內土壤不可壓縮。在塑性區內，采用體積鎖定的假設［25］來描述土壤的應力應變關系，即發生塑性變形后土壤的體積應變為常數η*，不隨平均應力的增大而變化。在這個假定下，靜水壓力將不影響土壤的后繼屈服。記ρ0和ρ*分別是土壤的初始密度和鎖定密度，則：

圖1 貫入儀及球形空腔膨脹模型Fig.1 Penetrometer and spherical cavity expansion model

土壤的屈服準則采用Mohr-Coulomb 準則，在球坐標系下可寫為［25］：

其中：σr，σθ分別為徑向應力和環向應力（以壓為正）；p=（σr+2σθ）/3 為平均應力；τ0，λ是與土壤內聚力C和內摩擦角φ有關的常數，且有［29］：

沖擊貫入過程中，質量、動量及能量守恒定律均需得到滿足。對于貫入阻力和貫入深度的解答，僅利用質量守恒和動量守恒關系即可得到。在球坐標系下，拉格朗日坐標表示的質量守恒方程和動量守恒方程可寫為：

其中：r是拉格朗日坐標，u為徑向位移（以向外為正）。空腔表面處的邊界條件為：

彈塑性界面處的Hugoniot 條件可表示為［24］：

其中，v是介質的徑向速度，下標1、2 分別代表彈塑性界面處的彈性區和塑性區。

利用以上本構關系和邊界條件求解控制方程（4）和（5）可解析地得到貫入儀端部表面的法向應力分布為［25］：

其中：Vz是貫入儀在豎直方向的速度；A，B是與τ0/E，λ和η*有關的參數。根據λ的不同取值，A，B的表達式存在如下三種不同的情況［25］：

式（12）中的第2 式表明，在貫入過程中空腔膨脹速度V與彈塑性邊界擴張速度c的比值始終保持為常數。隨著貫入深度的增加，貫入速度減小，空腔膨脹速度V也逐漸減小，對應的彈塑性區域的擴張逐漸變慢。

對于貫入儀相對土壤垂直貫入的情況，將式（9）在貫入儀端部積分，則得到貫入儀所受到的豎直方向阻力Fz為［25］

其中：

αs和βs是與土壤力學特性有關的系數。其中ψ=s/（2r）是貫入儀頭部的曲徑比，并且有θ0=arcsin［（2ψ-1）/（2ψ）］（見圖1（a））；μ為貫入儀與土壤間的滑動摩擦系數，二者之間的切應力可表示為στ=μσn。如果忽略貫入儀與土體間的摩擦（μ=0），并令：

式（18）和式（19）可進一步簡化為：

由式（9）和式（17）可以看到，在當前模型假設下，貫入儀端部的應力分布及貫入阻力由兩部分貢獻，一項是靜阻力項，另一項是流動阻力項（即速度平方項），而黏性效應和附加質量效應［14］未被考慮。

貫入儀在豎直方向的貫入速度Vz可由牛頓第二定律計算得到，注意到式（17），有：

其中，m為貫入儀的質量。積分上式即得到：

其中，V0=Vz（t=0）是貫入儀的初始貫入速度。貫入過程中貫入儀加速度大小為：

最終貫入深度P為：

3 土壤Mohr-Coulomb 參數的反演

3.1 Mohr-Coulomb 參數C，φ的反演方法

土壤Mohr-Coulomb 參數即土壤的內聚力C和內摩擦角φ。由于動態空腔膨脹模型刻畫了土壤力學特性參數對彈靶相互作用的影響，因此對C和φ的反演可以建立在動態空腔膨脹模型的基礎上。由式（3）可知：

上式表明，要反演內聚力C和內摩擦角φ，實際上就是要求得參數τ0和λ的值。由式（18）、式（19）（或對于不考慮摩擦時的方程（21）和（22））以及A、B的表達式可見，參數τ0、λ包含在阻力系數αs、βs中。一旦 求得αs和βs，聯立方 程（18）和（19）（或對于不考慮摩擦的情況，聯立方程（21）和（22）），即可求得參數τ0和λ的值。

在聯立方程（18）和（19）（或方程（21）和（22））求解參數τ0和λ時，需針對A、B取值的三種情況進行討論。由于參數A、B的表達式過于復雜，一般不能得到τ0和λ的解析表達式。此時，可采用如下幾何化方法進行求解。為簡便起見，忽略掉貫入儀與土體間的摩擦（即令μ=0），不影響對方法的討論。

（1）如果λ≠0 且λ≠3/4：

此時，由式（10）和式（11）可見，參數A、B均為（τ0，λ）的函數。根據式（21）和式（22），αs、βs也均是（τ0，λ）的函數。因此，由式（21）和式（22）可定義如下曲面：

此種情況下若存在τ0和λ的解，那么參數τ0和λ必是曲面Π1和Π2的公共零點。

（2）如果λ=3/4 或λ=0：

此時，由式（13）～式（16）可見，參數A、B僅是τ0的函數。因此αs、βs也僅是τ0的函數。于是，根據式（21）和式（22）可定義如下曲線：

此種情況下若存在τ0的解，則τ0必是曲線L1和L2的公共零點。

至此，對C、φ的反演就歸結為對阻力系數αs和βs的計算，它們可通過貫入阻力、貫入速度、貫入加速度以及最終貫入深度的測量結果計算得到。下面就針對這些情況進行具體說明。

3.2 基于阻力-速度關系進行的反演

在沖擊貫入過程中，最重要的物理量是貫入儀受到的貫入阻力Fz，它綜合反映了靶體介質的物理力學性質以及貫入儀與靶體材料之間的相互作用。動態空腔膨脹模型的核心就是揭示了貫入過程中的穩定階段貫入阻力Fz（t）與速度平方間的線性關系（式（17））。因此，我們首先利用這個關系來反演土壤的Mohr-Coulomb參數。

實際當中，貫入阻力Fz（t）可由加速度傳感器測量得到；而貫入速度Vz（t）則可由加速度時程曲線積分得出。因此，利用測量數據集｛Fz，i，Vz，i｝并基于式（17）進行線性擬合，即可得到參數αs和βs的數值，并進而求得內聚力C和內摩擦角φ。為此，構造誤差函數：

其 中，N是數據采樣點數目。得 到αs，βs的數值后，即可按照上節所述步驟進行C，φ值的反演。下面采用Forrestal 等人［25］報道的實驗結果來對這個反演方法進行詳細說明。

文獻［25］中，通過對不同采樣地點取樣的土壤進行三軸壓縮及單軸壓縮標定實驗，測得的土壤參數如表1 所示。需要說明的是采樣地點和現場貫入實驗地點的土層按照深度可分為三層［25］，第一層（Ⅰ：0～1.8 m）的彈性模量為120 MPa，鎖定體積應變為0.17；第二層（Ⅱ：1.8～3.0 m）的彈性模量為150 MPa，鎖定體積應變為0.10；第三層（Ⅲ：3.0～6.1 m）的彈性模量為210 MPa，鎖定體積應變為0.13。為了簡化模型計算，文獻［25］忽略了土壤的分層，將其視作均勻介質，并取三層土壤的平均彈性模量（160 MPa）和平均鎖定體積應變（0.13）作為土壤基本參數。本文按照同樣的假定進行計算。

表1 Forrestal 實驗土壤參數［25］Tab.1 Soil parameters of Forrestal′s experiments［25］

以文獻［25］中的貫入實驗3 為例，貫入儀質量為23.1 kg，彈體半徑r=47.6 mm，彈頭曲徑比ψ=3。忽略貫入過程中貫入儀與土體間的摩擦。由實驗測得的加速度數據計算得到貫入阻力和貫入速度后，給出的阻力-速度平方關系如圖2 所示。由圖2 可見，除初始階段貫入儀與土體碰撞造成的震蕩以及貫入結束階段的土體擾動外，在中間的穩定貫入階段，貫入阻力與貫入速度的平方之間呈線性關系，這表明動態空腔膨脹模型揭示了貫入過程的關鍵特征。根據式（17），對穩定貫入階段的阻力-速度關系進行線性擬合，結果如圖2 中紅線所示，擬合得到的阻力系數為αs=85.39 kN，βs=2.52 kg/m。圖中A、B兩點分別是擬合范圍的起點和終點，黑色圓點是參與擬合的數據。關于擬合范圍的選取，將在第4.2 節中進行詳細討論。

圖2 貫入阻力Fz和貫入速度平方Vz2間的關系［25］Fig.2 Relationship between penetration resistanceFzand penetration speed squareVz2［25］

由擬合得到的阻力系數αs、βs，可按照式（28）和式（29）畫出曲面Π1和Π2的圖像，如圖3。由圖3（a）、圖3（b）可見，曲面Π1與坐標平面Π1=0 有交線l1（曲面Π1零點的集合）；曲面Π2與坐標平面Π2=0 有交線l2（曲面Π2零點的集合）。將交線l1和l2畫在同一坐標系下，如圖4 所示，可見二者有公共交點（τ0=5.52 MPa、λ=0.58），顯然，這個交點就是曲面Π1和Π2的公共零點。按照第3.1節的分析，此時的τ0和λ就是所要求的解。此外，對于λ=3/4 和λ=0 的情況也應當進行校驗。按照式（30）和式（31）分別畫出曲線L1和L2的圖像（如圖5（a）和圖5（b）所示），可以看到，對于這兩種情況，L1和L2都沒有公共零點，因此不存在τ0和λ的解。綜合以上分析，將τ0=5.52 MPa 和λ=0.58 代入式（27），就得到C、φ的反演結果分別為：

圖3 曲面Π1和Π2的圖像Fig.3 Images of surfacesΠ1andΠ2

圖4 λ≠0 且λ≠3/4 時交線l1和l2的圖像及其交點Fig.4 Intersection linesl1andl2and their intersection point underλ≠0 andλ≠3/4

圖5 λ=3/4 或λ=0 時曲線L1、L2的圖像Fig.5 LinesL1andL2whenλ=3/4 orλ=0

對于實驗中的土壤試樣，Forrestal 等人［25］給出了三軸實驗下剪切強度τ與平均應力p之間的關系，如圖6，此處，剪切強度τ定義為試樣屈服時第一主應力σ1與第三主應力σ3的差［25］（τ=σ1-σ3）。注意到空腔膨脹模型的球對稱性，第一主應力σ1就是空腔膨脹時的徑向應力σr，第三主應力σ3就是空腔膨脹時的環向應力σθ（此時σ2=σ3=σθ），于是Mohr-Coulomb 準則（式（2））可以寫為：

根據上式，對圖6 中的數據進行線性擬合（如圖6 中紅線所示），可解得相應的C、φ值分別為C=3.73 MPa、φ=3.48°。

圖6 剪切強度與平均應力的關系［25］Fig.6 Relationship between shear strength and average stress［25］

對比式（36）所示的反演結果可見，對于內聚力C的反演結果與三軸實驗的結果較為接近，而對內摩擦角φ的反演結果卻與實驗值差別較大。這一定程度上與采用的擬合方法有關。接下來，我們分別基于貫入加速度數據和貫入速度數據對阻力系數進行擬合，并進一步討論擬合結果對方法的依賴。

3.3 基于加速度時程曲線和速度時程曲線進行的反演

貫入加速度與貫入阻力間僅相差m倍，原則上講，由貫入加速度進行阻力系數αs、βs的擬合與采用貫入阻力進行擬合沒有本質差別。但是，在利用阻力-速度關系進行擬合時，貫入速度是由加速度的數值積分獲得，由于進行了數值積分操作，數據處理過程與直接利用加速度時程曲線進行擬合是不同的，因此擬合結果也將產生差別。此外，由式（17）可見，貫入阻力Fz與貫入速度的平方Vz2間呈線性關系，阻力系數αs和βs通過線性最小二乘擬合即可得到；而由式（25）可見，貫入加速度az與時間t呈現非線性關系，αs和βs只能通過非線性最小二乘擬合得到。

基于加速度時程曲線［25］和速度時程曲線對阻力系數的擬合如圖7。根據貫入加速度表達式（25），采 用Levenberg-Marquardt 方法［30］對文 獻［25］中貫入實驗3 給出的加速度數據進行非線性最小二乘擬合，如圖7（a）中紅線所示。擬合的初值設置為根據阻力-速度關系擬合得到的數值，即αs=85.39 kN，βs=2.52 kg/m。由加速度曲線擬合得到的阻力系數分別為αs=89.22 kN、βs=2.69 kg/m。由此得到τ0=5.61 MPa、λ=0.72，Mohr-Coulomb 參數為：

圖7 基于加速度時程曲線［25］和速度時程曲線對阻力系數的擬合Fig.7 Fitting of the resistance coefficients based on acceleration and velocity curves

與三軸實驗測量值相比，C值的相對誤差為29.2%，與上一節中基于阻力-速度關系反演得到的結果接近，φ值的反演誤差仍然較大。此外，圖7（a）給出了按照文獻［25］中取τ0=10 MPa、λ=0 時計算得到的加速度時程曲線（圖7（a）中藍線），將其與紅線相比較可以看到，τ0=10 MPa 和λ=0 并不是測量數據的最小二乘解。

按照相同的非線性最小二乘擬合方法，利用貫入速度表達式（24）亦可對阻力系數αs和βs進行擬合，如圖7（b）所示。此時，同樣將擬合的初值設置為根據阻力-速度關系擬合得到的數值，即αs=85.39 kN，βs=2.52 kg/m。由速度曲線擬合得到的阻力系數為αs=101.48 kN、βs=1.87 kg/m。相應的τ0=7.76 MPa、λ=0.11，Mohr-Coulomb參數為：

將上式與三軸實驗測量結果對比，C值的相對誤差為2.14%，較前兩種方法的反演結果誤差大大減小，φ值的相對誤差也減小為9.77%。

總結起來，對比式（36）、式（38）和式（39）的反演結果可見，利用阻力-速度關系進行擬合的結果與利用加速度直接進行擬合的結果相近，而與利用貫入速度進行擬合的結果差別較大。更重要的是，相較于前兩者的擬合結果，利用貫入速度擬合得到的內聚力C和內摩擦角φ與三軸實驗給出的結果更加接近。以上結論并不是偶然的，它反映了貫入阻力系數對加速度、速度等物理量有不同的影響，這將在第4.2 節中進行更進一步的討論。此外，以上結論為土壤Mohr-Coulomb 參數反演方法的實際應用提供了重要參照。

3.4 基于深度-初速度關系進行的反演

阻力系數αs、βs的數值還可以通過對同一狀態的土壤進行多次貫入實驗后測得的貫入深度P和初始貫入速度V0來擬合。為此，根據式（26）構造誤差函數：

其中，N是貫入實驗次數。

可解得：

將式（44）代入式（43）并令：

實驗中，彈體質量為286 g，彈頭長度為140 mm，曲徑比ψ=8.54［31-32］。土壤的基本參數如表2 所示。

表2 何翔實驗的土壤參數［31］Tab.2 Soil parameters of He Xiang’s experiments［31］

表3 中列出了同一剛性彈以不同初速度多次貫入同一狀態土壤后測得的貫入深度結果。圖8為的圖像及貫入深度-貫入初速度關系。將表3 中的9 個實驗數據代入式（45）可畫出函數f的圖像（圖8（a）），并得到的零點為1.09×10-3s2/m2，聯立式（41）和式（44）可解得αs=276.92 N，βs=0.30 kg/m。取鎖定體積應變η*=0.13，可得到τ0=52.51 kPa，λ=1.16，Mohr-Coulomb 參數為：

表3 何翔貫入實驗結果［31］Tab.3 Results of He Xiang’s penetration experiments［31］

圖8 的圖像及貫入深度-貫入初速度關系Fig.8 Image of the functionand the relationship between penetration depth and initial speed

反演得到的內聚力C的相對誤差為58.74%，內摩擦角φ的相對誤差為3.07%。可見，利用9 個實驗數據整體擬合反演得到的內聚力C與三軸實驗測量結果相差較大，而得到的內摩擦角φ與三軸實驗測量結果一致。

此外，對表3中所列數據進行分組擬合。根據文獻［31］，表3 中所列實驗結果對應三種不同的貫入初速度V0，分別為150 m/s（序號1～3），200 m/s（序號4～6）和250 m/s（序號7～9），在每個速度等級下取1 個數據點，總計3 個數據點作為一組數據（總共27 組），來擬合αs、βs的值，結果見表4。分組擬合及反演得到的內聚力C的平均值及標準差為29.45±20.69 kPa，內摩擦角φ的平均值及標準差為25.00±17.88°。此時，分組擬合并反演得到的C的平均值較整體擬合結果有所提高，但相比61 kPa 仍差別較大；而對于φ的反演，整體擬合結果與實驗測得的30°更加接近。

表4 27 組數據的計算結果Tab.4 Calculation results of the 27 data sets

以上分析表明，利用貫入深度-初速度關系進行反演時，內摩擦角φ的反演結果更加可靠。

4 討論

4.1 參數的可解域

土 壤Mohr-Coulomb 參 數C、φ并 不是任意的，具有一定范圍，相應的阻力系數αs和βs也對應了一定范圍，稱為參數的可解域。只有擬合得到的αs和βs落在可解域范圍內，才能在動態空腔膨脹模型框架下進行C、φ值的反演。

圖9（a）和 圖9（b）分別給出了Forrestal［25］和何翔［31］等人報道的實驗結果對應的參數可解域。圖9（a）中給出了利用阻力-速度關系、加速度時程曲線和速度時程曲線擬合得到的αs、βs值，其均落在參數的可解域內。圖9（b）中的分組擬合值為第3.4 節中利用貫入深度-初速度關系對27 組數據擬合得到的αs、βs值；整體擬合值為利用9 個實驗數據擬合得到的αs、βs值。這些擬合結果也都在參數的可解域內，因此可以利用動態空腔膨脹理論來進行C、φ值的反演。

圖9 動態空腔膨脹模型下參數的可解域Fig.9 Solution domain determined by dynamic cavity expansion model

實際的地質環境往往十分復雜，數據噪聲也難以避免，根據測量數據擬合得到的αs和βs并不總能反演得到Mohr-Coulomb 參數C、φ的數值。此時，就需要對偏離模型的過程進行更細致的分析

參數的可解域是土壤Mohr-Coulomb 參數對阻力系數作出的限制，它反映了土壤物理力學特性對貫入阻力的影響。對于給定土壤，Mohr-Coulomb 參數總是確定的，因此在理想情況下，阻力系數的實驗值必定落在可解域內。但在實際當中，由于噪聲的影響，或簡化的模型未能涵蓋貫入過程的復雜細節，阻力系數的擬合值可能超出參數的可解域。此時，就需要對沖擊貫入的物理過程進行更加詳細的分析，并對相應的數據處理手段進行優化。

4.2 擬合區間的選取及對反演結果的進一步討論

由第3 節的分析可以知道，要準確地反演C、φ的數值，關鍵是要準確地擬合阻力系數αs和βs。第3.1 節至3.3 節利用貫入時程曲線進行擬合時，擬合數據是穩定貫入階段貫入阻力與貫入速度平方呈線性關系的部分（圖2 中A、B之間的數據）。數據擬合范圍對阻力系數的擬合及Mohr-Coulomb 參數的反演有重要影響，在這一節中就來對數據擬合范圍的選取進行討論。

為此選取擬合窗口起點為圖2 中的點A（對應時刻為t=4.84 ms），改變擬合窗口終點tend進行阻力系數αs和βs的擬合以及Mohr-Coulomb 參數的反演，結果如圖10 所示。由圖可見，不論擬合窗口的大小如何改變，基于阻力-速度關系反演的結果與基于加速度曲線反演的結果接近，而與基于速度曲線的反演有較大的差別；并且，對于內摩擦角φ，前兩者的反演結果始終大于后者。

圖10（a）顯示了兩個關鍵擬合范圍。第一個范圍對應三條擬合曲線的交點（tend=t2）。對于內聚力C的反演，當tend＜t2時，基于阻力-速度關系和加速度曲線反演得到的結果小于利用速度曲線反演的結果；當tend＞t2時，基于阻力-速度關系和加速度曲線反演得到的結果大于利用速度曲線反演的結果。在tend=t2的擬合范圍下，三種擬合方法得到的C的反演結果一致。此時反演得到的數值（C=4 MPa）比三軸實驗結果（C=3.73 MPa）大7.24%。另一個關鍵擬合范圍對應于圖10（a）中速度曲線擬合結果的最小值點tend=t1（這一點也是圖10（b）中速度曲線擬合結果的最大值點），此時，由速度曲線擬合并反演得到的內聚力C=3.81 MPa，內摩擦角為φ=3.14°，與三軸試驗結果一致。這一范圍就是在第3.1 節至3.3 節中擬合時所采用的范圍，圖10 中時刻tend=t1就對應圖2 和圖7 中的點B。

圖10 擬合區間大小對C、φ反演結果的影響［25］Fig.10 Influence of the fitting interval on the inversion results ofCandφ［25］

以上分析為反演方法的實際應用提供了一個可行方案：對于沖擊貫入實驗測得的數據，用速度曲線進行參數C（φ）的擬合，改變區間的大小進行搜索，使得C（φ）的擬合值最小（大）時對應的區間即為最終的擬合區間，相應的C、φ值即為最終的反演結果。

由圖2、7、10 可以看到，擬合區間終點tend一般并不能覆蓋整個貫入過程的數據段，圖10 中能夠反演得到Mohr-Coulomb 參數的最大tend為27.80 ms，這是由于參數可解域的限制。當tend大于27.80 ms 后，擬合得到的αs和βs已經超出了參數的可解域，這表明超過27.80 ms 后的貫入過程可能存在模型忽略掉的細節。在利用動態空腔膨脹模型進行參數反演時，假定了土壤介質的初始狀態是均勻的，這一假設保證了模型的球對稱性，而在現場貫入實驗中，土壤實際上是分層的。圖11 給出了文獻［25］中貫入實驗3 的加速度時程曲線并根據貫入儀進入不同土層（Ⅰ、Ⅱ、Ⅲ）的時刻進行了劃分。由第3.2 節的說明可知，這三層土壤的彈性模量和鎖定體積應變是不同的，特別是第三層土壤的彈性模量明顯高于前兩層；此外，第三層土層的厚度也幾乎是前兩層的兩倍，貫入儀大部分的運動時間在第三層土壤中。這表明，第三層土壤對貫入過程的貢獻應當比前兩層土壤的大，然而平均化的假定抹平了這些差別。當擬合區間為圖11 中的AB范圍時，貫入儀在三個土層中的運動時間相差不大，此時，平均化的假定基本適用；而當擬合區間范圍逐漸擴大，比如對于AC段或AD段，此時第三個土層的貢獻逐漸增大（但仍與前兩層可比擬），平均化的假定將產生誤差，此時可能需要考慮土壤分層的影響。總之，準確地反演土壤的力學參數，需要模型能夠準確地把握貫入過程的關鍵細節。

圖11 貫入過程的加速度時程曲線以及貫入儀進入不同土層對應的時間區間［25］Fig.11 Acceleration curve and the time interval corresponding to the penetration instrument entering different soil layers［25］

4.3 參數的敏感性

在求解沖擊貫入問題時，對于貫入阻力和貫入深度的求解屬于正問題。此時，彈性模量和鎖定體積應變等參數對于計算結果的影響不大［25］。然而，對于通過測量貫入阻力求解靶體介質力學特性參數這樣的反問題，反演結果往往具有參數敏感性。這種敏感性不僅與材料參數本身有關還與阻力系數有關。

圖12 和圖13 分別給出了文獻［25］和文獻［31］中，土壤彈性模量E和鎖定體積應變η*對C、φ值反演的影響。由圖可見，隨著彈性模量E的增大，反演得到的C值減小、φ值增大；隨著鎖定體積應變η*的增大，反演得到的C、φ值均增大。

圖12 利用速度曲線［25］進行反演時，參數E、η*對C、φ值反演的影響，其中αs=101.48 kN、βs=1.87 kg/mFig.12 Influence of parametersEandη* on the inversion results ofCandφ，respectively，whereαs=101.48 kN andβs=1.87 kg/m（using the velocity curve［25］）

由圖12（a）和圖13（a）可見，彈性模量E對于C、φ值反演結果的影響較小。在圖12（a）中，當彈性模量E的變化范圍為135 MPa 時，C值的變化大約僅有1 MPa，而φ值的變化約3°。在圖13（a）中，當彈性模量E的變化范圍為45 MPa 時，C值的變化不到2 kPa，而φ值的變化小于1°。這就是說，C、φ值的反演結果對彈性模量E不敏感。相反地，由圖12（b）和圖13（b）可見，鎖定體積應變η*對于C、φ值的反演結果卻有很大影響，即C、φ值的反演結果對參數η*是敏感的。

圖13 利用貫入深度-初速度數據［31］進行反演時，參數E、η*對C、φ值反演的影響，其中αs=276.92 N、βs=0.30 kg/mFig.13 Influence of parametersEandη* on the inversion results ofCandφ，respectively，whereαs=276.92 N andβs=0.30 kg/m（using the penetration depth-initial velocity data［31］）

此外，從圖12（b）和圖13（b）還可以看到，不同阻力系數導致根據相同范圍內的η*反演得到的C、φ值也不相同。這說明反演結果的參數敏感性也與阻力系數αs、βs有關。

5 結論

本文基于動態空腔膨脹模型，利用土壤沖擊貫入實驗中測得的阻力-速度關系、加速度時程曲線、速度時程曲線，以及貫入深度-初速度關系，建立了對土壤Mohr-Coulomb 參數（內聚力C和內摩擦角φ）進行反演的方法。其中，基于前三者建立的反演方法是利用單次貫入過程中的時程曲線進行的反演；而基于貫入深度-初速度關系建立的反演方法則是在對同一狀態的土壤進行多次貫入時采用。

在貫入過程中的穩定貫入階段，貫入阻力與貫入速度平方之間呈線性關系，準確地擬合線性系數（阻力系數）αs和βs是反演Mohr-Coulomb 參數的關鍵。本文揭示了利用動態空腔膨脹模型進行Mohr-Coulomb 參數反演時存在一個關于參數的可解域。對于測得的實驗數據，只有首先保證擬合得到的阻力系數落在可解域范圍內，C、φ值的反演才能在動態空腔膨脹模型的框架下進行。

本文分析表明利用阻力-速度關系進行反演的結果與利用加速度時程曲線進行反演的結果相近，而這兩者與利用速度時程曲線進行反演的結果差別較大。這本質上是由于阻力系數αs、βs分別對加速度和速度的貢獻不同。因此，這三種反演方法對C、φ值反演的準確性也是不同的。其中，利用速度時程曲線進行參數反演的結果最為準確。利用Forrestal 等人［25］報道的土壤原位貫入實驗結果進行的方法驗證表明，利用速度時程曲線對土壤內聚力反演的相對誤差為2.14%，對土壤內摩擦角反演的相對誤差為9.77%。

此外，Mohr-Coulomb 參數的反演具有參數敏感性，這不僅與材料參數本身有關還與阻力系數αs和βs有關。C、φ值的反演結果對鎖定體積應變η*敏感，而對彈性模量E不敏感。

本文建立的反演方法可為復雜地質環境下土壤Mohr-Coulomb 參數的快速確定提供基礎，并為大范圍內的地質力學信息勘測提供一個新的途徑。