水量均衡法

（一）基本原理

水量均衡法是根据水量平衡原理，建立均衡方程计算水量的方法，表达式为

∑Q补-∑Q排=ΔQ储 （3-1）

式中：∑Q补为均衡期内地下水系统各种补给量的总和（m3）；∑Q排为均衡期内地下水系统各种排泄量的总和（m3）；ΔQ储为均衡期内地下水系统内部储存资源的变化量（m3）。

（二）一般步骤

1.确定均衡区

根据地下水系统理论的要求，均衡区应是地下水系统边界所界定的空间范围，一般要求以地下水系统天然边界作为划分依据。由于水量均衡法属于集中参数系统，为了提高区域地下水数量评价精度，在实际计算时可以根据不同水文地质条件划分为不同级别的子区，分别计算各均衡要素，然后进行综合。例如根据给水度、降水入渗系数、地下水埋藏深度等条件，将均衡区划分为若干子区，分别计算各子区的储变量、降水入渗量和潜水蒸发蒸腾量，然后求和。

2.确定均衡要素

确定式（3-1）中∑Q补和∑Q排的组成，即确定地下水系统三维空间区域边界上的输入和输出量。从外界进入地下水系统的各种水量统称为补给项，系统输出的各种水量统称为排泄项。

一般而言，补给项包括：大气降水入渗补给量、地表水体渗漏补给量（河流、湖泊、水库等）、地下侧向流入补给量、越流补给量、凝结水补给量、地表水灌溉入渗补给量、地下水灌溉回归补给量、渠系渗漏补给量、人工回灌补给量等。

排泄项包括：潜水蒸发蒸腾量、地下水侧向流出量、地下水开采量、泉水溢出量、越流排泄量、向河湖排泄量等。

需要指出的是，不同的地下水系统与外部环境之间的水量交换关系不同，所以均衡要素的组成因不同地下水系统而异。在实际工作中，需要与研究区具体条件紧密结合，确定均衡要素的组成。

3.确定均衡期

地下水均衡计算是针对某一特定时间段进行的，称为均衡期。如前所述，在地下水的资源功能评价中，要求地下水数量评价的时间尺度为5～12年，以此为均衡期进行水量均衡计算。为保证水量平衡，各均衡要素计算和相关的资料的选取应采用统一的时间序列。

（三）均衡项计算方法

1.降水入渗补给量

降水入渗补给量确定方法包括：直接测定法、零通量面法、包气带达西定律法、氯质量平衡法、示踪法等。

（1）直接测定法

通常利用测渗仪或通过包气带蒸渗试验直接测定不同岩性、不同地表覆盖情况下的降水入渗补给量（Young et al.，1996）。我国于20世纪70年代末期开始在华北地区和西北地区建立了许多包气带试验场，开展了大量的实验研究。

（2）零通量面法（ZFP）

零通量面是Richards于1956年首先提出的，是包气带水分运移的分界面，其上土壤水分向地表运移，其下水分向地下水运移，将该面以下的水分运移速率作为地下水补给速率，利用该法需要测定包气带垂直剖面土壤水势和含水量。我国于20世纪80年代后期引入ZFP法（张光辉，1988；张惠昌，1988），目前该方法仍在应用（程辉等，2000；周金龙等，2003；李茜等，2006）。

（3）包气带达西定律法

达西定律法是干旱、半干旱地区常用的方法，需要测定包气带水力梯度和不同含水量下的渗透系数，计算公式如下：

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/d）；K（θ）为包气带水渗透系数（m/d）；H为包气带水侧压水头（m）；Z为垂直位置高程（m）。

（4）氯质量平衡法

该方法主要应用了氯的化学稳定性，其应用前提是（Kinzelbach，2002）：①由于包气带溶质输入和向饱水带的输出存在时间滞后，所以必须假定在此期间没有重要的气候变化；②没有额外的溶质加入，如肥料，同时也没有近期大气污染；③在ZFP之上和之下没有溶质储存的净变化，这种变化可能由于动植物引起或矿物沉淀/溶解和吸附/解吸附。在满足以上条件的基础上，可采用下式计算降水入渗补给速率：

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/a）；P为多年平均降水量（m/a）；CP为降水中氯离子浓度（mg/L）；Fd为氯离子干沉降量（mg/m2·d）；CS为零通量面以下包气带水的氯离子浓度（mg/L）。

（5）示踪法

利用人工或环境示踪剂，通过失踪剂峰面移动来计算降水入渗补给速率。常用的人工示踪剂包括：氚、溴、碘、染色剂等（Athavale，1988；Kung，1990；Flury，1994；Aeby，1998；Forrer，1999），环境示踪剂包括氚、氯-36 等受核爆影响的放射性同位素（Scanlon，2002）。降水入渗速率计算公式为

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/a）；Δz为示踪剂浓度峰面的运移深度（m）；Δt为峰面运移时间（a）；θ为体积含水率（无量纲）。

由于受地形地貌、地表覆盖、包气带岩性及厚度、降水强度及频率、包气带水分状况、地下水埋深等条件的影响，不同地带的降水入渗速率不同。而采用以上方法获得的数据仅是某一或某些条件下的实验结果，所代表的空间尺度有限，且不同的方法所代表的时间尺度也不相同（表3-2）。因此，在区域地下水资源评价中，往往根据研究区实际条件进行适当分区，选用不同方法求得不同分区的降水入渗系数（a），然后采用下式计算降水入渗补给量：

Qp=a·P·F （3-5）

式中：Qp为降水入渗补给量（m3/a）；a为降水入渗系数（无量纲）；P为多年平均降水量（m/a）；F为计算区面积（m2）。

表3-2 不同方法确定的降水入渗速率范围及时空尺度对比

注：表中数据根据Scanlon等（2002）整理。

目前，我国北方大部分地区已经通过包气带入渗试验、水位动态分析等方法，建立起不同地区降水入渗系数与地下水位埋深和包气带岩性之间的关系。

2.地下水与河流之间交换量

（1）断面流量差法

若均衡区有河流穿过，则在均衡区上、下游边界处各选一个测流断面监测流量，并确定断面之间的距离、测流时间间隔、河流水面宽度和水面蒸发量，然后采用以下公式计算：

Qr=（Q1-Q2）·Δt-B·L·E （3-6）

式中：Qr为测流期间河道渗漏补给量（m3）；Q1，Q2分别为河流上、下游断面的平均流量（m3/s）；Δt为计算时段（s）；B为河流水面平均宽度（m）；L为河流两断面间的距离（m）；E为测流期间的水面蒸发量（m）。

（2）渗流断面法（达西定律）

当河水与地下水有直接水力联系时，采用达西定律计算河道侧渗量，公式为

Qr=K·L·I·h·Δt （3-7）

式中：K为含水层渗透系数（m/d）；L为河道渗漏段长度（m）；I为河渠一侧地下水水力梯度（无量纲）；h为过水断面的厚度（m）；Δt为计算时段（d）。

h的取值应根据河流与地下水的关系而定。当河流一侧接受地下水补给，另一侧补给地下水时（图3-2a），h取值为河床到地下水位（河水位）的距离；当河流两侧都补给地下水时（图3-2b），h取值为含水层的整个厚度。

（3）基流分割法

在地下水补给常年性河流的地区，在枯水期河水流量几乎全部由地下水补给维持，这时的河水流量被称为基流量。把河流流量过程线上的基流量分割出来，即为地下水对河流的补给量（房佩贤等，1987；曲焕林等，1991；徐恒力等，2001）。图3-3为典型的单峰流量过程线，由起涨部分、峰值和退落部分组成。起涨部分的起点称为起涨点（图3-3中的a点），在退落部分，当降水影响消失时，河流量由地下径流组成，其起点称为地下水退水点（图3-3中的d点）。起涨点很好确定，而确定地下水退水点比较困难，一般有3种方法：经验法、退水曲线法和作图法。

图3-2 河流与地下水补排关系示意图

经验法，就是在过程线上的退落部分，找到曲线曲率最大的点即视为地下水退水点。

退水曲线法，认为从退水点开始，流量变化满足布西涅斯克方程（退水曲线方程）：

Q=Q0·e-kt （3-8）

式中：Q为从d点开始的任一时刻的河流量（m3/s）；Q0为d点的流量（m3/s）；k为衰减系数；t为以d点为起点的时间（s）。

由式（3-8）可知，从退水点开始流量呈等比级数递减。利用这个规律，在退落部分找到流量大体成等比递减的开始时刻，即为地下水退水点。

作图法，在退落部分按相等的时段，选取一系列流量，计算流量差（ΔQ），然后以ΔQ为纵轴，以时间（t）为横轴绘制曲线，将曲线的拐点所对应的时刻作为退水曲线的初始时刻，然后在图3-3中找到该时刻所对应的点即为退水点。

图3-3 河流流量过程线

由于河流与地下水之间的水动力关系不同，基流分割方法也不同。一般有两种情况（图3-4）：一种情况是河流与地下水有直接水力联系；另一种情况是二者之间不存在直接水力联系。

当河流与地下水无直接联系时（图3-4（a）），若不考虑地下径流峰值，则直接连接ad，其下方阴影部分的面积即为基流量（图3-5（a））；若考虑地下径流峰值，则可分别计算流量过程线起涨部分（af段）的平均流量（Q1）和退落部分（df段）的平均流量（Q2）。然后，再计算出起涨时段内大气降水形成的平均流量（Q′1）：

图3-4 河流与地下水关系示意图

区域地下水功能可持续性评价理论与方法研究

式中：P为af段的总降水量（m）；F为测站所控制的流域面积（m2）；taf为从a点到f点经历的时间（d）；v为径流系数，等于径流深度与降水深度之比。

同理，求出退落时段内大气降水形成的平均流量（Q′2）之后，在流量过程线上找到分别与 和 相对应的点，过两点作垂直t轴的直线，在两直线上分别减去Q′1和Q′2，可得到b，c两点（图3-5（b）），连接a，b，c，d，其下的阴影部分面积即为基流量。

图3-5 地下水与河流无直接水力联系时的基流分割图示

当河流与地下水有直接水力联系时（图3-4（b）），枯水季节地下水补给地表水，丰水季节河水位高于潜水水位，地表水补给地下水。此种情况下，在枯水季节（起涨点之前和退水点之后），河流流量全部为基流量。将起涨点对应的时间记为ta，退水点对应的时间记为td。进入洪水期后，河水开始补给地下水，但在ta之前进入河道的地下水与洪水一起从上游向下游流动。河流源头到测站的距离很容易测定，则可以求出河水从源头流到测站所用的时间（记为Δt），也就是说，到ta+Δt时刻河水全部由洪水组成，其对应的点记为b点，则地下水径流量应按 ab 逐渐减少；同样，在洪峰过后，河流源头首先有地下水进入河道，起始时刻为td-Δt，其对应的点记为c，则地下水径流量应按 cd 逐渐减少。基流分割如图3-6中的阴影部分。

图3-6 地下水与河流有直接水力联系时的基流分割图示

（4）示踪法

水中的氢氧稳定同位素常用来示踪地下水与地表水的相互交换量，通过河道水量均衡方程和质量平衡方程的联合求解，计算地下水与河流交换量（Scanlon，2002）。水量均衡方程和质量平衡方程如下：

Qup+∑Qin+Qgi=Qdown+∑Qout+Qgo+Er （3-10）

Qup·δup+∑Qin·δin+Qgi·δgi=Qdown·δdown+∑Qout·δout+Qgo·δgo+Er·δEr　（3-11）

式中：Qup，Qdown分别为上、下游断面河流量（m3/s）；Qin，Qout分别为测流断面间各支流的流入、流出量（m3/s）；Qgi，Qgo分别为测流断面间地下水流入、流出量（m3/s）；Er为测流断面间河道水蒸发量（m3/s）；δup，δdown分别为上、下游断面河水氢氧稳定同位素δ值（‰）；δin，δout分别为测流断面间各支流的流入和流出水的氢氧稳定同位素δ值（‰）；δgi，δgo分别为测流断面间地下水流入、流出水的氢氧稳定同位素δ值（‰）；δEr为河道蒸发水的氢氧稳定同位素δ值（‰）。

（3-11）式中的δEr通常难以测定，Krabbenhoft（1990）给出了计算公式：

区域地下水功能可持续性评价理论与方法研究

式中：δL，δα分别为地表水和大气水汽同位素含量；h为相对湿度；α′为水-气界面温度下同位素平衡分馏因子，等于1/α；ε为总分馏因子，ε=1000（1-α′）+Δε；Δε为动力学分馏因子，对于δD，Δε=12.5（1-h），对于δ18O，Δε=14.2（1-h）（Gonfiantini，1986）。

平原区的河流往往在上游地带渗漏补给地下水，在下游地带接受地下水补给，因此，在实际应用时，通常需要首先确定该两种情况发生的分界面，然后分段进行计算。

3.地下水侧向流入流出量

一般采用达西定律计算，公式为

Qg=v·B·M=K·J·B·M （3-13）

式中：Qg为地下水侧向流入流出量（m3/d）；v为地下水渗流速度（m/d）；B为过水断面宽度（m）；M为含水层厚度（m）；J为地下水水力梯度（无量纲）。

计算时，需要注意：①采用测流计或地下水示踪技术测定地下水流速时，所测定的流速是地下水的实际流速（va）。②采用达西定律计算时，需要将实际流速换算成渗透速度，即v=va·n（n为有效孔隙度）。③不同地段含水层的孔隙度和地下水水力梯度不同，不同时段地下水的水力梯度也不相同。因此，实际应用时，应分地段、分时段分别进行计算。

4.地下水与湖泊或水库的交换量

（1）水量均衡法

计算公式为

Qlr=Qgi-Qgo=ΔV-P·F+E·F-Qsi+Qso （3-14）

式中：Qlr为地下水与水库或湖泊的净交换量（m3/a），Qlr＞0，则地下水向湖泊排泄量大于湖泊向地下水的渗漏量，湖泊接受地下水的净补给；Qgi为地下水向湖泊的排泄量（m3/a）；Qgo为湖泊向地下水的渗漏量（m3/a）；ΔV为水体体积的年变化量（m3/a）；P为年降水量（m/a）；F为水体水面面积（m2）；E为水面蒸发量（m/a）。Qsi为地表水年流入量（m3/a）；Qso为地表水年流出量（m3/a）。

（2）示踪法

采用式（3-14），一般只能获得地下水与湖泊之间的净交换量，为了分别求取Qgi和Qgo，可以利用水体中的天然示踪剂建立质量均衡方程（Sacks，1998）：

Qgi·Cgi-Qgo·CL=ΔV·CL-P·F·CP+E·F·CE-Qsi·Csi+Qso·CL　（3-15）

式中：Cgi为补给湖泊地下水的示踪剂浓度；CL为湖泊水的示踪剂浓度；CP为降水的示踪剂浓度；CE为湖泊蒸发水的示踪剂浓度；Csi为补给湖泊地表水的示踪剂浓度；其他符号同式（3-14）。

联合求解（3-14）和（3-15）两个方程，即可获得Qgi和Qgo。常用的天然示踪剂是水中的氢氧稳定同位素（Sacks，1998；Scanlon，2002），这时需要知道蒸发水的氢氧稳定同位素值，其计算方法见式（3-12）。

5.越流量（包括补给和排泄）

计算公式为

Qy=F·K·J （3-16）

式中：Qy为越流量（m3/d）；F为计算面积（m2）；K为弱透水层的垂直渗透系数（m/d）；J为弱透水层上下含水层间的水力梯度（无量纲）。

6.凝结水补给量

可根据均衡试验场地中渗透仪的观测资料求得，但是计算时应注意将观测中冬季潜水冻结层融化的水量扣除。

7.渠系渗漏补给量

根据渠系衬砌状况，选用实测或经验系数计算。若渠道没有任何衬砌，其渗漏补给量与河道渗漏补给量计算方法相同。若渠道有衬砌，则可采用如下公式计算：

Qci=r·（1-η）Qc·Δt （3-17）

式中：Qci为渠道渗漏量（m3）；r为渠道渗漏修正系数（无量纲）；η为渠系有效利用系数（无量纲）；Qc为渠道过水量（m3/s）；Δt为计算时段（s）。

8.田间灌溉入渗补给量

（1）入渗系数法

计算公式为

Qsi=β·Qs·F·N （3-18）

式中：Qsi为田间灌溉入渗量（m3）；β为入渗系数（无量纲）；Qs为灌溉定额（m3/m2）；F为灌溉面积（m2）；N为灌溉次数。

（2）水量均衡法

根据水均衡原理，用灌溉量减去排放量、蒸发量和其他消耗量计算。

（3）地中渗透仪法

在田间专门设置地中渗透仪，直接测定灌溉水渗漏补给量。

9.潜水蒸发蒸腾量

（1）蒸发系数法

计算公式为

Qe=E·c·F （3-19）

式中：Qe为潜水蒸发蒸腾量（m3/a）；E为水面蒸发量（m/a）；c为潜水蒸发系数（无量纲）；F为计算面积（m2）。

（2）经验公式法

通常利用经验公式求出潜水蒸发强度（ε），然后按下式计算：

Qe=ε·F （3-20）

式中：ε为潜水蒸发强度（m/a）；F为计算面积（m2）。

潜水蒸发强度一般采用柯达夫-阿维利扬诺夫公式计算，即

区域地下水功能可持续性评价理论与方法研究

式中：λ为植被修正系数（无量纲）；h为潜水水位埋深（m）；h0为潜水蒸发的临界深度（m）；θ为无量纲指数，因气候和土壤而异，取值1～3，一般可以取1；其他符号同前。

10.储变量计算

计算公式为

ΔS=μ·F·ΔH （3-22）

式中：F为计算面积（m2）；ΔH为水头变化（m）；μ为给水度（潜水）或储水系数（承压水）。

（四）若干问题说明与适用条件

1.参数的获取

以上介绍的方法中涉及各种参数（渗透系数、导水系数、给水度、储水率、储水系数、孔隙度、垂向渗透系数、越流系数、降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等），按“全国地下水资源及其环境问题调查评价技术要求系列（一）”的要求获取。

2.参数分区

水量均衡法是一种集中参数系统的方法，在地下水数量评价时往往难以满足计算精度要求，尤其是在区域地下水数量评价中。地下水系统是一个复杂的非均质系统，各种参数（渗透系数、导水系数、给水度、储水率、储水系数、孔隙度、垂向渗透系数、越流系数、降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等）是空间位置的函数，有些还是时间的函数（如降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等），所以为了提高计算精度，需要综合考虑各种参数的时空变化特征，在空间上将地下水系统划分成若干子区块，在时间上划分成若干时段，在子区块和各时段可以认为各种参数是一个相对稳定的数值，然后分别计算各个区块和时段的水量，最后集成总量。理论上讲，划分的区块越小、时段越多，计算精度就越高，但是工作量就越大。在实际应用时，应根据评价区所拥有的资料状况和计算精度要求进行适当的划分。

3.点参数的区域化

通过各种方法获得的参数，大多是点源数据。在区域地下水数量评价时，需要将点源数据转化成区域数据。点源数据的区域化，常采用的方法是Kriging插值法。Kriging法是一种最佳空间估计法，其本质是最佳无偏估计，是对空间分布的数据求线性最优、无偏内插估计的一种方法（Gress N.A.C.，1990，1991；Dentsch C.V.，1992；Gelhar L.W.，1993）。常用的软件中都有Kriging插值功能模块，如Surfer、MapGIS等，也有一些文献中给出了计算程序（徐士良，1995；Dentsch C.V.，1992）。

4.适用条件

水量均衡法方法原理明确，计算公式简单，计算精度高低可调，适应性强。但是，在补给、排泄条件复杂的地区，涉及的均衡要素较多，某些均衡要素难以准确测定或求取成本和工作量较大，计算精度不如数值法高。

水量均衡法既可用于区域地下水数量评价，也可用于局域地下水数量评价或水源地评价；既可评价地下水补给资源量，又可评价可采资源量，是最常用、最基本的地下水数量评价方法，其成果是其他方法的验证依据之一。在补给、排泄条件简单，地下水系统边界比较清楚，水均衡要素容易确定的地区，应用效果较好，评价结果精度高。

  对于求解无约束优化模型 通常会有下面的一种迭代步：

  通过某种搜索方法确定步长因子 ，使得   这实际上是目标函数 在规定的一个方向上移动所形成的单变量优化问题，也就是所谓的 “线搜索” 技术，令 这样，搜索式等价于求步长 使得

   精确线搜索的基本思想： 首先确定包含问题最优解的 搜索区间 ，然后采用 某种插值或分割技术缩小区间 ，进行搜索求解。

  由于不少精确线搜索算法都是针对单峰函数建立起来的，下面介绍一种确定搜索区间并保证具有 近似单峰性质 的数值算法—— 进退法

算法1：进退法

S1： 给出 , ,

   计算 ,

S2： 如果 则转 S4

   如果 则转 S6

S3： 令 , 如果 则转 S3

   如果 则转 S6

   , 转 S7

S4： 令 , 如果 则转 S7

   如果 则转 S5

   令 , 转S4；

S5： 令 , 如果 , 则令 , 转 S7

   如果 , 则令 , 转 S7

   令 , 转 S5.

S6： 令 , 转 S2

S7： 令 , 即停。

例1： 利用进退法求解极值区间实例。取初始点 ，步长 ，用进退法求函数 的极值区间。

  黄金分割法又称 0.618 法，其基本思想是通过试探函数值的比较，使得包含极小点的搜索区间不断缩小，直至区间内的值足够接近极小值为止。该方法仅需要计算函数值，适用范围广，使用方便。下面我们来推导 0.618 法的计算公式。

  设   其中 是搜索区间 的单峰函数。

  假设在第 次迭代时搜索区间为 , 为确定下一次迭代区间，我们取两个试探点 , 计算得到 和 , 存在两种情况：

  （1）若 , 则令

  （2）若 , 则令

  我们要求两试探点满足如下条件：

  （1） 与 长度相同，即

  （2）区间长度的缩短率相同，即 .

  从而可得   由于两区间长度一致，不妨假设新的迭代区间为 为进一步缩小区间，取新的试探点 , 得 若令 则 这样，新的试探点就不用重新计算。得到区间的缩短率为

算法2：黄金分割法（0.618法）

S1： 选定区间 及精度 , 令 , 计算试探点

  

  

S2： 若 ,则停止计算

   否则，当 时转 S3；当 时转 S4；

S3： 令 , 转 S5

S4： 令 , 转 S5

S5： 令 , 转 S2.

例2： 利用黄金分割法求解极值实例。利用黄金分割法求下面函数的极小值

  基本牛顿法是一种是用导数的算法，它每一步的迭代方向都是沿着当前点函数值下降的方向。其基本思想是：用 在探索点 处的二阶 Taylor 展开式得到的二次函数 来近似代替 ：   其中， 可以认为是 的近似。因此，求函数 的极小值点近似于求解 的极小值点， 函数 应该满足一阶必要条件 ：   即   上式即为牛顿法的迭代公式，当 时，对于区间内的 都成立；反之当 时，牛顿法可能收敛到极大值点。

算法3：基本牛顿法

S1： 给出 精度 , 令

S2： 若 , 停止, 极小值点为 ；

S3： 令

S4： 令 , 转 S2.

例3： 取初始点 ，用牛顿法求 的任一极小值点。

  线搜索技术是求解许多优化文体下降算法的基本组成部分，但精确线搜索往往需要计算很多的函数值和梯度值，从而耗费较多的计算资源。特别是当迭代点远离最优点是，线搜索方法通常不是十分有效和合理的。因此，既能保证目标函数具有可接受的下降量又能使最终新城的迭代序列收敛的非精确线搜索越来越流行。常用的有 Walf 准则与 Armijo 准则 。

  牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性，但基本牛顿法初始点需要足够“靠近”极小点，否则有可能导致算法不收敛。这样就引入了 阻尼牛顿法（也称全局牛顿法） ，阻尼牛顿法最核心的一点在于可以修改每次迭代的步长，通过沿着牛顿法确定的方向一维搜索最优的步长，最终选择使得函数值最小的步长。

  阻尼牛顿法是 基于 Armijo 准则的搜索 ，满足 Armijo 准则：   一般的，可取 为 或更小的值。由于 时下降方向，不等式右边是关于 的线性减函数。因此只要 不取的太小，不等式可以保证新迭代点 的函数值较 有一定量的下降。

算法4：阻尼牛顿法（全局牛顿法）

S1： 给出 精度 , 令

S2： 计算 , , 若 则转 S4

   若 则停止

   令初始值

S3： 令 , 如果 , 则转 S3

   令 , , 转 S2；

S4： 令 , 如果 , 则令 , 置

S5： 如果

则转 S6；置 ； 转 S5

S6： 令 , , 转 S2.

例4： 取初始点 ，用全局牛顿法求函数 的极小值点。

M=10%产生M行N列的随机数矩阵

N=8

miu1=1%第一个分布的参数

sigma1=2%第一个分布的参数

miu2=6%第二个分布的参数

sigma2=1%第二个分布的参数

R = 0.2*normrnd(miu1,sigma1,M,N)+0.8*normrnd(miu2,sigma2,M,N)

单点的概率全是0，那你取出来的随机数算什么？

若干个随机数要满足统计分布，是要按区间统计的

另外我不知道你要做什么就是了。

你如果想按一定的概率密度来产生随机数，你最好用反函数法之类的来弄。

比如产生一个x.^2分布的随机数，不过这些要归一化。

============================================

首先，我知道我的是错的了。如下图就可知

M=1000%产生M行N列的随机数矩阵

N=1

miu1=1%第一个分布的参数

sigma1=2%第一个分布的参数

miu2=6%第二个分布的参数

sigma2=1%第二个分布的参数

R = 0.2*normrnd(miu1,sigma1,M,N)+0.8*normrnd(miu2,sigma2,M,N)

x=-5:0.001:15

y1=normpdf(x,miu1,sigma1)

y2=normpdf(x,miu2,sigma2)

subplot(2,2,1)

plot(x,y1)

subplot(2,2,2)

plot(x,y2)

subplot(2,2,3)

y3=0.2*y1+0.8*y2

plot(x,y3)

subplot(2,2,4)

dx=0.5

xx=-5:dx:15

yy=hist(R,xx)

yy=yy/M/dx

plot(x,y3)

hold on

bar(xx,yy)

=======================================

正确做法，我还没弄出来，继续中。。。。

============================================

_____________________新的尝试

下面的结果我觉得可能可以接受。

思路：基于反变换法

Matlab下面有

p=normpdf(x,miu,sigma)是求出x处的概率密度。

p=normcdf(x,miu,sigma)是求出X<x的累积概率密度(就是从负无穷大到x处的概率密度的积分)

我给定一个区间，这个区间外的概率我认为是0（这一点不够严谨，理论上应当是从负无穷到正无穷）

我这里取的是-10:15，其间我取了25000个点，求出这些点的累积概率值(两个的加权和y3)，记这个为F(x)，根据反变换法，

F(x)=u,其中u是一个0到1的均匀随机数。只要求出它的解x0,那么x0就满足所给定的概率密度分布。这里我用的是插值。用

（y3,x）来插值出u所在的位置

声明，这里有一些地方不够严谨，严谨应当用解析的方法来做反变换。

%%%%%下面是程序

M=1000%产生M行N列的随机数矩阵

N=1

miu1=1%第一个分布的参数

sigma1=2%第一个分布的参数

miu2=6%第二个分布的参数

sigma2=1%第二个分布的参数

x=-10:0.001:15

y1=normpdf(x,miu1,sigma1)

y2=normpdf(x,miu2,sigma2)

y3=0.2*y1+0.8*y2

y1=normcdf(x,miu1,sigma1)

y2=normcdf(x,miu2,sigma2)

y=0.2*y1+0.8*y2

u=rand(N,M)

R=interp1(y,x,u,'linear')

dx=0.5

xx=-10:dx:15

yy=hist(R,xx)

yy=yy/M/dx

bar(xx,yy)

hold on

plot(x,y3,'r*')

---