地形變旋轉張量探討＊

劉序儼 黃聲明 林巖釗

(福建省地震局,福州 350003)

地形變旋轉張量探討＊

劉序儼 黃聲明 林巖釗

(福建省地震局,福州 350003)

在分析地形變組成的基礎上,指出旋轉與應變張量矩陣同出而異名,兩者皆為位移梯度矩陣的不同組合。旋轉張量矩陣為位移梯度矩陣與其轉置矩陣之差的 1/2,為一反對稱矩陣,由于該矩陣的對角線元素為零,旋轉張量矩陣已退化為一向量,該向量的方向與模表征了地殼的剛體旋轉角位移。給出了位移梯度矩陣、旋轉與應變張量矩陣的普適表達式,而旋轉張量矩陣的諸元素皆為位移分量偏導數的減組合,要想根據旋轉張量矩陣計算旋轉角位移,必須首先求得位移分量的坐標函數式,且計算得到的旋轉角位移又不是一直接觀測量,為了克服求取位移分量坐標擬合式的困難,又能對地表旋轉成分進行直接觀測,提出了一種簡單且行之有效的觀測與計算方法。

地殼形變;位移梯度矩陣;旋轉張量矩陣;旋轉觀測;計算方法

1 引言

對于前者,dr與 dr′之間關系可用下式表示:

對于后者,顯然式 (1)是不適用于圖 2所示情況,此時須用一個地形變張量矩陣取而代之,對圖 2中的情況[1,2]:

式中,F為地形變張量矩陣,對于二維平面情況,F為二階矩陣,對于三維空間,F為 3階矩陣,根據文獻[1,2],任何一種地形變可由平移、旋轉和應變 3部分組成,因此,式(2)中的地形變張量矩陣 F又可分解為：

式中,I為單位矩陣,表示剛體平移;Ω為旋轉張量矩陣,表示剛體旋轉;ε為應變張量矩陣,表示地塊形狀的改變。將式(3)代入式(2),則得表示地形變的普適表達式：

式(4)的幾何含義就是如圖 2的那樣對 dr進行以下操作:先把 dr平移使其A與A′重合,然后把 dr旋轉一個β角,使 dr與 dr′重合,最后使 dr伸縮,使 dr等于 dr′。同樣,這種操作也適用于平面上的正方形和空間中的立方體,不過除了平移、旋轉和伸縮外,正方形和立方體的相鄰兩直角邊的夾角也可能發生改變,即剪應變。既然地形變是由平移、旋轉和應變 3部分組成,宛如顏料中的紅、黃、藍 3種原色可以調出百色一樣,平移、旋轉與應變也是地殼變形的 3種“原色”,它們表征了地殼質點在外力激勵下做出的響應。

圖1 平行情況Fig.1 Parallel case

圖2 不平行情況Fig.2 Unparallel case

旋轉與應變皆由質點位移所引起,它們分別為位移梯度矩陣的不同組合,旋轉與應變張量矩陣之和即等于位移梯度矩陣。設其質點位移向量為 u,E為位移梯度矩陣,則有：

2 位移梯度矩陣及其性質

設M與N為在ξ坐標系中相鄰近的兩點。在外力作用下,M移至M′,N移至 N′。M的位移為u(r),N的位移為 u(r+Δr),其中,Δr=(ΔS1,ΔS2, ΔS3)T為M點的位置向量的增量(圖 3)。

圖3 地殼質點位移示意圖Fig.3 Sketch of crustal particle displacement

在圖 3中,由于M與 N是非常接近的兩個點,對M點位移所引起的N點位移向量 u(r+dr)進行臺勞級數展開,僅取一次項可得:

由上式又可得:

根據向量值函數導數定理[3],由式 (5)可得:

式中,ΔU1、ΔU2、ΔU3分別為位移向量 u的全微分在M點單位活動標架上的分量,則 E的表達式為：

E稱為M點處的位移梯度矩陣或位移向量 u的雅可比矩陣。由張量理論,位移向量為一階張量,位移梯度矩陣為二階張量矩陣。

位移梯度矩陣 E給出了位移向量的空間變化率,刻畫了地殼介質的形變。這種“形變”體現在介質內部質量體積元位移的相互作用上,這種介質內部相互作用產生的“形變”,在彈性力學中稱為“應變”,位移的空間變化率即為應變[4]。

由M點的位移梯度矩陣 E及轉置矩陣 ET所構成的兩種特定組合,分別給出由M點位移在其鄰域所產生的地殼介質的應變和旋轉張量矩陣[1,2],前者可表達為:

式中,

后者可表示為:

式中,

其中,

前者的對角線元素分別表示在M點處正交坐標架軸上的線應變,而非對角線諸元素則分別表示在M點處由彼此相互正交的坐標軸所構成的平面上的剪應變,在三維空間中,應變張量矩陣僅有 6個獨立變量,在二維情況,僅有 3個獨立變量。雖然在不同的正交曲線坐標系下應變張量矩陣的元素各不相同,但它們都是相似矩陣,都能刻畫在M點處的所有地形變信息。無論由哪一種坐標系下的應變張量矩陣均能求出M點處的地形變不變量,這些不變量與坐標系的選擇無關,這些不變量隱藏在該矩陣之中。它們分別為該質點處的主應變及其主方向,在由彼此正交的兩個主方向構成的平面上的剪應變為零,主方向上的應變即為主應變,主應變之和等于該矩陣的跡,為該質點處的體應變,矩陣的行列式等于主應變的乘積,該矩陣皆可由主方向的單位向量所構成的變換矩陣轉換為以主應變為對角線元素的對角矩陣,主應變與主方向皆可由矩陣的特征值方程求得,以上這些不變量正體現了應變張量矩陣的幾何物理性質。對于后者,由于為一反對稱矩陣,此時該矩陣已通化為一向量,該向量表征了該質點處的剛體旋轉方向,而該向量的模則為剛體旋轉的角位移。M點處的地形變應變與旋轉成分是由M點處的位移空間形變率即位移梯度所產生的,因此M點處的應變與旋轉則表征了位移梯度矩陣的性質。

3 旋轉張量矩陣及其性質

從式 (13)可看出,旋轉張量矩陣為反對稱矩陣,因此式(13)中的對角線分量皆為零,此時,旋轉張量已退化為一個旋轉向量。式 (13)中各分量ωij(i≠j)表示為位移所引起的旋轉向量的模在 ij平面上的分量,具有明確的幾何意義。例如,ω12表示在M點處的單位標架中由 e1和 e2單位向量所構成的平面上(以下稱為 12平面)所發生的角度轉動量(以弧度為單位),在單位時間的轉動量稱為角速度。ω12值的正負由下述約定給出:如果ω12為正值,那么,根據右手法則,四指按反時針由 e1轉向 e2,大拇指所指的方向就代表在 12平面上所發生的轉動向量的方向[5],該轉動的大小與方向可表示為ω12e3;反之,ω21代表在 12平面上按順時針方向轉動,轉動向量的方向指向 -e3,該轉動的大小和方向可表示為 -ω21e3,因ω21=-ω12,則ω21e3=-ω12e3,說明 12平面上的轉動向量與約定的時針方向無關,其他兩個平面上的旋轉向量亦可依此類推。最后可得由式(13)所表示的旋轉向量ω表達式為

設:cosα、cosβ與 cosγ分別為向量ω的方向余弦,則ω的單位向量為

式中,

旋轉向量ω又可表示為

式中：ω0表示由位移所引起的單位轉動向量;ω表示由位移所引起的與ω0相正交平面上的旋轉量,其正負根據右手法則給出;ω23、ω13、ω12分別為ω在23、13、12平面上的分量,而ω本身也就是旋轉向量的模。

不同坐標系下的旋轉張量矩陣亦如應變張量矩陣一樣,也可通過這兩個坐標系的轉換矩陣 c進行轉換,其表達式為[1]:

式中,c=cB,其中 cB與 cA分別為由M點處的直角坐標系的單位標架轉換到該點處的 B與 A坐標系的單位標架的轉換表達式[5],且 ccT=I,故旋轉張量矩陣如同應變張量矩陣一樣也是相似矩陣,旋轉分量與的轉換亦同應變張量分量的轉換一樣是可轉換的。

考慮到旋轉張量矩陣ω為一個由式 (13)所表示的一個已退化的轉動向量,當然ω =(ω23,ω23, ω23)亦可由坐標系轉換公式得到以下相互轉換表達式:

式中,c=cA由于 c為正交矩陣,ccT=I,不難證明：

故不同坐標系下的旋轉張量都能給出該質點處地殼整體旋轉方向及其旋轉角位移,這正是該矩陣所隱藏的幾何物理不變量。

顧及到式 (15),則式 (23)又可寫為：

4 地表面旋轉張量的運動方程

旋轉張量矩陣與應變張量矩陣都是同出而異名,即它們都是位移梯度矩陣的不同組合。前者為減組合,后者為加組合。根據數字濾波理論,前者為高通濾波器,后者為低通濾波器。在無地震情況下,地殼噪聲一般表現為地脈動,其振幅都是較小的,因此,地殼旋轉分量的數值一般比平移和應變要小。

地殼旋轉如同應變一樣,都是由地殼剪切應力所引起的。在這里,我們僅討論地表旋轉問題,因為我們目前還沒有這樣一個三維形變觀測系統。根據牛頓第二定律與物體轉動的力矩平衡原理,繞垂直于地表面的 Z軸旋轉的旋轉運動方程為[4]:

式中,θz是繞 Z軸的轉動角,Iz為地殼體積微元dxdydz繞 Z軸的轉動慣量,Gz為體積力力矩密度的Z軸分量,Gzdxdydz為該 Z軸分量繞 Z軸的旋轉力矩,Iz=ρdxdydz[(dx)2+(dy)2]平均,其中ρ為地殼密度,σxy與σyx分別為地表面上的 2個剪應力。當體積微元 dxdydz趨于零時,轉動慣量 Iz趨于零的速度更快,因而上式變為:

在小應力情況下,上式又可近似寫成:

同理,可得繞 X軸與 Y軸轉動的運動方程:

式中,Gx與 Gy分別為體積力力矩密度的 X軸與 Y軸分量繞 X軸與 Y軸的旋轉力矩。

當上式 Gi(i=x,y,z)=0,則得到剪切應力張量的重要關系式:

文獻[7,8]則直接從力矩平衡方程出發,在假定Δxi(i=1,2,3)→0前提下,同樣得到了σij=σji。式(30)表明剪切應力是對稱的。在線性小應變彈性波傳播理論中,Gi=0是允許的。應力張量的對稱性“消隱”了轉動運動方程存在的形式,因此,在一般文獻中對轉動運動方程都不予討論[4]。張量的對稱性本身就呈現了“旋轉成分”。

既然旋轉為組成地形變 3種成分中的一種,為一種獨立成分,據筆者所知,到目前為止,雖然還沒有一種專門的觀測系統對其進行觀測研究,但旋轉作為地形變的一種獨立成分,如何對它進行觀測一直成為學者關注的焦點,其中,地震旋轉儀的研制與試測就是顯著的一例[9]。該旋轉儀的本體為一完全對稱的圓環形旋轉慣性擺。雖然該儀器經過 1年的試測記錄到了一些地震旋轉波,但因地震位移波與旋轉波同時到達,位移波會對旋轉波產生極強的干擾,使旋轉波被淹沒在位移波中難以分辨,如何抑制地震 P波波動使之不會對旋轉波產生影響,是一項非常關鍵的而又十分不易解決的難題,終因種種原因遭遇到重重困難,最后也只好將試測樣機束之高閣而作罷,但不管怎樣,研制者卻開創了中國地震旋轉儀研制的先河。迄今為止,還沒有找到一種行之有效的觀測系統與計算方法,在這里,本文提出一種簡單易行的觀測方案與計算方法。

5 旋轉張量計算面臨的困境

在取得諸測站的位移觀測資料后,是否就可直接按式(11)與(13)進行地應變與旋轉形變分析呢?回答是否定的,其原因是我們無法知道位移的坐標函數表達式,因此也就無法得到位移諸分量對該測站坐標的偏導數值。根據彈性力學理論基礎,地殼介質這種彈性體中所有的物理量都是連續的,即是說,密度、位移、應變、應力都被假定為空間點的連續變量,并且假定空間點變形前后應該是一一對應的,位移向量為空間點的單值函數并具有所需的各階連續偏導數[1],在一般情況下,我們無法知道位移的坐標函數表達式,但天體起潮力引起的測站位移是個例外。根據固體潮理論,在采用古登堡-布倫地球模型的基礎上,天體起潮力引起的位移可根據該測站處的天體起潮力計算出來,從而可得到天體起潮力所導致的測站位移在球坐標系下以測站地心向徑ρ、地心緯度φ與經度λ為變量的位移坐標函數表達式,最后根據式(11)與(13)進行應變與旋轉張量計算,以得到固體地球應變與旋轉張量的理論值[10-12],把應變固體潮觀測值進行維尼迪柯夫調和分析,就可在振幅與相位方面對兩者進行分析,以取得可靠的地形變信息。

在一般不知位移坐標函數表達式的情況下,如何解決應變張量計算呢?為此,許多學者進行了探討,例如,文獻[13]建議以若干個觀測點作為支撐點來確定某個內插函數,在進行內插時,如果沒有給觀測值提供一種假定的動力學模型數據的話,采用純粹的幾何模型內插方法必須滿足觀測點的分布密度能保證相鄰觀測點之間的線性內插達到足夠的精度以滿足研究的需要,以及所選定的幾何內插函數能確保相鄰觀測點之間存在近似線性內插關系,內插函數可采用樣條逼近技術或擬合推估方法求得。文獻[14]認為,在我們所得到的是一些離散點的位移,而不是點位的函數時,則無法解析地求出偏導數,此時必須用數值估計,有兩種方法可供采用。一種是求一個小區域的導數,此時需要采用有限元方法,另一種是求一個點上的導數,此時可用有限差分法。因此,利用應變張量矩陣公式計算應變張量遇到了建立位移分量的坐標線性擬合公式這個難題。在利用旋轉張量矩陣公式計算旋轉張量時也遇到了同樣的難題,如何繞開這個難題,直接利用位移觀測資料進行旋轉張量計算只有另覓他途。

6 旋轉觀測與計算方法探討

地殼旋轉作為地形變的一種獨立成分,除了文獻[1,2]從理論上給予證明以外,更主要的是要以觀測事實予以證明,因為自然科學不像思辨哲學那樣僅是在形而上學層面上進行思考,而自然科學本身是一種實證科學,對任何一種理論都必須經過觀測或實驗加以驗證。如何對地形變旋轉成分進行觀測,本文提出一種觀測方案 (圖 4)。圖中OA與OB相互垂直,分別位于 x軸與 y軸上,OA=OB,分別在A、O、B進行位移觀測。

圖4 旋轉觀測示意圖Fig.4 Sketch of rotation observation

圖 4中Ux(·)與 Uy(·)分別為 A、O、B點的位移分量,其中,Ux(O)與 Ux(A)的差異會使 OA發生脹縮,而Uy(O)與 Uy(A)的差異則會使 OA發生偏轉,類似地Uy(O)與Uy(B)之差異會使OB脹縮, Ux(O)與Ux(B)會使 OB發生偏轉。圖中 OA邊旋轉角位移向量ωOA為:

類似地,OB邊的旋轉角位移向量ωOB為:

以上計算是依照右手法則來確定旋轉向量方向的,旋轉方向定為為與地表垂直向上方向的一個單位向量,如果式 (31)與 (32)中右邊的系數為正,表明旋轉為逆時針,反之為順時針方向。最后,兩邊位移所引起的O點處的地殼旋轉位移ω為上述兩式之平均值：

上式正是由式 (14)所表示的二維地表應變張量矩陣中的ωyx的旋轉向量表達式。

采取這樣的方案不但可以對地表旋轉進行觀測,而且可以進行計算,從而可以彌補這方面的空白。

至于三維地殼旋轉張量觀測,我們必須在空間6個不同方向進行位移觀測。目前,我們還不可能在地殼洞體中建立這樣一個觀測系統,因此,也就無法取得其他兩個旋轉張量。如何對地殼旋轉進行完備觀測與計算,仍然是一個今后值得重視的問題。

在這里,值得提出的是,本文給出的地形變旋轉張量觀測與計算公式(34),是受旋轉張量式 (14)的啟發,式(34)不過是式 (14)在圖 4中等邊直角三角形上的近似展開罷了,觀測與計算原理是正確的。由于福建省地震局 GPS網點布設時,沒有考慮地旋轉觀測方案,因此,無法以觀測資料實例對該方法的可行性與實效性進行驗證。不過我們可以對式(34)計算地旋轉張量的精度進行估計。

設:Δx=Δy=l,設 GPS單點觀測位移的中誤差為σu,設σω為ω的計算中誤差,根據誤差傳播定律,不難得出:

在目前 GPS單點觀測精度最高可達σu=1 mm,設 l=100 m=105mm,則按式(35),此時,σω= 10-5(弧度)=2″,如取 l=1 000 m=106mm,則σω=0.2″,如 l=10 km=107mm,則σω=0.02″,此時ω的測定中誤差相當于傾斜固體潮最大幅度。因此,采用此種方法是無法勝任旋轉固體潮觀測的,但可以滿足對大震的同震旋轉觀測,此種旋轉觀測系統是可以對大震所引起的同震旋轉作出響應的,因為大震的同震旋轉量較大。

7 認識與討論

文獻[4]認為應力矩陣的對稱性“消隱”了應力轉動方程存在的形式,但不是意味著旋轉成分的消失,而是認為應力張量的對稱性本身就呈現了“旋轉成分”,旋轉成分的存在是一種客觀事實,不過其量值比正應變與剪應變值要小罷了,大地震所造成的煙筒或鐵軌扭曲就是旋轉成分存在的一種例證。究其原因,旋轉是由不同距離處的質點位移之差異所造成的,即是由質點的不同空間形變率所造成的。本文從地形變構成出發,闡述了應變與旋轉張量矩陣皆為位移梯度矩陣的不同組合,并指出位移梯度矩陣刻畫出了空間形變率,并給出了應變與旋轉張量矩陣的普適表達式。由式(11)與 (13)可發現,要想根據位移觀測資料進行應變與旋轉張量分析,首先要給出位移分量的坐標函數式,除了應變固體潮外,我們一般是無法獲知位移分量的坐標函數的表達式,此時要從建立動力學或幾何模型并采用線性內插方法來取得位移分量的坐標擬合式,在此基礎上,才能根據式(11)與(13)進行應變與旋轉張量分析[10,11],從這個意義上說,根據式 (13)對位移分量的坐標擬合式取偏導而取得的旋轉并不是一種直接觀測量,而且存在擬合誤差帶來的影響。對此,本文給出了一種繞開式(13)進行旋轉分析的方法,提出了一種在等邊直角三角形 3個頂點進行位移觀測及其計算方法,該方法簡單易行,可避免采用式 (13)要確定位移分量坐標函數式的困難。在目前的地形變觀測系統中,唯缺失旋轉觀測系統,只有建立了旋轉觀測系統,我們才真正具有了對地形變 3種成分進行觀測的一個完備的觀測系統。

1 王敏中,等.彈性力學教程[M].北京：北京大學出版社, 2002.(WangMinzhong,et al.The course of elastic mechanics[M].Beijing：BeijingUniversity Publishing house,2002)

2 米恩斯W D.應力與應變[M].丁中一譯,王仁校.北京：科學出版社,1982.(Means W D.Stress and strain[M]. Springer-VerlayNew York,Inc,1976)

3 楊永發,徐勇.向量分析與場論[M].天津：南開大學出版社,2006.(Yang Yongfa and Xu Yong.Analysis of vector and field theory[M].Tianjin：Nankai University Press, 2006)

4 牛濱華,孫春巖.固體彈性介質與地震波傳播[M].北京：地質出版社,2005.(Niu Binhua and Sun Chunyan.Solid elastic medium and seismic wave trans mit[M].Beijing：Geologic Press,2005)

5 劉序儼,等.正交曲線坐標系的應變張量轉換[J].大地測量與地球動力學,2008,(2)：71-76.(Liu Xuyan,et al. Conversion of strain tensor matrices between t wo orthogonal curvilinear coordinates[J].Journal of Geodesy and Geodynamics,2008,(2)：71-76)

6 GreenbergM D.Advanced engineeringmathematics(The 2nd edition)[M].Beijing：Publishing House of Electronics Industry,2004.

7 Lay T,et al.Modern global seis mology[M].San Diego,New York：Boston：Academic press,1995.

8 Jaeger J C.Elasticity fracture and flow[M].Lodon：Methuen amp;CO.LTO.New york：John wi-leyamp;sows.inc,1964.

9 蔡乃成,付子忠.旋轉地震儀的研制[J].地震學報,2009, 31(3)：347-352.(Cai Naicheng and Fu Zizhon.Manufacture of rotation seis mograph[J].Acta Seismological Sinica, 2009,31(3)：347-352)

10 北京大學地球物理系,等.重力與固體潮教程[M].北京：地震出版社,1982.(Depart ment of Geophysics of BeijingUnivisity,et al.Course of gravity and earth tide[M]. Beijing：Seis mological Press,1982)

11 梅爾基奧爾 P.行星地球的固體潮[M].杜品仁,等譯.北京：科學出版社,1984.(Melchior P.The tides of the planet earth[M].Translated byDu Pinren,et al.Beijing：Science Press,1984)

12 Vanicek P and Krakiwsky E J.Geodesy：The concepts[M]. New York：Elsevier Science Publishing company,I NC. 1986.

13 Altiner Y.Analytical surface deformation theory for detection of the earth’s crust movement[M].Berlin.Heidelberg.etc.Springer-Verlag,1999.

14 瓦尼切克 P.四維大地測量定位[M].黃立人,等譯,陳鑫連校.北京：地震出版社,1990.(Vanicek P.Four-dimensional geodetic positioning[M].Springer International 1987)

RESEARCH ON ROTATION TENSOR OF CRUSTAL DEFORMATION

Liu Xuyan,Huang Shengming and Lin Yanzhao
(Earthquake Adm inistration of Fujian Provine,Fuhzou 350003)

On the basis of the analysis of component combination of crustal deformation,it is pointed out that strain and rotation tensor matries both are different combination of displacement gradient matrix,and the rotation tensormatrix is a half of diffrence of both displacement gradieutmatrix and it is transpose one and anti-symmetrix one.Because its diagonal elements of rotation tensormatrix is zero so it turns into a kind of vector in degradation, the direction and modulus of the vector characterize the crustal rigid rotation.

crustal deformation;displacement gradieatmatrix;rotation tensormatix;rotation observation;calculation of rotation

1671-5942(2010)05-0057-07

2010-06-21

中國地震局老專家基金

劉序儼,男,1939年,研究員,長期從事固體潮與地殼形變研究.E-mail：xuyanliu@126.com

P315.72+5

A