各向异性是地壳内部普遍存在的一种结构特征，地壳各向异性结构特征的研究，对地下结构的应力分布、介质物性特征、块体运动方向、地质构造背景及地震预测(郑秀芬等，2008)的研究都具有重要意义.

地壳各向异性的地震学研究方法主要有横波分裂法(刘希强等，1998)、Pn波法、面波法、噪声成像法(王琼和高原，2012)、反射波法以及接收函数法(徐震等，2006).其中Pn波法是最早开始研究的，它的优势在于垂向分辨率好，Hess(1964)曾用Pn波法成功解释了洋底存在各向异性.横波分裂法能够精确地给出平均各向异性程度，在各向异性的横向分辨率上具有明显优势，目前该方法已在国内外广泛使用(刘希强等，1998；罗艳等，2004；高祥林，2004；常利军等，2006，2008；吴晶等，2007；张建利等，2012；张辉等，2012).接收函数方法通过分离地壳、上地幔中相关界面的震相对地球内部结构进行研究，特别是当地壳和上地幔中存在各向异性的分层结构时，在接收函数的切向分量具有较为明显的响应.近年来，随着地震观测台站的迅速增加，接收函数的分辨能力和稳定性也越来越高，这为用接收函数做各向异性反演提供了条件.

使用接收函数研究各向异性主要有三种方法：

第一，沿袭剪切波分裂的做法，计算moho面Ps转换波的偏振方向和分裂时间;

第二，通过加权叠加的方法突出各向异性特征，直接找出偏振方向;

第三，利用反演搜索，找到各向异性弹性参数.

第一种方法的优点在于方法稳定可靠，已有很好的验证，不足点在于不能得到分层各向异性，而第二种方法的优点在于计算快速简便，但缺点是只适用于各向异性对称轴水平的情况，因此本文按照第三种方法的思路，调整反演策略及先验约束来改善反演的效果.

目前，国内外利用接收函数研究各向异性已有不少实践.Levin和Park(1997)等利用远震接收函数研究了乌拉尔地槽的地壳各向异性.Savage(1998)用理论接收函数研究了新西兰SN20台站下地壳各向异性和倾斜界面特征.Zandt、Sherrington、Ozacar等(Sherrington et al.，2004; Ozacar and Zandt, 2004; Frederiksen et al.，2003)利用接收函数研究了青藏高原下方地壳内部各向异性的分布特征.Peng和Humphreys(1997)利用理论接收函数对内华达西北地区的Moho倾斜和地壳各向异性进行了研究. 徐震等(2006)利用接收函数资料，解释了哀牢山—红河断裂带地区地壳的各向异性 .田宝峰等(2008)利用模拟退火算法反演了华北太行山区地壳的各向异性.

本文首先对影响接收函数的因素做了单变量分析，系统总结了每个参数对接收函数的影响规律.然后结合邻域算法(Neighborhood Algorithm)对理论接收函数做了反演数值试验，并采用不同失配函数进行了对比测试.数值试验结果表明对复杂模型需要采取分步反演的策略才能得到预期的效果.此外，还根据反演的数值试验确定了邻域算法参数的选取.最后，本文以青藏高原东北缘三个台站实际的接收函数资料为例，对本文的反演方法和程序进行了应用.

六方晶形对称结构有一个唯一快轴方向和慢轴方向，并且具有独立的快慢轴速度(Weiss et al.，1999).体现在弹性系数矩阵上面需要五个独立的参数.Levin和Park(1997，1998)根据Backus(1965)的推导给出了这五个弹性参数的定义为

远震P波传播到台站下方时可以近似看作是平面波，在遇到台站下方的间断面入射时，可以得到反透射系数的δ脉冲序列，因此我们使用δ脉冲作为震源计算理论接收函数.由于介质各向异性结构对接收函数的影响主要反映在不同方位上，即接收函数按照不同方位角的排列可以直观反映各向异性的特征，因此在正演模拟中取定射线参数，而计算不同方位角的理论接收函数.

为了快速有效的模拟接收函数，本文采用Frederiksen等(2003)提出的基于震相追踪的射线方法.该方法基于层状模型，在同一层中同时考虑介质倾斜(包括走向strike和倾角dip)、各向异性程度(pctP、pctS)、对称轴与水平面的夹角(Plunge)以及其水平面投影与正北方向夹角(trend)，如图 1所示.

图 1

Fig. 1

图 1 各向异性快轴分解示意图 OF为实际各向异性快轴(或者慢轴)方向，OG为OF在水平面上的投影.ON代表正北方向，OE代表正东方向.其中，plunge为水平投影(OG)到OF的夹角，trend为正北方向(ON)到OG的夹角. Fig. 1 Anisotropy decomposition schematic diagram OF is the fast(or slow)axes of anisotropy, OG is the projection of OF onto a horizontal plane. ON is the North orientation, OE is the East orientation. Plunge is the angle between OG and OF, trend is the angle between ON and OG.

该方法采用震相追踪，可以方便的对所研究的问题中的震相进行个性化定制，排除不必要的震相干扰.因此具有快捷的计算速度，在应用于反演中时可大大节约时间.

目前波形反演中使用的计算失配函数的方法主要有对数法、互相关法和范数法.对数法主要用于震相振幅的比较；互相关法适用于形状匹配；而范数法适用于整体的比对，结合不同的范数类型，可以得到不同的效果.另外，在这些基本方法上结合不同的数学变换，如希尔伯特变换、小波变换等，可以衍生出许多新方法(Bozdağ et al.，2011).但是考虑到我们研究的问题以及所使用的正演方法，本研究中只使用了互相关法和范数法计算适配函数.具体的公式为

计算出失配函数的值，就可以调用反演算法来迭代求解.波形反演是典型的非线性问题，常用的方法有遗传算法、模拟退火算法、邻域算法以及人工神经网络.Sambridge(1999)等人提出的邻域算法在上述反演方法中具有其独特的优势.该算法在模型空间内随机采样，然后对模型空间进行剖分，利用迭代搜索的方法确定最优解.该方法具有优异的拟合模型空间的能力，当把该方法和模拟退火或遗传算法结合时可以加速收敛.而单独使用该方法，全局收敛的效果更佳.不同于模拟退火和遗传算法，邻域算法跳出局域最优的解决办法不依赖于一些概率常数，如遗传算法里的变异率和交换率.另外，该算法是自适应的，能够自动的跳出局部最优点.目前虽然没有一种算法能够完全摆脱局部最优而找到全局最优解，但邻域算法的这种自适应性在此方面上有着直观的效果.这也是本研究采用该算法的一个原因.此外，该搜索算法具备很好的并行性，Rickwood和Sambridge(2006)发布了该方法的并行反演代码，这也为我们的工作提供了较方便的实现途径.

在实际资料中，造成接收函数切向分量的能量不为零的原因主要有各向异性、倾斜界面、横向不均匀和仪器故障以及仪器安放错误.排除安装和仪器故障，其余都是由介质性质引起的.当模型基于层状介质假设时难以反演横向不均匀性，但是由于台站下方的区域相比整个板块来说很小，这样每一层的横向不均匀就可以忽略，即认为是均匀的.因此，本文仅研究各向异性和倾斜界面的影响.为此，我们做了如下的正演模拟数值试验.

倾角和走向是描述界面倾斜的基本物理量，我们分别对这两个量在接收函数上的影响进行测试.

我们在两层各向同性模型中的第二层加入不同程度的倾斜，计算了不同方位角的理论接收函数，结果如图 2所示.结果表明当倾角发生正负变化时(即倾斜面沿走向方向发生反转)，R分量图案将关于后方位角180度处(由走向决定)镜像反转，T分量在反转的基础上反号.当倾角不变号时，R分量图案变化较大，主要体现在具体震相幅度值正负反转的位置(后方位角).T分量图案变化虽然不大，却很有特征，呈现两段式分布.

同上所述，我们在两层各向同性模型中的第二层加入10度倾斜，然后改变走向.从图 3可以看出，图案随走向的变化存在周期性的规律.无论是R分量还是T分量都只有明显的移动而无振幅的变化.它们随走向每增加1度，朝后方位角变大的方向移动1度，当超过360度时会回到0度.显而易见，这是因为走向会改变由倾角造成的入射角变化的空间排列.

各向异性的对称轴可以分解成trend和plunge两个角度，如图 1所示.

我们在水平两层模型中的第二层加入P和S波各5%的各向异性，其他层保持各向同性，假设Plunge为45度，然后改变trend的大小.从图 4的结果可以看出，R分量和T分量的基本图案不变，两个分量的图案都会随着trend夹角的变化而移动.移动规律都是随trend每增加1度，就向后方位角增大的方向移动1度，到达360度时就从0重新开始.

同样，我们在水平两层模型中的第二层加入P和S波各5%的各向异性，其他保持各向同性，假设trend为0度，然后改变plunge的大小.那么，从图 5可以看出，图案的变化是巨大的，T分量尤为显著.当plunge接近0度时，T分量图案呈四段式分布，但当plunge向90度变大时，图案的四段式分布被破坏，几个直达波明显的向两段式过度，而多次波还能维持四段式分布.但当plunge为90度时，各向异性的特征消失，表现的如同各向同性一般，该特征也体现在了R分量上.

根据以上数值试验结果，我们得到各参数对接收函数影响的如下规律：

(1)从上面的数值研究可以看出在相同的层厚、速度和密度情况下，决定图案总体形状的是介质分界面的倾角和各向异性轴的plunge, 且plunge起到的作用会比倾角更大；

(2)介质分界面的倾斜走向和各向异性轴的trend只会在其他条件相同的情况下影响图案的排列(移动)顺序，而不会从根本上改变图案的性质；

(3)当倾斜和各向异性混叠在一起时，问题就变得非常棘手，甚至有些情况下已经不能直接从图案中找到先验信息.但是我们依然希望从中发现一些有用的信息，这就取决于我们的先验信息和反演手段；

(4)当plunge接近90度时，各向异性特征消失.

通过以上正演分析，可以看出当同时存在倾斜和各向异性时，很难直接从观测资料中确定介质各向异性的参数.此时，反演是确定各向异性结构的一种有效选择.下面分三种情况来测试反演的有效性，进而确定反演时应采取的策略.

图 2

Fig. 2

图 2 界面倾角(dip)对接收函数的影响 图中第一行为R分量，第二行为T分量，两分量从左到右分别对应倾角-30、-20、-10、10、20、30度. Fig. 2 The influence of interface dip to the receiver function The first row is the R components, and the second row is the T components. From left to right in each row, the dip take -30，-20，-10，10，20，30 degrees.

图 3

Fig. 3

图 3 界面倾斜走向(strike, 向东倾斜为0度，向西倾斜为180度)对接收函数的影响 图中第一行为R分量，第二行为T分量，两分量从左到右分别对应strike为0、60、120、180、240、300度 Fig. 3 The influence of interface strike to the receiver function The first row is the R/ components, and the second row is the T components. From left to right in each row, the strike take 0，60，120，180，240，300 degrees

图 4

Fig. 4

图 4 trend对接收函数的影响 图中第一行为R分量，第二行为T分量，两分量从左到右分别对应trend的0、60、120、180、240、300度. Fig. 4 The influence of anisotropy trend to the receiver function The first row is the R components, and the second row is the T components. From left to right in each row, the trend take 0，60，120，180，240，300 degrees.

图 5

Fig. 5

图 5 plunge对接收函数的影响 图中第一行为R分量，第二行为T分量，两分量从左到右分别对应plunge的0、20、40、60、80、90度. Fig. 5 The influence of anisotropy plunge on the receiver function The first row is the R components, and the second row is the T components. From left to right in each row, the plunge take 0，20，40，60，80，90 degrees.

图 6

Fig. 6

图 6 两层模型互相关法反演数值试验结果 (a)~(f)分别为P波各向异性程度(pctP)、S波各向异性程度(pctS)、trend、plunge、波速比(k)和平均横波速度(vs)随深度的搜索结果，其中每一条颜色的线代表一次搜索得到的最优模型，颜色代表失配函数值(msfit)的大小，黑线代表反演最终得到的模型.(g)为失配函数值随迭代次数的变化，横坐标为迭代次数，纵坐标为失配函数值(misfit).其中的红线代表每次迭代的失配函数值的最小值，蓝线代表每次迭代的失配函数值的平均值. Fig. 6 Two layer model inversion result with cross correlation method (a)~(f)are search result for P wave degree of anisotropy(pctP)，S wave degree of anisotropy(pctS)，trend, plunge, velocity ratio(k)，average S velocity(vs)，respectively. Each line in the figure st and s for a optimal model and its color st and s for misfit value, black line st and s for the final model.(g)is misfit value change with iteration, x-axis is iteration, y-axis is misfit value. The red line is the minimal value of misfit in each iteration and the blue line is the mean value of misfit in each iteration.

假设初始模型为水平两层模型，如表 1所示，仅在第二层中加入各向异性(包括程度和方向).失配函数则分别使用互相关、归一化一范数、归一化二范数、一范数、二范数方法(其中归一化是为了避免在实际资料中观测的振幅和正演模拟的振幅大小不匹配).反演的参数包括层厚(thick)、平均S波速度(v_s)、波速比(k)、各向异性程度(pctS、pctP)、各向异性轴的方向(trend、plunge).根据前面正演数值试验的研究得知，T对plunge的变化比较敏感.而且T分量可以很好地拟合出trend的变化.因此，此处选用的R、T、Z三分量权重分别为1:5:0.

表 1 (Table 1)

表 1 两层各向异性反演结果 Table 1 Two layer anisotropic inversion results Layer Method Thick(km) k vs(km/s) pctP pctS Trend(°) Plunge(°) 1 0 38.982 1.767 3.618 0 0 0 0 1 41.175 1.757 3.809 0 0 0 0 2 40.104 1.760 3.708 0 0 0 0 3 49.146 1.549 3.709 0 0 0 0 4 46.782 1.586 3.687 0 0 0 0 r 40.000 1.763 3.710 0 0 0 0 2 0 50.788 1.776 4.090 4.635 4.861 89.897 46.840 1 50.436 1.848 4.326 5.243 4.071 89.536 44.003 2 51.071 1.808 4.187 4.796 6.357 89.772 48.315 3 51.862 1.777 4.194 5.214 3.903 90.043 43.187 4 59.528 1.623 4.190 4.450 4.339 89.591 44.982 r 50.000 1.833 4.200 5.000 5.000 90.000 45.000 Layer代表对应的层，Method中0、1、2、3、4分别代表互相关、归一化一范数、归一化二范数、一范数、二范数法，r代表真实值，k为波速比，vs为平均S波速度，pctP、pctS分别代表P和S的各向异性程度，Trend与Plunge合起来表示各向异性快轴的方向

表 1 两层各向异性反演结果 Table 1 Two layer anisotropic inversion results

反演结果如表 1所示.可以清楚地看到，各向异性对称轴方向的反演接近模型预设值，而各向异性程度和层厚的反演相对较差.造成厚度反演相对较差的原因有两点：(1)选用的权重比侧重于T分量，而R分量对层厚的约束更佳；(2)此处使用的是默认的震相列表，多次波仅考虑第一层，其他多次波没有得到计算.但由于本研究更多的是希望得到各向异性的信息，层厚和平均速度可以根据别的方法如接收函数H-K搜索得到，从而在搜索空间内加入限制，或者如后面实例中所做的那样调整震相列表，考虑更多震相.

根据失配函数值变化特征(如图 6g)，可以看出在当前参数下，程序迭代3000步后结果基本收敛，且反演结果稳定.而从模型空间的覆盖情况来看，如图 6(a~f)，也已达到计算的需求.从图 7中我们能够得出，互相关方法在各向异性方面具有明显的最优效果.另外还可以看出，范数法对S波各向异性程度的约束不好.

Erickson(2002)研究表明，当连续的层中加入各向异性，反演(基于遗传算法)结果将会变差.为了检验本方法对复杂模型的反演效果，建立一个三层水平模型，在其中的第二、三层中加入各向异性(如表 2所示).反演参数同3.1.1的两层模型.失配函数此处仅使用互相关、归一化二范数和二范数.而R、T、Z三分量权重依然是1:5:0.

表 2 (Table 2)

表 2 三层各向异性反演结果(表中括号内为二次反演结果) Table 2 Three layer anisotropic inversion results Layer Method Thick(km) k vs(km/s) pctP pctS Trend(°) Plunge(°) 1 0 26.489 1.780 3.938 0 0 0 0 2 25.678 1.732 3.782 0 0 0 0 4 24.995 1.738 3.736 0 0 0 0 r 25.000 1.763 3.710 0 0 0 0 2 0 21.134 1.858 4.492 4.806 (5.015) 8.514 (4.990) 91.140 51.435 2 21.478 1.808 4.151 5.028 (4.988) 8.734 (4.923) 90.596 45.125 4 22.388 1.706 4.272 5.296 4.182 87.135 44.209 r 20.000 1.833 4.200 5.000 5.000 90.000 45.000 3 0 28.791 1.833 4.417 5.049 (5.973) 7.264 (6.086) 29.877 40.215 2 31.752 1.766 4.181 5.468 (6.066) 9.517 (5.935) 29.002 34.706 4 28.430 1.980 4.197 6.257 8.006 26.992 39.029 r 30.000 1.833 4.200 6.000 6.000 30.000 35.000 各参数如表 1所述.

表 2 三层各向异性反演结果(表中括号内为二次反演结果) Table 2 Three layer anisotropic inversion results

从表 2分析可知，三种方法对层厚和对称轴方向的反演都很好，平均S波速的反演结果也在可接受范围内，S波各向异性程度反演的误差较大，P波的各向异性程度接近真值.此外，从图 8分析可以得出互相关方法对各向异性对称轴方向的反演更稳定一些.由于此反演主要关注的是各向异性参 数，但各向异性程度(pctS)反演的依然不理想(如图 8所示)，未能达到本研究的目的，因此我们考虑二次反演.

图 7

Fig. 7

图 7 两层模型反演结果与预设模型参数的相对误差(模型的第二层) 图中不同颜色的条柱分别代表层厚(think)、波速比(k)、平均S波速度(vs)、P波各向异性程度(pctP)、S波各向异性程度(pctS)、trend和plunge与真实值的相对误差，横坐标代表不同的失配函数方法，定义与表 1一致，纵坐标代表相对误差的百分比. Fig. 7 Relative error between two layer result(the second layer) and its initial model The different color bar is different relative error for thick, velocity ratio, average S velocity, P wave degree of anisotropy(pctP)，S wave degree of anisotropy(pctS)，anisotropy trend, anisotropy plunge, respectively. Its x-axis is misfit method(according to table 1)，y-axis is percentage relative error.

图 8

Fig. 8

图 8 三层模型一次反演结果与预设模型参数的相对误差(左图为第二层的结果，右图为第三层的结果) 图中不同颜色的条柱分别代表层厚(think)、波速比(k)、平均S波速度(vs)、P波各向异性程度(pctP)、S波各向异性程度(pctS)、trend和plunge与真实值的相对误差.图中横坐标代表不同的失配函数方法，定义与表 1一致，纵坐标代表相对误差的百分比. Fig. 8 Relative error between three layer result(the left figure is the second layer result and the right figure is the third layerresult) and its initial model The different color bar is different relative error for thick, velocity ratio, average S velocity, P wave degree of anisotropy (pctP)，S wave degree of anisotropy(pctS)，anisotropy trend, anisotropy plunge, respectively. Its x-axis is misfit method(according to table 1)，y-axis is percentage relative error.

图 9

Fig. 9

图 9 台站分布及各向异性快轴反演结果 图中红色三角表示台站，台站所在位置处不同颜色粗的线段代表不同层各向异性快轴的方向，其长度代表S波的各向异性强度.棕色细线代表断层，黑色细线代表块体边界. Fig. 9 Station distribution and anisotropy fast axis reversion result The red triangle is station, the different thick line in the station st and s for anisotropy orientation of different layer, and its length is S wave anisotropy degree. The brown thin line st and s for fault, the black thin line st and s for block boundary.

在同样的模型下，二次反演只考虑各向异性的程度，其他参数不参与反演，在实际处理资料时，可以沿用一次反演的结果.结果如表 2括号中所示，可以看出，我们采用的两种方法都很好的收敛到了真实解.

从前面的正演数值试验得知在各向异性和倾斜混叠的情况下是很难用经验区分出各自信息的.下面采用分步反演的方法，尝试解决该问题.

为了简单起见，只考虑3.1.1中的两层模型，且层厚、平均波速、各向异性程度还有走向均不参加反演，只反演各向异性快轴方向trend、plunge以及界面倾角dip, 其他的参数可在一次反演中确定.

结果如表 3所示.可以看出，存在倾斜的情况下，三种方法依然能够很好的反演出各向异性的对称轴方向，但是对于倾斜的反演结果误差较大.不过可以看到，互相关方法对倾斜还是有很好的分辨能力，当然这牺牲了一点trend的精度，但还在接受范围内.

表 3 (Table 3)

表 3 各向异性与倾斜混叠反演结果 Table 3 Results of anisotropic and inclined mixed stack inversion Layer Method Trend(°) Plunge(°) Dip(°) 2 0 87.888 45.324 6.650 2 90.322 44.557 5.251 4 89.652 46.508 4.777 r 90.000 45.000 7.000

表 3 各向异性与倾斜混叠反演结果 Table 3 Results of anisotropic and inclined mixed stack inversion

青藏高原作为目前板块运动最活跃的区域之一，一直备受地学研究者的关注.根据Clark和Royden(2000)提出的地壳流假说，青藏高原在南部受到印度板块的推挤，在东部分别受到鄂尔多斯块体和四川盆地的阻挡，造成高原隆起，同时地壳流被侧向挤出，在青藏高原的东北缘和东南缘形成比较发育的地壳流，这种构造背景下可能造成地壳中各向异性结构的存在.本文将上述反演方法，应用于青藏高原东北缘三个固定台站(如图 9所示)，研究各台站下方地壳各向异性结构特征.

通过实际资料的分析，T分量接收函数随后方位角变化呈现出明显的规律性，说明台站下方地壳中存在有各向异性或倾斜的结构特征.为了使反演结果更具客观性，参数的搜索边界选定一个较宽的范围.将初始搜索模型构建成三层各向异性结构，前两层的厚度搜索设在10~20 km, 第三层设在10~25 km, 地壳最大厚度的搜索边界可达65 km.如果三层均考虑各向异性和倾斜，这样反演问题就变得过于复杂，反演空间的维数(参数个数)过大，Tarantola(2004)曾讨论过，反演参数过多时，在反演迭代次数不能足够大时，高维反演空间可能是空的，这可能导致反演结果的不准确性.因此，本文采用分步反演的策略，在第一次反演中不考虑各向异性的程度(设为5%)，这样每一层参加反演的参数就减少到6个，分别是层厚、trend、plunge、倾角、波速比(k)还有平均S波速度(v_s).每一个参数都设在一个较宽的搜索范围内，以确保反演的客观性.

为了提高正演结果和观测资料的可比性，对不同震相到时先进行了初步估计，在本文所考虑的时长范围内共需计算27个震相，如图 10所示.这些震相主要包括Moho界面的一次透射转换波和第一、二层之间多次的透射和反射波.

图 10

Fig. 10

图 10 真实数据反演中所使用的震相路径示意图 图中实线代表P波，虚线代表SV和SH波，共计27个震相. 包括Moho界面的一次透射转换波和第一、二层之间多次的透射和反射波. Fig. 10 Schematic diagram of wave phase used in inversion The solid line is P wave phase, the dotted line is SV and SH phase. There are totally 27 phases in used, include first transmission phase from Moho, transmission and refraction phase between first and second layer.

考虑到分层的合理性，在初次反演的基础上，进行了二次反演.根据正演数值试验结果，当plunge接近90度时，所得到的图案与无各向异性时是一致的.因此当一次反演结果中plunge接近90度时，认为该层是各向同性的.另外，还可以根据一次反演的速度分布合并相应的层.这样在二次反演中反演的参数个数就会有所减少，此时可以加入各向异性程度的反演.

各台站反演结果如图 11~13所示，各层各向异性的快轴方向也显示于图 9.可以看出，所有台站反演结果中，在第三层存在低速(S波)层，且波速比(k)较高.该结果较符合地壳流的物性特征.根据前人的研究，一般认为地壳中的各向异性主要存在于上地壳和下地壳，而本研究也在下地壳发现了各向异性层，且其各向异性快轴方向与前人推测的地壳流方向一致，该结果从地震学上得到了地壳流的直观地震学证据.本研究在上地壳中未发现明显的各向异性特征，而在该地区中地壳却普遍发现各向异性，而且各向异性快轴方向与断层方向近似垂直，并与印度板块挤压青藏高原的方向一致，初步推断，中地壳的各向异性可能来自于块体的相互推挤作用导致的应力积累.

图 11

Fig. 11

图 11 合作台反演结果 (a)为观测接收函数，左小图为R分量，右小图为T分量；(b)为反演最终的模型对应的接收函数，左小图为R分量，右小图为T分量；(c)为反 演模型部分参数的搜索结果，vs为每层S波平均速度，k为每层波速比，pctS为每层S波各向异性程度，trend为每层各向异性快轴方向在水平面投影和正北的夹角.(d)为反演搜索过程中，misfit随迭代次数的变化，蓝线代表每次迭代中misfit的平均值，红线代表每次迭代中misfit的最小值. Fig. 11 Hezuo inversion result (a)Observed receiver function, the left is R component, the right is T component.(b)Receiver function for final model, the left is R component, the right is T component.(c)Inversion result for parts of model parameters. vs is average S velocity, k is velocity ratio, pctS is S wave degree of anisotropy, trend is angle between the projection of fast axis of anisotropy on the surface and north orientation.(d)Misfit value changes with iteration, the blue line st and s for the average misfit values, the red line st and s for minimal misfit values.

图 12

Fig. 12

图 12 迭部台反演结果 (a)为观测接收函数，左小图为R分量，右小图为T分量；(b)为反演最终的模型对应的接收函数，左小图为R分量，右小图为T分量；(c)为反 演模型部分参数的搜索结果，vs为每层S波平均速度，k为每层波速比，pctS为每层S波各向异性程度，trend为每层各向异性快轴方向在水平面投影和正北的夹角.(d)为反演搜索过程中，misfit随迭代次数的变化，蓝线代表每次迭代中misfit的平均值，红线代表每次迭代中misfit的最小值. Fig. 12 Diebu inversion result (a)Observed receiver function, the left is R component, the right is T component;(b)Receiver function for final model, the left is R component, the right is T component;(c)Inversion result for parts of model parameters. vs is average S velocity, k is velocity ratio, pctS is S wave degree of anisotropy, trend is angle between the projection of fast axis of anisotropy on the surface and north orientation;(d)Misfit value changes with iteration, the blue line st and s for the average misfit values, the red line st and s for minimal misfit values.

图 13

Fig. 13

图 13 岷县台反演结果 (a)为观测接收函数，左小图为R分量，右小图为T分量；(b)为反演最终的模型对应的接收函数，左小图为R分量，右小图为T分量；(c)为反演模型部分参数的搜索结果，vs为每层S波平均速度，k为每层波速比，pctS为每层S波各向异性程度，trend为每层各向异性快轴方向在水平面投影和正北的夹角.(d)为反演搜索过程中，misfit随迭代次数的变化，蓝线代表每次迭代中misfit的平均值，红线代表每次迭代中misfit的最小值. Fig. 13 Minxian inversion result (a)Observed receiver function, the left is R component, the right is T component;(b)Receiver function for final model, the left is R component, the right is T component;(c)Inversion result for parts of model parameters. vs is average S velocity, k is velocity ratio, pctS is S wave degree of anisotropy, trend is angle between the projection of fast axis of anisotropy on the surface and north orientation;(d)Misfit value changes with iteration, the blue line st and s for the average misfit values, the red line st and s for minimal misfit values.

4.1  本研究通过系统的数值试验，分析了各向异性和倾斜结构对接收函数的影响，总结了不同参数对接收函数的影响规律.结合邻域算法，利用接收函数反演了地壳的各向异性结构参数.根据数值试验结果，确定了各向异性结构分步反演的方法.并将该反演方法应用到了青藏高原东北缘固定台站下方地壳各向异性结构的反演中.结果显示，在青藏高原东北缘三个台站下方的下地壳存在明显的各向异性低速结构，该结构和下地壳流模型较为一致.

4.2  此外，数值试验和实际资料处理结果表明，在六方晶形对称结构的地壳各向异性模型中，诸因素是独立影响接收函数图案的.其中各向异性和介质倾斜的混叠很难用常规方法区分开，但根据我们的反演能够稳定的给出各向异性轴的方向.加上一些先验信息(层厚、波速比)的给定，以及二次反演等策略，最终可以将各向异性参数反演到可接受范围之内.

4.3  目前，横波分裂法是各向异性分析的主要方法，而本方法的优势在于能够确定各向异性的深度分布特征.在青藏高原东北缘，张辉等(2012)运用横波分裂法给出了该地区地壳的平均各向异性分布，王琼等(2013)给出了该地区上地幔的各向异性分布，经对比发现本文中的第三层各向异性与王琼等的较为接近，可以认为是受到了地幔物质的影响，但第二层与第三层明显解耦似乎受断层影响较大，由于计算深度不同，与张辉等的地壳各向异性存在差异.在未来的工作中我们将尝试引入横波分裂法的信息，联合接收函数对反演做约束，使结果更加真实可靠.

致 谢 本文采用了Rickwood和Sambridge(2006)的邻域算法代码和Frederiksen(2003)的正演代码，图件绘制中使用了matplotlib和GMT(Wessel and Smith, 1995)，地震波形观测数据由甘肃省地震台网中心提供，对此表示感谢.另外，感谢甘肃省地震局监测中心为本研究提供了计算平台.