第七章 土壤特性的空间变异性2

第七章 土壤特性的空间变异性

土壤特性在空间分布是非均一的，例如在平面上土壤的质地、剖面上土层的厚度，以及土壤含水率、土壤水分运动参数、土壤水基质势（负压）、含盐量等在同一时刻，即使是相距很近的点，其数值也是不同的。这种土壤特性在空间上分布的差异性称为土壤特性的空间变异性。

为了探讨土壤各种因素的变化规律，必须进行田间试验，布设观测点和取样点，由于土壤特性空间变异性的存在，观测和取样点数目不宜过少，但因受人力、物力的限制，也不宜过多，这样，就存在确定合理取样数目、对未观测点进行估值、利用土壤特性的变异规律对田间土壤水分运动进行分析等问题。以下将分别对这些问题作简要介绍。

第一节 土壤特性的变异性和合理取样数目

一、土壤特性的变异性分析[55

，63，68，69]

土壤特性受随机因素的影响，在一定的空间内进行多次试验所取得的数值可能是不同的，因而存在一定的偶然性。如将土壤特性看作是一个随机变量Z ，在空间上变化看作是独立的变化，则其变化特征可由其概率密度函数产P （x ）表示。

X 为土壤特性x 的可能取值，发生这一事件的概率为p （x ）。发生随机变量X 的取值小于或等于X 的事件的概率为

?∞

-=≤x

dx x p x X p )()( （2-7-1）

概率P （X ≤x ）称为累积概率，概率密度函数P （x ）为一非负函数，P （x ）≥0，且 P （－∞〈x 〈∞〉=

?

∞

∞

-dx x p )(。累积概率函数P （X ≤x ）。有时写成：

)()(x P P x F ≤=

F （x ）称为随机变量X 的分布函数。

随机变量有多种分布形式，对于土壤特性最常见的有以下两种。 1．正态分布

在这种情况下概率密度表达式为

2

2

2)(ex p[21

)(σσπm x x p --=

] (2-7-2) 式中δ、m ——常数

随机变量X 服从正态分布，记为X ～N （m ，δ2）。在m=0，δ=1时的正态分布，称之为“标准正态分布”，记为 N （0，l ）。

??

∞

-∞

---==x x

m x dx x p x F 2

2

2)(ex p[21

)()(σσπ] (2-7-3) 2．对数正态分布

在数据取对数后服从正态分布的，其概率密度的表达式为

2

2

2)(ln ex p[21

)(**--*=

σσπm x x x p ] x>0 式中δ*、m*—一对数正态分布的特征常数。

判断随机变量X 是否属于正态分布或对数正态分布，可以根据试验或观测值或对观测数据取对数计算累积概率，并点绘于正态概率纸上，如观测参数值与相应累积概率呈直线则为正态分布。如对数值与累积概率呈直线关系，则为对数正态分布。若随机变量X 为正态分布，式（2－7－2）中m 和δ。分别为随机变量X 的均值和方差；若随机变量为对数正态分布，则m*和δ*为Inx 的均值和方差。如随机变量X 有一个容量为N 的样本：x 1，x 2，…，x 。，其分布属正态分布，该随机变量总体的均值和方差可用样本的均值x 和方差S 2进行估计：

∑==

N

i i

x

N

x 1

1

(2-7-5)

1

1

2

-=

N S 21

)(x x

N

i i

-∑= (5-7-6)

同理，随机变量为对数正态分布时，m*和δ*也可用x*和 S*2进行估计。 标准差（方差δ2的平方根）与均值之比称为变差系数C v ，：

m

C v σ

=

(2-7-7)

土壤特性的变差系数Cv ；可以反映土壤特性变异性的大小。

二、合理的取样数目

在根据一定容量N 的样本分析土壤特性时，样本的均值x n 和方差与总体的均值m 和方差是有一定差别的。如将x n 也作为随机变量，则取样数N 越大，x n 越接近于m ，x n 的方差越小。土壤特性分析的要求一定，即样本的均值x n 和总体的均值m 之差必须小于或等于一定精度μ时，则达到这一精度要求的取样数目N 必须使发生小平或等于这一精度μ的事件的概率达到所要求的置信水平，即

L n p m x p =≤-}|{|μ (2-7-8)

在取样数目足够多时，根据概率统计原理可知，随机变量N

m

x u N /2

σ-=

为标准正态分布（即均值为0，方差为1）。根据正态分布特点

在P L 已知时，可自正态分布双值分位数表查得满足置信水平P L （显著水平a=1一P L ）时

的

N

m

x N /2

σ-值μα即

P u L a N N m x p =??

????????≤-/2

σ (2-7-9)

自式（ 2－ 7－9）可求得满足置信水平 P L 和一定精度μ要求的取样数目：

2

2???

? ??=μσα

u N (2-7-10) 例如，置信水平P L =95％时，μα查得为1.96，则

2

84.3???

?

??=μσN

(2-7-11)

若取μ=km （k 可取5％、10%等），由于Cr=δ/m ．所以取样数：

2

84.3??

?

??=m C N V (2-7-12)

当k 取10%．Cv=0.l 时，合理取样数N=4；

Cv=l.0时，合理取样数N=40

在实际工作中，总体方差是未知的，须用样本方差S 2代替。由概率统计原理可知，随机变量t=（xn 一m ）／√S 2/N 服从t 分布。满足置信水平P L ，（显著水平a=1—P L ）的t αv ，可自t 分布函数累积概率表2－7－1查得。

P t L a N N S m x p =??

????????≤-ν,2

/

（2-7-13）

式中t αv 一—当显著水平a=l 一p L ，自由度v=N —1时的t 值。

自式（2－7－13）可得达到精度μ要求的取样数目：

2

2,???

? ??=μναS t N （2-7-14）

由于t 为取样数 N=v ＋l 的函数[70]须通过试算求得。

例如，根据试验资料，土壤特性标准差为S=0. 05，要求试验精度为μ=0.01，显著水平α=0.1。达到以上要求的取样数N 值为

2

,2

2

2501.005.0ναναt t N =??

? ??=?

70,678.11,1.0==-N t N

第二节 土壤特性的空间结构和Kriging 内插法

在土壤特性空间变异性分析中，常将土壤特性参数在空间上的变化看作是随机的、互相独立的。但是，实际上在一定的范围内各点的参数值存在着一定的相关关系，只有在取样点间的距离超过一定数值时，各点取样才能认为是独立的。以下将分析土壤特性空间分布的相关性，亦即土壤特性的空间结构。

一、土壤特性空间分布的自相关分析

在进行土壤特性参数测定时，如沿某一方向以等间距Δx 布置测点，总测点数为N 。测点位置分别以x 1，x 2，x 3…，x n 表示。

如图2－7－l ，测得各点参数值为Z （X 1），Z （X 2），…，Z （x N ）。两组共N －1对相对应的参数系列，对这两组参数值进行相关分析，所得的相关系数称为间距为h=ΔX 的自相关系数。如取间距h=2ΔX ，将形成两组共N —2对相互对应的参数系列，可求得间距为h=2ΔX 的自相关系数，其余类推，可得间距为任一间距h 时的自相关系数r （h ），其表达式为

()()()[]

()[]()[]

h x Z D x Z D h x Z x Z C h r ov ++=

, （2-7-15）

()()[]()[]μμ-+-=

∑-h x

Z x Z N

h C i

N

i i

ov 1

1

（2-7-16）

式中D[ Z （X ）]、D[ Z （X+h ）]——分别为随机变量Z （X ）和Z （X 十h ）的方差；

C 0v [Z （X ），Z （X 十h ）]——这两个随机变量的协方差，可简写为C 0v （h ）；

μ——参数Z （X ）的均值。

在两组系列足够长时

()[]()[]2σ=+=h x Z D x Z D

(2-7-17)

可简化为

()()

2

σ

h C h r ov =

(2-7-18)

自相关系数r （h ）随间距h 而变化，故又称自相关函数。h →0，r （h ）=1。r （h ）随 h 的增大而减小，在不存在相关关系时，r （h ）=0，参数值是互相独立的。若土壤特性参数Z （x ）的测定不是沿一个方向，此时间距

h 为矢量，仍可沿不同方向采用式（2－7一15）进行自相关分析，只是式中为矢量。

二、土壤特性空间分布的半方差分析

反映土壤特性空间结构的另一指标为半变异函数或半方差，进行空间变异分析时有以下两项基本假设。 1．均值稳定

认为土壤特性参数Z （x ）的均值E[Z （x ）]存在，且为常数Δ，即

()[]()[]?=+=h x Z E x Z E

(2-7-19)

根据这一假定可推论Z （x ）与Z （X 十h ）的协方差存在，且为有限值，即

()()[]()[]()[]{}?-+?-=+h x Z x Z E h x Z x Z C ov ,

()()[]2

?-+=h x Z x Z E

=()h C ov

(2-7-20)

2．D[Z （x ）—Z （X 十h ）]存在，且为有限值

两系列Z （x ）与Z （X 十h ）对应值之差的方差风D[Z （x ）—Z （X 十h ）]仅与h 有关，记为2γ（h ），即

()()()[]()()[]{}

2

2h x Z x Z E h x Z x Z D h +-=+-=γ (2-7-21)

由于γ（h ）= D[Z （x ）—Z （X 十h ）]/2，故称半方差；γ（h ）随h 而变化，有时亦称

为半变异函数。

将式（2－7－21）右端展开，并利用式（2－7－19）和式（2－7－20），可导出：

()()()h C C h r ov ov -=0

其中 ()()[]{}

2

0?-=x Z E C ov

（2-7-22）

将式（2－7－18）代入式（2－7－22），并利用γ（0）=l ，则可导出半方差与自相关系数理论上关系：

()()[]h r h -=12

σγ （2-7-23） h ＝0，r （0）＝1，γ（0）=0随着h 的加大，半方差γ（h ）也随之增大。当h ≥a 时，r （h ）

=0，则，r （h ）=δ2（或 S 2）。因此，由半方差图可以判断出该参数空间分布的相关距离 a 。 在根据实测资料进行土壤特性变异性分析时，半方差值可自下式计算：

()()()[]{}

2h x Z x Z E h +-=γ

()()[]2

1

21∑=+-=N

I i i h x Z x Z N (2-7-24)

在求得不同间距h 时的半方差γ（h ）后，常将γ（h ）与h 的关系用经验公式表示。例如，雷志栋等根据商丘县大吴庄5．2km 2范围内 46眼观测井 1979年 10月 1日地下水位观测资料，求得γ（h ）与h 关系点据如图2－7－2所示，并将γ（h ）与h 关系概化为以下经验公式[67]： ()h a

C h 1

=

γ a h <

()1C h =γ

a h ≥

(2-7-25)

式中α—一为相关距离（如图2－7一2中所示）；

C 1——取观测值的方差（图2－7－2）。

三、Kriging 最优内插法

研究变量Z （x ）有N 个测点。x 1，x 2，x 3，…，x N 。为已知，需求在x 。点Z （x 。）的估计值。假定Z （x 。）*的估计值是N 个有效数据的加权平均值：

()()∑=*

=N

i i i x Z x Z 1

0λ

(2-7-26)

为了找出一组能达到最优估计的权因子λ，必须满足无偏估计和方差最小两项基本要求 （1）无偏估计，即无系统误差，可表示为

()()[]

000=-*

x Z x Z E

(2-7-27)

将式（2-7-26）代人式（2－7－27），有：

∑==N

i i

1

1λ

(2-7-28)

（2）估计方差最小，即

()()[]

()()

[]{}min 2

0000=-=-*

*x Z x Z E x Z x Z D

(2-7-29)

式中方差引入了线性组合()()∑=*

N i

i i x Z x Z

λ0为了区别一般方差，故称之为Krige 方差。

根据空间变异分析理论．式（2－7－29）经推导可得：

()()[]

()011

0,x x x Z x Z D i j N

i N

j i o γλλ∑∑==*

-=-

(2-7-30)

式中：λ（x i ，x j ），λ（x i ，x 0）分别代表向量的一端扫过x 0和x i ，另一端独立地扫过点x i 时半变异函数（半方差）值，二者均可由试验资料采用式（2－7－24）求得的λ（h ）半方差图或与公式（2－7－25）相类似的经验公式求得。

为使估值最优，所选取的系数λi 应满足式（2－7－30）和式（2－7－28）。为此，选用拉格朗日乘数法，令

()()()[]

??? ??---=∑

=*

n i i i x Z x Z D L 10012,λμμλ

(2-7-31) 式中μ一待定的拉格朗日乘子。

取L 对λi ，μ的偏导数，并使其为零。由δL/δλi =0得：

()()j

i

j

i

N

j j

x x x x ,,1

γμγλ=+∑=

（2-7-32）

由0=??i

L

λ，得： 11

=∑=N

i i

λ

（2-7-28）

于是可得一个N+1个未知数（N 个λ，一个μ），有N 十1个方程的方程组，写成矩阵形

式为：

求解矩阵（2-7-33），即可求得各项系数λi 及μ。利用式（2－7－20）和式（2－7－32）可得估计方差：

()()[]

()μγλ+=-=∑=*

N

i i i x x x Z x Z D 1

000,

第三节 田间土壤水分特征曲线与导水率的标定

土壤水分特征曲线——土壤水基质势和含水率之间的关系ψ～θ和非饱和土壤导水率与含水率之间的关系k ～θ，都是分析土壤水分运动问题的基本参数，这些参数都具有明显的空间变异性。所谓“标定”即是通过对每一点选取适当的比例系数（或标定系数），将空间变异的ψ（θ）和k （θ）关系，标定为对各点土壤均适用的ψ～θ和k ～θ关系。这样就可用统一的ψ～θ、k ～θ关系和随各点变化的标定系数a 代替各点均不相同的ψ～θ和k ～θ关系。以下简要介绍标定的基本概念和方法。

一、田间土壤水分特征曲线的标定

（一）土壤水基质势ψ的标定

土壤水基质势ψ（常以水头h 表示）是由土壤介质对土壤中水分的吸持作用而引起的。这种吸持作用常近似地简化为水气界面的毛管作用，可表示为

r

σ?2-

= （2-7-35）

式中 δ——水的表面张力系数；

r ——毛管半径或弯液面的曲率半径。

若在两不同点上的土壤，其孔隙大小、形状及其分布在几何上相似，如图2－7－3所示，如土壤中水分状态也是相似的，则

θθθθθθ??2

121====或 （2-7-36）

式中 θ、θ*——分别为各点的含水率与饱和度；

θ、θ——相应θ、θ的平均值。

处于几何和状态相似的各点土壤，尽管含水率或饱和度相等，但弯液面的曲率半径r 并不相同。因此，在同一含水率和炮和度情况下，各点的基质势ψ1，ψ2，… 可能是不同的。现对各点分别取一特征长度λ，如图2－7－3所示。由于处于相似状态，因而：

??==

2

2

1

1

r r λλ

自式（2-7-35）可得：

1

1

112r λσ

?λ-=

2

2

222r λσ

?λ-=

……………………

由于各点表面张力系数δ是相同的，则可导出： λ1ψ1=λ2ψ2=…=λ*ψ*

式中λ*——各点特征长度的平均值或标定值；

ψ*——各点基质势的平均值或标定值。

对任一点i ，可以写出：

??λ

λ?=i i

令a i =λi /λi 称为标定系数，则

?

?α?=i i

(2-7-37)

式（2－7－37）即为基质势的标定公式，该式表明，当各测点土壤相似时，在某一饱

和度θ下．将各测点的基质势ψi 乘以该点的标定系数a i ；后所得的值均相等，该值即为在此饱和度下的基质势的标定值ψ0。各点的标定系数αi =λi /λ不同，但在同一点上αi 则是一个不随含θ而变化的常数。

（二）田间土壤水分特征曲线的标定 以上根据土壤的几何相似性说明了土壤水基质势标定的概念。实际上，田间各点土壤几何上并不完全相似，也不可能直接量测土壤各点的特征长度，从而计算各点的标定系数，田间土壤水分特征曲线必须根据实测资料进行标定。

以下介绍由实测资料进行田间土壤水分特征曲线标定的方法。为应用方便，含水率θ均用饱和度Θ代替。

在田间就地测定或取原状土样室内测定不同饱和度下的土壤水基质势时，设有N 个现

场测点或土样，编号为I=1，2，…，N 。各点实测的饱和度Θ和相应基质势ψ资料的数目可以是不同的，分别记为J1，J2，…，JN ，故总记有Eni=1 Ji 对(ψ,θ)实测数据。在任一点i 处，任一状态j 的饱和度记为Θij ，相应的基质势的实测值记为ψij 。

田间土壤水分特征曲线标定的任务，就是要确定一条适用于各点的标定水分特征曲线ψ～Θ和田间各点的标定系数 a i （I= 1，2，…，N ）。为了使用和求解的方便，一般事先确定标定水分特征曲线的经验关系式，目前通常采用的是：

()()()[]

θθθθ?n n a a a -+?+-+-=111?221 (2-7-38)

其中 a 1, a 2, …, a N 为待定常数， n 一般取4即可。这样，标定问题就变成根据实测资料求标定系数a 1, a 2, …, a N 和待定常数a 1, a 2, a 3, a 4的问题。

定义方差SS 为

()2

11

?∑∑==-=N i L

j ij i ij SS ?α?

(2-7-39)

式中—ψij 按式（2－7－38）由Θij 计算所得的值。

当各测点的土壤在几何上完全相似时，由于ψij ,=αi ψij ，式（2－7－39）表示的方

差值SS=0。实际上，各测点的土壤并不完全相似，因此标定系数和待定常数的选取，必须使SS 值有最小值。具体进行标定时，可用迭代法或一次法。

1．迭代法

迭代法是Warrick （1977）提出的[64]，求解步骤如下：

（1）先估计一组由式（2－7－38）表示的ψm ～θ关系中的待定常数值a 1, a 2, a 3, 和a 4。令标定系数a 1=a 2=a 3=…= a n =1（未标定），由式（2-7-39）计算SS 值，记为SSA 。 （2）根据已知待定常数求标定系数。在迭代求解过程中，若不加限制必然会得出待定常数和标定系数全为零的结果。尽管此时SS=0，但毫无意义。为此，令所求的标定系数满足下列约束条件：

N a a a N =+?++21

(2-7-40a)

上述标定系数中的N —l 个可视为独立变量，而余下的一个，如a n ，则可表示为其他标定系数的函数。即

121--?---=N N a a a N a

(2-7-40b)

这样，已知待定常数求标定系数的问题便归结为满足式（2－7－40a ）的约束条件下，求以式（2－7－40）表示的方差SS 为最小值的问题，即所求标定系数应满足：

1,2,1,0-?==??N P SS

P

，α (2-7-41)

将成（2－7－40a ）代入上式，并利用由式（2－7－40b ）得出的关系：

()1;0;1-=??≠=??=??P

N P i P P

P i αααααα 方程式（2－7－41）可展开为

()()0?2?21

1

=---∑∑==P

N

J j J j Nj N Nj pj pj p pj

?α?

??α?

(2-7-42)

即 0??1

1

1

21

2

=+--∑∑∑∑====Jn

j Nj Nj J j pj Pj JN N Nj

N J j pj

p P P a ????

??α

1,,2,1-?=N P

(2-7-43)

上式也可简写为 N P N p p B B A -=-αα

1,,2,1-?=N P

其中

N P B A Jn

j Nj J j Pj Pj p Jn

j Nj

J j Pj

P P P

?=???

?

?

??

==∑∑∑∑====,2,1?1211212

???

?? (2-7-45)

由实测资料，用式（2－7－45）可求得各A p 值。在待定常数a 1, a 2, a 3, a 4已知的情况下，根据实测资料的θpj 值，由式（2-7-38）可计算相应的ψpj 值，进而由式（2－7一45）亦可求得各B p 值，有了以上关系之后，则可进一步导出求标定系数的公式。

以P=1代入式（2－7－44），得：

N N B B a a A -=-111

(2-7-46)

由式（2－7－44）和式（2－7－46）可得到由a 1计算a P 的公式：

1,,2,1,111-?=--=

N p A B B A a a P

p

p

(2-7-47)

将式（2－7－47）应用于式（2－7－40a ），则有：

∑

==--N

p P

p

N A B B A a 1

111

由上式可解出a 1

()()

()

∑∑∑=-=-=-+=

N

P P N p N p p p P A A A B A B N a 1

111

11

1111 (2-7-48a)

如利用A N =1，上式还可表示为

()()

()

∑∑∑-=-=-=++-+-+=

1

2

1111

11

111111N P P N p N p p p P N A A A A B A B B B N a (2-7-48b)

标定系数a 1由式（2－7－48a ）求得后，利用式（2－7－47）便可求解马a 2 ,a 3 ,…,a N .然后由式（2－7－39）求得相应的方差SS ，记为SSB 。

（3）根据已知标定系数求待定常数。标定系数已知后，式（2－7－39）表示的方差SS 便是待定常数a 1, a 2 ,a 3和a 4的函数。根据SS 取最小值的条件，分别对a P （P=1，2，3，4）求导，并使之等于0，不难得出求解待定常数方程和待定常数应满足的条件

?????

???????=????????????????????????45352515432144434241343332312423222114131211b b b b a a a a b b b b b b b b b b b b b b b b （2-7-49）

其中 ()()∑∑

==--=N

i J j ij p ij

k

ij

Pk i

b 11

2

11θθθ

（2-7-50）

4,3,2,1,4,3,2,1==k p

()∑∑

==-=N i J j ij

p

ij

P i

b 11

51θθ

（2-7-51）

4,3,2,1=p

根据实测资料及已知的标定系数，由式（2－7－50）和式（2－7－51）可计算方程组（2－7－49）中的各项系数及常数，解代数方程组（2－7－49）便可求得标定水分特征曲线经验表达式中的待定常数a 1, a 2 ,a 3和a 4。

为了推求标定水分特征曲线和各点的实测水分特征曲线的方差值，需将a 1=a 2=…=a N =1代人式（2-7-39），所得SS 值记为SSA 。

（4）重复以上（2）和（3）步骤，直到前后两次所得SSB 值之差满足要求为止。最后所得标定系数a 1, a 2 ,…，a N 和待定常数A 1，A 2，A 3 ,A 4即为所求。

2．一次法

一次法是一种相对简单的标定方法[67]，在利用这种方法进行土壤水分特征曲线标定时，首先a 1=a 2=…=a N =1（未标定），根据实测结果由式（2－7－50）和式（2－7－5l ）可计算各系数项和常数项的值，解方程组（2－7－49）便可求得待定常数a 1, a 2 ,a 3和a 4，由此可得标定水分特征曲线ψm ～θ[式（2－7－38）]。按式（2－7－39）计算SS 值，此即标定水分特征曲线与各实测水分特征曲线之间的方差值，记为SSA 。

已知标定水分特征曲线ψm ～θ后，各测点标定系数a i 的确定应使得式（2－7－39）的SS 取最小值。如果各测点标定系数相互独立，则该点的标定系数a i 可根据δSS/δa i =0求得以下的计算式：

()

∑∑∑===?

?

?

???-=

i

i

J j ij

J j ij P p ij p ij i S a a 1

2

14

1

1?

θ?

（2-7-52）

N i ，?=,2,1

将已求得的待定常数和标定系数代人式（2－7－39）计算SS 值，此即各实测水分特征曲线标定后与标定水分特征曲线的方差值，记为SSB 。

在河南、河北、甘肃和北京等地取得的12个土样（砂性土）测其水分特征曲线（低吸力），每个水分特征曲线各有17～21个实测数据，总计有206个人，ψmij ～θij 测点，结果如图2－7－4所示。分别应用一次法和迭代法（5次）对这12个上样的水分特征曲线进行标定，结果如图2-7-5和图2-7-6。

二、田间非饱和土壤导水率的标定

土壤孔隙中的水流可视为粘性流，从理论上说可用Navier —Stokes 方程来描述。但是，由于土壤孔隙几何形状的复杂性，实际应用该方程是十分困难的。因此，描述多孔介质中流体流动的基本定律通常采用达西定律。将达西定律的表达式重写如下： 式中V ——渗透流速； grad ψ——水势梯度； k ——导水率。

直接用水动力学方法根据土壤中孔隙的大小和形状求得导水率k 的精确表达式是不可能的，一些研究者将土壤孔隙理想化，利用毛管模型，分析导水率和孔隙几何特性的关系。设毛管为圆形，半径为r ，水势梯度grad ψ沿管轴方向，此时管内流速v 与水势梯度grad ψ的关系由Poiseulle 方程给出：

?μ

grad r v 82

-=

(2-7-54)

式中μ——流体的粘滞系数。

比较式（2－7－53）和式（2－7－54），可知导水率：

μ

82

r k =

(2-7-55)

将此结果应用于非饱和土壤水分的流动，式（2－7－55）中的r 可理解为水气界面的曲率半径，考虑到土壤孔隙实际状况和理想的圆形毛管的差别，将式（2－7－55）应用于非饱和土壤水导水率时应乘以一个修正系数η，即

2

8r k μ

η=

(2-7-56)

当土壤中各处孔隙的几何形状和水分状态完全相似时，存在λ2

1/r 21=λ22/r 2

2。因此，当

各点的粘滞系数μ和修正系数η均相等时，不难得出：

2221

211

??λλλk

k k =?== (2-7-57)

式中k ——土壤含水率为θ（或饱和度Θ）时，土壤中各点导水率的平均值或标定值；

λ——特征长度的平均值或标定值。 如前所述记a i=λi /λ为i 点的标定系数，则

k a k i

i

?2

= (2-7-58)

以上分析表明，几何相似的土壤在同一含水率θ（或饱和度Θ）下，各点的导水率k i 是不同的，但除以标定系数的平方a 2i 多后，所得各点对应的标定导水率k 值是相等的。这样，在土壤完全相似的条件下，土壤中各点互不相同的导水率函数k i （θ）或k i （θ）按式（2-7-58）标定后，可得到完全相同的标定导水率函数k （θ）或k （θ）。不仅如此，对土壤中任一

点的基质势和导水率进行标定时，两者所用标定系数的值相同。由于田间各点土壤并非完全相似，因此非饱和土壤的导水率标定的结果不能得到单一的k ～θ曲线，同一测点基质势和导水率的标定系数也不会正好相等。为此，必须根据各点的实测资料，求出标定的导水率函数k （θ）和各测点的标定系数。

标定导水率函数k （θ）可用多种经验函数表示，例如指数函数和幂函数等。在用幂函数k=ηθβ

时。则有：

θβηln ln ln +=k

（2-7-59）

式中η,β——待定常数。

导水率标定的基本关系式（2－7－58）可改写为

i i a k k ln 2ln ln -=

（2-7-59’）

如田间观测点的编号为i=1，2，…，N ，每个点有实测不同饱和度θ时的导水率k ，各观测点测得的不同饱和度时的导水率的数目分别为J 1，J 2…，J N ，定义方差SS 为

()

2

11

ln ln 2?ln ∑∑==-+=N i J j ij

i ij i

k a k SS (2-7-60)

式中：lnk ij 是由实测的θij 用式（2－7－59）计算的结果。

在田间条件下应根据方差 SS 值最小的条件来确定标定导水率函数表达式中的待定常数和各点的标定系数，根据式（2－7－60）将SS 对a i 求导，得：

??????

?

?

?-=∑∑==i

J j J j ij

ij i J k k a i

i

2?ln ln exp 11

(2-7-61)

如标定导水率函数的表达式采用式（2－7－59），则将式（2－7－59）代入式（2－7－60），并将SS 分别对η和β求导，并使之等于0，可得待定常数η和β的计算式为

()()-???????-???? ??=∑∑∑∑∑∑∑∑∑=========N j N i J j ij i N i J j ij N i J j N i J j ij i ij ij i i i i J a k 111221111112ln ln ln ln ln exp θθθθη()()()???

????-???? ??∑∑∑∑∑∑∑∑∑=========N j N i J j ij i N i J j ij N i J j N i J j ij i ij i i i i

J a k 11122

111111

22

ln ln ln ln θθθ(

)()-

-?

??

? ??=

∑∑∑∑∑∑∑∑∑=========N

j N

i J j ij i N

i J j ij N i J j N i J j ij

i

ij i

i i

i

J a

k 11122

111111

2

ln ln ln ln θθθ

β()∑∑∑∑∑∑∑∑========-?

??

? ???

??

?

??N

j N

i J j ij i N

i J j ij N

j N

i J j i ij ij i i

i

i

J a k J 1112

2

111112ln ln ln ln θθθ

对导水率函数的标定也可用一次法或迭代法。在采用一次法时首先令a 1=a 2=…=a N =1（未标定），然后由式（2－7－62）和式（2－7－63）求得标定导水率函数中的待定常数η和β,按式（2－7－60）计算方差，记作SSA 。求得标定导水率函数k ～θ后，由式（2－7－61）可得到各点的标定系数a 1,a 2,…,a N 。将结果代入式（2－7－60），计算方差值，记作SSB 。

在采用迭代法标定导水率函数时，首先假定一组待定常数η，β.已知k ～θ后，由式

（2－7－61）求得各标定系数的计算值，记为a’i 。由于用迭代法计算时必须满足式（2－7－40a ）所示的约束条件，因此，对计算的a’i 。值需进行修正，即

∑='

'

=

N

i i

i i N

1

ααα （2-7-64）

将以上求得的各项待定常数和标定系数代入式（2－7－60），可计算SS 值记为SSB ，

此即各实测导水率函数标定后与标定导水率函数的方差。

根据第一次迭代求得的标定系数，由式（2－7－62）和式（2－7－63）求得以式（2－7－59）表示的标定导水率函数中的待定常数η和β,若令a 1=a 2=…=a N =1，由式（2－7－60）计算方差值SS ．此即标定导水率函数与各实测导水率函数的方差值，记为SSA 。

按上述两个步骤反复迭代，直至前后两次的SSB 值之差小于规定值为止。

同前，导水率k 和饱和度的关系也可采用一次法和选代法标定，两种方法标定结果相近

[67]

，表明一次法不仅计算简便，也有足够的精度。

第四节 土壤水分入渗的标定

由于土壤水分运动参数的空间变异性，在土壤水分人渗的过程中的任一时刻，田间各点的含水率分布和入渗量都是不同的。为了研究田间土壤水分入渗问题，首先必须选择适当的标定系数和待定常数，求得标定的土壤水分运动参数。同时，也要对时间、空间坐标和状态变量（如入渗量等）进行标定或无因次化，从而导出标定的土壤水分运动基本方程或标定的入渗方程。然后，对标定方程求解，得出所求问题的标定解，最后将田间具体点的标定系数回代到标定解中求得该点的具体人渗解，或土壤含水率或入渗量的空间分布。

以下举例简要介绍入惨问题的标定方法和应用。

一、入渗基本方程的标定及其数位解法

在初始含水率分布均匀，地下水埋藏很深，地表有薄层积水，水分向均质土壤垂直入渗时的定解问题表示为

z

k

z k z t C

??-??????????=???? （2-7-65）

0,0,≥==z t a ?? （2-7-65a ） 0,0,>==t z b ??

（2-7-65b ）

0,,>==t L z a ??

(2-7-65c)

式中 C ——容水度，等于d θ/d ψ

ψa 、ψb ——分别为初始和边界的土壤水基质势； L ——计算土层厚度； 其余符号同前。

基质势和导水率的标定关系同前：

2?,?a k k a ==??

容水度C 的标定关系为：

a C C

=? (2-7-66)

当已知ψ～θ后，C 可用下式计算：

θ

?θd d C

s

?1

?= (2-7-67)

式中 θs 一饱和含水率。

坐标Z 的标定关系取Z=Z/α加；时间t 的标定关系取为

3

?ta t

= (2-7-68) 利用上述变量的标定关系，则垂直入渗定解问题式（2-7-65）转化为下列标定定解问题：

z

k z k z t C ??-??????????=????????? （2-7-69a ）

0?,0?,??≥==z t a ??

（2-7-69b ） 0?,0?,??>==t z b ??

（2-7-69c ）

0?,??,??>==t L z a ??

(2-7-69d)

其中，ψa =αψa , ψb =αψb ,。L=αL 。在已知ψ～θ、k ～θ、C ～θ，以ψa 、ψb 值（假设田间并点θa ，θb 均相同），可用数值方法求解式（2－7－69），得出标定ψ（z 、t ）。这一标定解对田间各点均是适用的，由ψ～θ即可转化为θ（z 、t ）。

二、入渗公式的标定和入渗量的空间分布

（一）入渗参数和入渗方程的标定

在前面第五章中介绍了几种描述土壤水入渗过程的入渗方程，以下仅以PhiliP 入渗方程为例，说明入渗方程的标定的原理与方法。 Philip 入渗方程为 At St I +=2

1

(2-7-70)

A t S i +???

? ??=212

(2-7-71)

式中 I ——累积入渗量（cm ）；

i ——入渗率（cm ／min ）； t ——时间（min ）；

S ——入渗参数（cm ／min 1/2）

A ——稳渗率（cm ／min ）。

()t k d t Z I a b

a

+=?θθθθ,

(27-72)

式中 θa 、θb ——分别表示含水率的初始值和地表值；

k a ——相应于含水率民的导水率。

将philip 垂直入渗级数解前两项Z ＝η1（θ）t 1/2＋η2（θ）t 代入式（2-7-72），即可得式（2-7-70），其中： ()θθηθθd S b

a

?=1

(2-7-73)

()a k d A b

a

+=?θθηθθ2

(2-7-74)

由于η1（θ）即是同样定解条件下的水平渗吸解，并可表示为

()2

11-=ct θη

式中x ——距离坐标。

应用标定关系x=x/α和t=t/α3则式（2－7－73）可改写为

θθθ

d t z

a

S b

a

?=??2

1 今 S=∫θ

b

θ

a xtd θ

，S 称为 S 的标定值，在应用试验资料求标定系数 a 时，一般取 S 为各

点S i 值的平均值S ，于是

2

??

?

??=S S a

(2-7-75)

由于地表入渗率i 是在Z=0处土壤水分通量，故

???

? ?

?-??-==10

Z b Z

k i ?

k b 为相应于含水率θb 的导水率，当t 足够大时，δψ/Δz=0，I=k b ，由式（2－7－71）可知，当t 足够大时，I=A ，因此，可近似为A=k b 。

由于 k b = k b α2， A= k b α2，A= k b ，A 称为A 的标定值，若取A 为各点试验值A i 的平均值，则

2

1??

? ??=A A a

(2-7-76)

将Z=Z/α，t=t/α3，k=k/α2代入式（2－7－22）可得：

()Ia t k d t z I a b

a

=+=???,??θθθθ

(2-7-77)

将式（2－7－75）、式（2－7－76）和式（2－7－77）代人式（2－7－70），则可得到标定的入渗公式：

t A t S I

???2

1

+= (2-7-78)

在求得入沙参数的平均值S 和A 时，即可求得标定的入渗曲线I 一t ，如已知I 一t 和M 某

一点的标定系数a ，则可自I=Ia 和t=ta 2算出该点的入渗曲线I ～t 。

（二）由试验资料计算标定系数

根据入惨试验结果，首先计算各点的入渗参数S i 和A i 值，然后再根据式（2－7－75）和式（2 –7-76）计算．各点的标定系数a i 值，记为 a si 和 a Ai 。由于a si 和 a Ai 是基于相似理论而求得的．这与实际情况并不完全相符，因而a si 和 a Ai 并不相等。标定系数的取值可以采用a si 和 a Ai 的调和平均值a hi ，即

Ai

Si Ai

Si hi i a a a a a a +=

=2

(2-7-79)

或按最小二乘原理取值，记为a opti ，即每一点标定系数的取值应使该点实测的入渗曲线I ～t 经标定后与标定曲线I ～t 之间的方差ss 最小，方差为

()

2

11

?∑∑==-=N

i J j ij

i ij i

I a I SS (2-7-80)

式中 N ——入渗试验的点数；

J i ——第i 个试验点的实测的I ～t 关系中资料的“对”数；

I ij ——相应t ij 由式（2－7－78）计算所得的值。

若认为各点的标定系数是相互独立的，则要求任一点i 的标定系数a i 应满足下式方差SSI 最小：

()

∑=-=i

J j ij

i ij I a I SSI 1

2

? (2-7-81)

根据入渗试验资料，用逐次逼近法可求得SSI 值最小时的标定系数a opti 。

例如，图2－7－7是在山东茬平进行的入渗试验资料点绘出I ～t 关系[66]，经按最小二乘法原理，使SSI 最小，求得标定系数a ，其标定入渗曲线I ～t 如图2－7－7b 所示，经标定后各次观测点与标定曲线I ～t 相当吻合。