区域古构造应力场研究
1.古构造应力的测定
构造活动的强度反映了构造应力的大小，通过类比可以近似估算出古应力的大小。声发射法是一种较精确的计算古应力的方法(丁，2000，2001)。在边界条件已知的情况下，数值模拟方法是确定古应力场最有效的方法，不仅可以确定盆地内详细的古应力值，还可以确定古应力方向。为了进一步认识济阳坳陷北部带馆陶组的构造特征和形成机制，有必要扩大区域范围，以济阳坳陷为对象进行研究。
(1)古应力方向的确定
①区域构造和动力机制反映的应力场方向。深部地幔物质的上涌控制着浅部地壳的伸展和减薄，产生张应力场。在地幔物质不断上升的过程中，如果上升流强度过大，地壳隆升大于地壳伸展减薄，在区域构造挤压环境发生时会遭受短期抬升剥蚀；之后，如果以地幔上涌导致的地壳变薄和盆地裂陷伸展为主，盆地在短暂隆升后会发生张性裂陷。这种动态不平衡可能在地幔上涌过程中多次形成，导致盆地裂陷过程中多次短期构造抬升，直至地壳下岩石圈达到热能积累-释放平衡的极限，岩石圈热松弛、冷却、增厚，盆地进入整体坳陷阶段。渤海湾盆地新生代盆地演化具有这一特征。由于地壳减薄与上地幔温度上升之间的地壳平衡关系处于持续的不平衡状态，东营运动正处于由伸展裂陷向坳陷的过渡期。这种持续的不平衡状态的积累，使盆地的地壳减薄和伸展裂陷运动在新近纪末达到极限，也达到了地壳下岩石圈热能积累-释放平衡的极限，使盆地在短暂隆升后开始收缩，使岩石圈增厚。
渤海湾盆地在坳陷早期(馆陶组沉积早期)受热平衡控制，沿古近系北北西向仍有微弱的残余伸展减薄。由于伸展的继续，盆地边界的断裂带有强烈的火山活动。馆陶组沉积时期形成的火山岩主要为玄武岩(曾广策，1997)，尚林断裂带、高青断裂带、八面河断裂带均有火山活动，以平静溢流为主(宗，1999)。玄武岩被认为是拉伸环境的产物。曾广策等(1997)估算沉积期馆陶组的拉张速度约为0.30cm a- 1+0，明显低于沙三段强烈裂陷期的0.45cm a- 1+0。
图2-20渤海湾新近纪构造背景分析(根据侯桂亭等，2001，修改)
在新近纪(24Ma以来)，由于印度板块向欧亚板块的汇聚运动，中国西部向东逃逸(侯桂亭，2001)，向东分散，造成不同地块之间的滑动应力场。同时，太平洋板块退到欧亚俯冲带，倾角增大，日本海盆地开始扩张(乔汉生，1999)，推挤渤海湾盆地。这样在渤海湾盆地边界形成走滑剪压构造应力场，北北东向边界承受右旋走滑剪压应力，NWW对边界起调节平衡作用，主要是左旋走滑剪切(图2-20)。
②构造特征和断层活动反映的应力场方向。馆陶组构造主要继承了古近系构造格架，断裂活动前期强烈，后期减弱。早期以NEE和近东西向断层活动为主(图2-21)，断层落差较大，一般在50 m以上，而北西向断层基本不活动或仅轻微活动，落差通常小于20 m，表明北北西向地区的张性构造应力场约为NNW340，应与NEE向断层近垂直，对该类断层有利。NWW走向的断裂也有一定的活动性，如陈南断裂、紫珍断裂南部的断裂、东部与仁丰断裂相邻的NWW走向的断裂。但这些断层的落差远小于NEE向断层的落差，说明区域拉应力相对较弱，拉应力场方向与断层走向的夹角是断层活动强度的控制因素。
图2-21济阳坳陷馆陶组沉积早期断裂系统分布图
与馆陶组早期相比，馆陶组晚期断裂活动性明显减弱(图2-22)，北西向断裂活动性明显增强。如成南断裂东部早期断落差只有30 m，晚期断落差90 m，表明区域应力场偏转，最小主应力可能从早期的NNW 340变为晚期的NNE 10°，岩石圈。
明化镇组沉积末期，最小主应力方向变为近南北向，挤压作用进一步增强。东营凹陷和惠民凹陷的断裂活动基本停止，沾化凹陷东部和埕北地区形成了许多东西向和NEE附近的小断裂。应力场逐渐向现今应力场方向变化。
(2)古应力的测定
虽然没有可靠的方法来确定古应力的大小，但通过类比可以大致估算出应力变化的范围。之前分析过，在馆陶组沉积早期，由于岩石圈热平衡，济阳坳陷可能存在微弱的拉应力，最大拉应力分量不应超过4 MPa(馆陶组岩石抗拉强度约为2 MPa)，平均拉应力应在1 MPa以内，否则会产生较大的张张量，导致强烈的断裂活动。印藏碰撞和太平洋板块俯冲的长程作用对该区挤压作用相对较小，一方面应力传递过程要衰减和消耗；另一方面，馆陶组是近地表堆积，构造环境相对稳定。盆地边界的压应力主要起约束作用，盆地内不可能有太大的压应力，最大主应力应小于同一深度的垂直主应力。凸边由于地形高差的影响，应该是压应力相对集中的区域。同理，馆陶组沉积后期和明化镇组沉积时期的区域压应力也不会太大，声发射测试结果可以反映出来。样品取自新2井，目前深度为2065.5米，沙河街组二段岩心声发射测试显示，馆陶组末最大主应力为2.2 MPa，明化镇组末最大主应力为3.6 MPa(徐建春，2004)。显然，馆陶组的压应力比沙河街组弱。
图2-22济阳坳陷晚期断裂系统示意图
(3)通过数值模拟得到古应力场。
①数值模拟的基本原理。应力场的数值模拟通常采用有限元法，其基本思想是将所研究的连续体简化为由有限元组成的离散模型，然后用计算机进行数值求解(徐建春，2004)。将一个地质体离散成有限个连续的单元，这些单元通过节点连接，每个单元被赋予实际的岩石力学参数。求解研究区域内的连续场函数转化为求解有限离散点(节点)处的场函数值，基本变量为位移、应变和应力。根据边界应力条件和节点平衡条件，建立并求解以节点位移为未知量、以整体刚度矩阵为系数的方程组，通过插值函数得到每个节点上的位移，进而计算出每个单元中的应力和应变值。随着单元数的增加，越接近实际地质体，解就越真实和精确。有限元法的基本步骤是(刘泽荣，1983):
A.结构或对象的离散化；
b .在单元格中选择场变量插值函数；
c .进行单位计算，得到单位特征矩阵和数组；
d .进行全局分析，组装全局矩阵和数组，建立全局方程；
E.计算单元内部的场变量。
②模拟的简化条件。济阳坳陷内部地质构造条件极其复杂，边界断层凹凸不规则。在应力场模拟中必须做一些必要的简化，以便于了解应力作用方式和边界条件。
A.济阳坳陷总面积(包括陆地、海域和外围)约2.70 × 104km2，其长度(约260 km)和宽度(约194 km)远大于其厚度(馆陶组厚度一般为300 ~ 900m)，可将三维应力应变问题简化为平面应力问题。
图2-23济阳坳陷及其外围主要边界(应力场模拟目标区)
B.以济阳坳陷为模拟目标区，其东边界由郯庐断裂西支组成，南边界为鲁西隆起，西边界与北向主要边界断裂(如紫珍断裂)末端切割，北边界为成宁隆起(图2-23)。这种边界极不规则，可以扩展为规则形状，主要地质边界可以作为远场应力边界，便于应力施加和边界条件约束。这样实际应力场模拟的边界就变成了:成宁隆起北部北北西走向的断裂为北边界，郯庐断裂为东边界，鲁西隆起为南边界，西边界为兰聊断裂北北东向的走滑断裂带(图2-24)。
C.由于难以知道新近纪实际应力场的大小，外力在模拟中通过类比施加一个虚拟值。只要虚拟外力符合客观地质规律，模拟应力场和实际应力场就可以足够接近，而应力场的变化趋势是完全一致的。
③应力场模拟检验标准。数值模拟的结果需要一定的测试标准来衡量和评价，从而判断模拟的有效性和准确性。但古应力场离我们太远，目前还没有有效的方法对古应力场进行评价，对应力场的模拟结果也没有成熟可靠的检验标准。本书认为，一个成功的古应力场模拟至少应该达到以下两点:
A.模拟的应力场应与模拟区域的构造强度有很好的对应，模拟结果能从应力场的角度解释构造活动的性质和强度。
b .模拟应力值应符合客观地质规律，应充分考虑应力值与岩石强度、埋深、构造特征等因素的关系，确保应力场模拟结果的有效性和准确性。如果可能，应参考具有相似构造背景的实测应力值来模拟应力场。
图2-24延伸到主构造边界后应力场模拟的实际范围。
2.馆陶组沉积早期古应力场的数值模拟。
(1)地质模型
地质模型的建立是应力场模拟的第一步，也是最关键的一步。好的地质模型不仅容易计算边界的远场应力和约束，还能反映内部的构造特征。建模时以济阳坳陷为应力场模拟的目标区域，北边界延伸至成宁隆起以北的北西向断裂带，西边界沿兰聊断裂北北东向延伸至走滑断裂带，从而得到应力场模拟的规则边界区域，便于应力施加和位移约束。
以济阳坳陷馆陶组沉积早期的断层分布图(图2-9)为基础，将主要的凸起和凹陷边界投影到这张底图上，数字化提取断层、凸起和凹陷边界的坐标位置，然后输入Ansys软件建立实体模型(图2-25、图2-26)。
图2-25馆陶组沉积早期应力场模拟的立体地质模型图中，后面带数字的字母表示区域实体的数量。
图2-26在馆陶组沉积早期应力场模拟的凹陷区域立体地质模型中，字母后接数字表示区域实体的数量。
(2)选取力学参数和单元划分
不同构造单元在平面上力学性质不同，一般断裂带正常沉积地层强度较弱，凸起区老地层强度强于凹陷区新地层(刘泽荣，1983；陈波，1998)，大型断裂带的岩石强度比小型断裂带的岩石强度弱。根据实际模拟经验，这种力学性质的差异对应力场的分布影响不大，但对应力的大小影响显著。根据实际情况，对应力模拟的不同区块赋予不同的力学参数，简化抽象出四种类型:馆陶组沉积的疏松砂岩和泥岩、隆起(或隆起)区前古近系相对致密坚硬的岩石、弱化的一级断裂带和二级断裂带。其中，济阳坳陷一、二级断裂带被赋予一种类型的力学参数，三、低级断裂带被赋予另一种力学参数。力学参数参考一些实验测试(宋树军，2003；游明清，2003)，如表2-4所示。
表2-4沉积期馆陶组应力场模拟力学参数表
确定力学参数后，可以通过网格划分形成有限元模型。选取具有4个节点的平面四边形单元，在Ansys的图形用户界面(GUI)下改变力学参数选项。依次选取力学参数一致的区域面积，对实体模型进行网格化，共划分出89271个四边形单元，其中济阳坳陷87288个单元。
(3)确定加力方式和边界条件。
由于新近纪印藏碰撞效应，中国大陆向东滑动，太平洋板块也向欧亚大陆移动，限制了中国大陆向东移动，使中国大陆整体处于压应力状态。郯庐断裂带和兰聊断裂带表现为弱右旋走滑性质，济阳坳陷存在弱北北西向拉应力。因此，确定馆陶早期外力为西、东边界的压应力和右旋剪应力，南、北边界的左旋剪应力。在南部和北部边界的四个顶点上施加位移约束，以避免刚体运动。此外，NNW340向济阳坳陷及边界断层施加弱张应力，自西向东为自镇-阳新断层、埕南断层和渤南隆起南边界断层，走向近NE和NEE(图2-27)。经过反复试验和计算，最终从57个方案中确定边界应力的大小:西边界施加0.15 MPa右旋剪应力和0.35MPa压应力；东边界(郯庐断裂带)右旋剪应力0.08 MPa，压应力0.10 MPa；南边界(鲁西隆起)施加0.10 MPa的左旋剪应力；北边界(成宁隆起北西向断裂)施加0.10 MPa的左旋剪应力；向紫镇-阳新断裂、城南断裂和渤南隆起南边界断裂施加NNW 340±0.15 MPa的张应力(附图2-27)。外力总体上表现为弱压缩、弱拉伸和弱走滑，反映了坳陷早期的近地表应力场特征。
(4)模拟计算和后处理
利用Ansys软件，通过弹性平面有限元程序对每个模拟方案进行计算，调用后处理模块对计算结果进行绘制和处理，生成反映应力场特征的等值线图和方位矢量图(图2-28和图2-29)，进行分析比较，确定下一个方案的实施，直至出现一个方案。本方案实施前后，无论如何修改模拟方案，如果模拟结果比本方案差，则停止模拟，输出本最佳方案的模拟结果，如图2-30应力场模拟流程图所示。
图2-27馆陶组沉积早期应力场模拟及其边界条件
图2-28济阳坳陷馆陶组水平最小主应力等值线的正值为张应力。
图2-29济阳坳陷馆陶水平最大主应力等值线的负值为压应力。
(5)计算结果分析
馆陶组沉积早期的最小主应力全部为拉应力(图2-28)，在济阳坳陷呈现“东西分带、高低交替”的拉应力分布格局，自西向东交替出现两个高应力区和两个低应力区。第一张应力集中区位于惠民凹陷西北部；第二张应力集中区大致沿仁丰断裂和石村断裂穿过彬县隆起，到达沾化凹陷西部，然后向北延伸至车镇凹陷。应力值可达0.1 ~3.7 MPa。紫珍断裂与伊东断裂之间的西南向带状区是第一个低张应力区，应力值小于500Pa。石村断裂东部，绕过东营凹陷中央隆起带，陈家庄凸起东部，沾化凹陷东部和埕北广大地区构成第二低张应力区，应力值为0.01 ~0.05 MPa。
最大主应力全部为压应力(图2-29)，与拉应力正好相反，在济阳坳陷呈现“南北分带，依次增大”的应力分布格局。一般沿东北方向，自西向东，从惠民凹陷到车镇凹陷，再到城北地区，压应力高，应力值为0.3 ~10.2 MPa，陈家庄凸起也在该区；在本区南部，东营凹陷中央隆起以北，与沾化凹陷构成第二高的压应力区，应力值为0.2 ~0.34 MPa。中央隆起以南的东营凹陷和董卿凹陷是压应力最小的地区，应力值一般小于0.2 MPa。
整体水平剪应力表现为逆时针左旋剪切(图2-31)，右旋剪切只发生在局部区域。
图2-30古应力场模拟流程图
图2-31济阳坳陷馆陶组沉积早期水平剪应力等值线图(逆时针负向)
最大主应力(压应力)方向在惠民凹陷西部为北西向，中部变为近南北向，东部在东营凹陷、沾化凹陷和车镇凹陷逐渐偏转至近东西向，在埕北地区进一步偏转至NEE(图2-32)。最小主应力(拉应力)在惠民凹陷西部开始由北北东向偏移，在惠民凹陷中部向近东西向偏移，东部逐渐过渡到近南北向和北北西向(图2-32)。将馆陶组沉积早期的应力场模拟结果与断层落差图对应起来，可以很好地解释馆陶组沉积早期的断层活动特征。图2-29:高张应力区对应高断层落差区，说明应力场的分布控制着断层的活动性。应力场对断裂活动的控制不仅体现在断层落差与拉应力的一一对应关系上，还与应力场方向的变化和差异应力有关。在惠民凹陷西部，由于拉应力和压应力较大(图2-29，图2-30)，形成了较大的差应力，达0.6 ~ 4.0 MPa。同时，最小主应力的NNE向与这些近NE向和NEE向的断层夹角较大，有利于断层活动。紫珍断裂和尚林断裂西段落差较大，均大于70m；惠民凹陷中部地区为低张应力区，张应力方向变得接近东西向，与断层走向近平行，故断层活动较弱，如尚林断层东部的断层，最大落差不超过20m；仁丰断层、石村断层、南林断层、贝林断层以东，拉应力方向变为近南北向和北北西向，与主断层近垂直相交，拉应力值相对较高，断层活动性较强，特别是车镇凹陷，断层落差普遍较大；石村断层以东、陈家庄凸起以东、沾化凹陷以东和城北地区，虽然号称应力方向和断层走向也是近垂直的，但张应力值普遍较低，差异应力值较小，断层活动性较弱。
图2-32济阳坳陷馆陶组沉积早期水平最小主应力方位图
图2-33济阳坳陷馆陶组沉积早期水平最大主应力方位图
3.馆陶晚期古应力场的数值模拟。
馆陶组晚期的断层活动性比早期弱，许多断层停止活动。显然，沉积早期的馆陶组地质模型不能用于后期的应力场分析，需要建立新的模型。
(1)地质模型
建模时仍以济阳坳陷为模拟目标区，向外延伸，东边界至郯庐断裂，南边界至鲁西隆起，西边界沿兰聊断裂成宁隆起西部的北北东向走滑断层，北边界至成宁隆起北部的北西向走滑断裂带，作为远场应力作用边界和约束边界。由于馆陶组沉积后期和沉积初期的地质模型基本相同，但内部断层不同，因此建模时只需修改部分实体，删除图2-25所示实体模型中不必要的断层带，最终生成新的地质模型，如图2-34所示。
(2)选取力学参数和单元划分
考虑到晚期与早期差异不大，沉积物主要堆积在地表附近，因此应力场模拟的力学参数与早期一致。后期断层数量大大减少，图2-36中的实体模型需要重新网格化，生成有限元模型。选择平面4节点4边单元，改变Ansys图形用户界面(GUI)下的力学参数选项，依次选择力学参数一致的区域区域(图2-37)，将实体模型网格化，划分为145359个四边形单元，其中济阳坳陷143356个单元，如图2-38所示。为了便于施加不同方向的边界外力，在边界区域的节点上增加了437个“表面效应单元”，这样实际参与运算的单元就有145819个。
(3)确定加力方式和边界条件。
馆陶后期试验了20种施力方案，最满意的方案是:西边界施加0.25 MPa右旋剪应力和1.0 MPa压应力；东部边界(郯庐断裂带)施加0.25 MPa的右旋剪应力和0.65 MPa的压应力；南边界(鲁西隆起)左旋剪应力0.15 MPa，压应力1.0 MPa；在北部边界(成宁隆起北部NW向断层)施加0.15 MPa的左旋剪应力，位移约束在模型的四个顶点(图2-35)。与馆陶早期挤压走滑相比，外力加强了，区域拉应力消失了。从应力场模拟结果看，拉应力主要出现在局部，是应力平衡的结果，反映了早期坳陷向全面坳陷期逐渐过渡的近地表应力场特征。
图2-34沉积晚期和明化镇期馆陶组应力场模拟立体地质模型
(4)模拟计算和后处理
模拟计算和后处理过程与馆陶早期相同，如图2-30所示。与早期应力场模拟结果不同的是，垂直主应力始终为零，两个水平主应力可能在零附近变化，这样垂直主应力可能是三维空间中的最大、最小或中间主应力，因为An-sys软件是按照大小输出三维空间中的应力值，而只有平面中的最大和最小主应力才是期望的结果。因此，需要调用Ansys软件的APDL (Ansys参数化设计语言)编写相应的后处理模块，才能输出有效的仿真结果。早期应力场模拟的结果表明，平面内的最大和最小主应力正好是三维空间的最大和最小主应力，所以这一步是不必要的。
(5)计算结果分析
济阳坳陷在馆陶组沉积晚期进入坳陷期，应力场以挤压为主。最小主应力一般由馆陶组沉积早期的拉应力变为压应力(图2-36中的压应力区)，仅在惠民凹陷与东营凹陷的过渡带、车镇凹陷东部、沾化凹陷东北部和埕北地区有微弱的拉应力，拉应力值小于0.1MPa。拉应力主要受盆地边界条件制约，是应力场调整平衡的结果，反映了局部特征。
图2-35沉积晚期馆陶组应力场及边界条件模拟
图2-36济阳坳陷馆陶组沉积晚期水平最小主应力(压应力为负)等值线图。
沉积晚期馆陶组最大主应力分为东西向和南北向两个区块。惠民凹陷东部自西向东构成第一压应力低值区，应力值为0.5 ~1 MPa。东营凹陷、沾化凹陷和埕北地区构成第二个压应力低值区，中部陈家庄凸起为压应力高值区，应力值约为1.57 ~2.36 MPa，将其分为两个区块(图2-37)。
图2-37沉积晚期济阳坳陷馆陶组水平最大主应力等值线图(压应力为负)
平面剪应力在馆陶组沉积晚期主要表现为顺时针右旋走滑，在盆地中部尤为突出，剪应力在0.3 MPa 0.3 ~2.3 MPa之间。在盆地的东西两侧，平衡边界向右走滑，表现为较弱的左旋剪应力场，左旋剪应力值相对较小，一般小于0.71 MPa(图2-38)。
最大主应力方向以阳新凹陷-林范家-彬县-博兴线为界，形成两个“漩涡”。应力场西部从惠民凹陷西南角的NNW转向惠民凹陷北部的NNE，再转向惠民凹陷东南部的SEE，形成一个顺时针方向的应力场“漩涡”。东部地区最大主应力方向逐渐由NE向NEE过渡，形成一个逆时针方向的“漩涡”(图2-39)。最小主应力方向表现出由西南向东北收敛的特征。从西南向东北，最小主应力方向从近东西向逐渐过渡到近南北向，从东南向东北，最小主应力方向从北北西向逐渐过渡到近南北向(图2-40)。最小主应力的低值区正好对应断层活动强度高的区域，这就决定了惠民凹陷东部馆陶组断层活动虽然总的趋势是早期强于晚期，但晚期比早期强。沾化凹陷北部和车镇凹陷东部以及埕北地区存在弱压应力至弱拉应力，最大主应力方向与断层走向几乎一致，有利于断层活动的发生，主要断层落差较大。
图2-38沉积晚期济阳坳陷馆陶组水平剪应力等值线图(逆时针负向)
图2-39沉积晚期济阳坳陷馆陶组水平最小主应力方向图
图2-40沉积晚期济阳坳陷馆陶组水平最大主应力方位图