J. Jpn. Soc. Soil Phys.
土壌の物理性No. 155, p.35 ∼ 44 (2023)

団粒構造が発達した土壌の水分特性曲線の回帰手法の改良— dual-van Genuchten モデルのパラメータ決定の自動化に向けた取り組み —関 勝寿 1 · 岩田幸良 2 · 柳井洋介 3 · 亀山幸司 2

Improved regression of the water retention curve for highly aggregated soils:
Toward the automation of parameter determination for dual-van Genuchten models
Katsutoshi SEKI1 , Yukiyoshi IWATA2 ,Yosuke YANAI3 and Koji KAMEYAMA2

1. はじめに
Abstract: The van Genuchten (VG) model is widely used
to represent soil water retention curves. Durner’s dual-VG
model, which combines VG functions linearly, is suitablefor expressing retention curves of highly aggregated soils
like Andisols, which are commonly found in Japan’s arable
land. However, we encountered difﬁculties in determining
appropriate model parameters for the dual-VG model using the regression method. To overcome this difﬁculty, we
improved the algorithm for nonlinear regression using thedual-VG model. This was accomplished by separating thewater retention curve into subcurves that represent regionsof high and low water content. Each subcurve was individually ﬁtted with VG models, after which the entire curve
was ﬁtted with the dual-VG curve. The new algorithm was
implemented in SWRC Fit, which is the web-based application to determine model parameters automatically andits accuracy was validated using measured data from agricultural soils throughout Japan. The results demonstrated
that the proposed method generated retention curves witha goodness of ﬁt comparable to the global solution, while
maintaining signiﬁcantly low computational costs. The developed method is expected to contribute to improve theresults of numerical simulation with well-aggregated soil
in arable land by accurately estimating water holding capacity and hydraulic conductivity.
Key Words : Water retention curve，nonlinear regression，
dual-VG model

水分特性曲線は土壌の保水性を評価するための基本的な物理量である．また，この曲線は間隙形状を反映していることから，この曲線から不飽和透水係数を推定できるため，土壌中の水分移動に関する数値シミュレーショ
ンを実施する上で重要なパラメータである．土壌水分動態の解析は，適切な灌水スケジューリングの検討や硝酸態窒素等の肥料成分の動態解析，農地や森林を含む広域な水 · 物質循環の評価等にとって重要である．そのため，土壌の水分特性曲線パラメータのデータベースを整備することが求められているが，土壌の化学性のデータベースに比べ（農研機構, 2021），物理性のデータベースは特に日本国内においてはほとんど整備されていない現状にある．
水分特性曲線は一般的には土壌を実験室内に持ち込んで測定する．飽和した試料を段階的に排水させ，所定のマトリックポテンシャル h（cm）に調整した試料の体積
含水率 θ を測定することで，h と θ の関係に関する複数点のデータを得る．これらの離散的なデータから両者の連続的な関係を得るために，いくつかのモデルが提案されている．その中で最も広く使われている代表的なもの
が以下の VG モデル（van Genuchten, 1980）である（不飽和土壌の圧力水頭はマイナスの値だが，表記の利便性から −1 を掛けて正値として扱っていることに注意）．

θ = θr + (θs − θr )[1 + (α h)n ] n −1
1

1 Natural

Science Laboratory, Toyo University, 5-28-20 Hakusan,
Bunkyo-ku, Tokyo 112-8606, Japan. Corresponding author. 関勝寿, 東洋大学自然科学研究室
2 Institute for Rural Engineering, National Agriculture and Food Research

Organization, Kannondai 2-1-6, Tsukuba, Ibaraki 305-8609, Japan.
3 Institute of Vegetable and Floriculture Science, National Agriculture andFood Research Organization, Kannondai 3-1-1, Tsukuba, Ibaraki 3058519, Japan.
2023 年 6 月 19 日受稿 2023 年 10 月 3 日受理(1)

ここで，θs は飽和体積含水率，θr は残留体積含水率，α（cm−1 ）と n は土壌構造により決定されるフィッティン
．
グパラメータである（α > 0，n > 1）

VG モデルは間隙径に対する間隙量のピーク（間隙径36

土壌の物理性第 155 号

(2023)

分布）が 1 つであるような間隙構造を仮定しており，団

適切な初期値を与えるためのアルゴリズムの一つとし粒構造が発達した黒ボク土のように間隙径分布のピークて，SWRC Fit では単純な式のフィッティング結果をよが 2 つ（団粒間間隙によるピークと団粒内間隙によるり複雑な式のフィッティングにおける初期値とするこピーク）ある土壌では，VG モデルを 2 つ足し合わせるとで，段階的により複雑な式のフィッティングをするこ

dual-VG モデル（Durner, 1994) がより適しているとされている（坂井 · 取出, 2007）．dual-VG モデルは 2 つのVG モデルを重み係数 w で線型結合したとを基本設計としている（Seki, 2007; 関, 2019）．その基本設計に従って， VG モデルのパラメータと w = 0.5
を初期値として dual-VG-CH モデルのフィッティングをし，その次に dual-VG-CH のパラメータを初期値として

θ=
[
]
1
1
(2)
n1 n1 −1n2 n2 −1θs w[1 + (α1 h) ]+ (1 − w)[1 + (α2 h) ]という式で表される．ここで，θs ，θr ，w，α1 ，n1 ，α2 ，n2はフィッティングパラメータであり，0 < w < 1，n1 > 1，

n2 > 1 である．
dual-VG-CH モデル（Seki et al., 2022）は，dual-VG モ
デルにおいて α = α1 = α2 とした

dual-VG モデルのフィッティングをするという方法（CH法）を SWRC Fit では従来採用していた．しかしながらこの方法では dual-VG モデルで十分に良い回帰が得られないことがあった（結果は 3 節に記載）．その原因は 2
つあり，1 つはパラメータ値の範囲に制限を与えていなかったため n1 あるいは n2 が不自然に大きい値になってしまうこと，もう 1 つは初期値として dual-VG-CH モデ
ルのパラメータを使っているため α1 = α2 となる初期値近傍での局所最適解を探索している可能性があることが考えられた．この中で 1 つ目の問題に関しては，SWRC

θ=
[
]
1
1
(3)
n1 n1 −1n2 n2 −1θs w[1 + (α h) ]+ (1 − w)[1 + (α h) ]Fit にパラメータの上限を設定する機能を実装する，軽微な修正で対応できると考えられる．一方，2 つ目の問題に関しては，初期値を設定するアルゴリズムの改良と，その有効性の検証が必要となる．

という式である．ここで，CH は Common Head の略でそこで本研究では，水分特性曲線から dual-VG モデル

あり，パラメータ α1 と α2 を共通化することを意味してのパラメータを決定するアルゴリズムについて，従来のいる．パラメータは θs ，w，α ，n1 ，n2 の 5 個であり，dualVG モデルよりも 1 つ少なくなる．Seki et al.（2022）は，
黒ボク土のように団粒化が進んだ土壌では dual-VG モデ
ルが適しているが，砂丘砂土壌では dual-VG-CH モデル手法である「CH 法」に加えて，新しく「分割法」を開発し，日本の農地土壌および UNSODA データベースの水分特性データを用いて，改良したモデルによる精度向上の効果を検証した．

を使うことでパラメータを減らしながらも広い圧力範囲の水分特性曲線を的確にあらわせることを例示した．さ2． 方法

らに Seki et al.（2023）では多くの UNSODA データベース（Nemes et al., 2015）の土壌において dual-VG-CH モ
デルの有効性が確かめられている．

2.1 試料測定された水分特性のデータから VG モデルや dual-

北海道から沖縄まで，13 地域の試験圃場から土壌を

VG モデルのパラメータを決めるためのツールの一つとして，RETC（van Genuchten et al., 1991）がよく使われているが，RETC ではパラメータの初期値の入力が必要である．特に dual-VG モデルはパラメータが多いため，採取した．採取した圃場の位置を Fig. 1 に示す．島尻

適切な解を得るためには初期値を適切に与えることが重においては，後述（3 章）のように水分特性曲線が二峰要であり，そのためには，初期値を土性ごとの平均的な性を示すものが多く，供試土壌は団粒構造が発達したも値から与える方法が有効である．しかしそのためには，のが多い．

マージを除く土壌は全てアスパラガスの圃場から採取した．供試土壌には一部真砂土のような砂質土等も含まれるが，黒ボク土，褐色森林土，低地水田土，灰色低地土

ユーザーがいくつかの土壌についてモデルパラメータを団粒構造が発達した赤黄色土として，宮古島のサトウ

すでに持っていることや，適切な初期値を選択するため

キビ畑の作土層からサンゴの化石である琉球石灰岩を母に経験と勘を要するなど，水分特性曲線に精通していな材とする島尻マージを採取した．アスパラガス圃場につい研究者や生産者にはハードルが高いという問題があっいては，高さ 40 ∼ 60 cm 程度の畝を立ててその上にアた．そこで Seki（2007）が開発した SWRC Fit では，Web

スパラガスを定植する枠板式高畝栽培（池内 · 佃, 2020）

ブラウザに水分特性のデータを入れることで，適切な初の圃場において，畝部ならびに畝の下の土壌を採取した．

期値を与える等のユーザーの経験を要する作業を必要となお，塩尻と諫早については一般的なアスパラガスの栽

せずに，VG モデル，dual-VG モデル，dual-VG-CH モデ

培方法である平畝栽培の圃場の各土層からも土壌を採取

ルなどの水分特性モデルのパラメータを自動的に決める

した．畝部の採取深さは 10，30，50 cm の 3 つを基本とことができるように工夫をした．

した．

論文：団粒構造が発達した土壌の水分特性曲線の回帰手法の改良

37

Fig. 1 試料採取地．
Soil sampling sites.

採取深度毎に 100 cm3 円筒サンプラーを用いて未攪乱お，アスパラガスは多量の灌水を実施することから塩類土試料を 3 つずつ採取し，水分特性曲線の試験に供試し集積等の問題は生じていないと考えられたことと，試料た．13 地域のそれぞれから複数の地点と深度におけるは蒸留水で毛管飽和したことから，土壌の水ポテンシャ

試料を採取したことから，水分特性曲線の試験に供試し

ルはマトリックポテンシャルと同等と扱い，測定されたた試料数は合計で 98 となった．また，同じ深さから攪

圧力値（負圧）を水柱高に換算して対数をとることで pF

乱土を採取し，粒度試験等に供した．

に換算した．

比較のために，不飽和土壌水文データベース UNSODAの実験室の水分特性（脱水過程）の中から，8 点以上の

h と θ 関係の測定点がある試料 565 個のデータを解析した．

攪乱土を用い，ピペット法で粒度試験を実施し，国際法（Murano et al., 2015）で土性を分類した．

2.3 非線形回帰

1.8 ∼ pF 3.0 を加圧板法で水分調整した（寺沢, 1972）．
各 pF 段階において試料の重量を測定し，試験終了後にSWRC Fit では Python ライブラリの unsatﬁt（Seki etal., 2023）を利用して回帰をしている．本研究では unsatﬁt によって一括して計算を実行した．unsatﬁt では
Python の scipy.optimize.least squares パッケージを利用して，信頼領域 Reﬂective 法による非線形回帰によって

炉乾燥して乾燥重量を求めることで重量法により各マト

パラメータの最適化をしている．これはニュートン法のリックポテンシャルに対応する体積含水率と試料の乾燥ようにパラメータの初期推定値から目的関数の勾配に密度を計算した．

よって目的関数を最小化するパラメータを探索する手法

2.2 測定方法
未攪乱土試料を毛管飽和し，pF 1.5 までを砂柱法，pF

攪乱土を 50 cm3 円筒容器に詰めて毛管飽和した後，であり，一定の範囲内からパラメータを決定することが加圧板法で pF 3.5 に水分調整した試料について鏡面冷できる．ここで θ の実測値と推定値の RMSE（二乗平却露点計（Meter 社，WP4T もしくは WP4C）で試料の均平方根誤差）を目的関数とした．なお，パラメータに水ポテンシャルを測定した．試料は鏡面冷却露点計用のまったく制限を与えないと，n1 または n2 が 100 を超え

ステンレス容器に移し替え，各試料につき 3 反復で試験るような不自然な値へと発散して不自然な形のグラフがを実施した．水ポテンシャル測定後に試料の重量を測定得られることがあり，n1 と n2 の上限値を設定することし，1 ∼ 2 時間程度風乾させた後，再度水ポテンシャルが望ましいことが予備的な研究により明らかになった．

を測定し，試料の重量を測定した．この操作を水ポテン

不自然な曲線形とならないような上限値を検討した結シャルが pF 4.2 に相当する値になるまで繰り返し実施し果，本研究では n1 と n2 の最大値を 8 に設定することとた．試験終了後，試料を炉乾燥し，重量法で含水比を求した．

め，乾燥密度を掛けることで体積含水率を計算した．な

非線形回帰の方法については，「はじめに」に記したよ土壌の物理性38

第 155 号

(2023)となることから，高水分領域において θ2 = θm と近似す

VG1 (1/α=29.4, n=1.72)

れば θ1 と一致し，低水分領域において θ1 = θm と近似

0.4

すれば θ2 と一致する．そのことから，（θs ，w，α1 ，n1 ，α2 ，n2 ）は dual-VG モデルの良い初期パラメータとな

θ

θm

0.2

る．CH 法では α1 = α2 を初期値とするため，α1 と α2が近い値で収束しやすいが，分割法では α1 と α2 がそ

VG2 (1/ α=1629, n=1.46)

れぞれ高水分領域と低水分領域から選ばれるため，それぞれまったく異なる初期値からのパラメータの推定とな0る．ここで，Fig. 2 において 1/α1 と 1/α2 は高水分領域

1

10

100
1000
h(cm)

10000

と低水分領域それぞれの変曲点に近いサクションに相当している．
この他，上記の CH 法と分割法によって得られた近似

Fig. 2 大仙黒ボク土の水分特性曲線における dual-VG モデ
ルの初期パラメータ決定のための分割法の説明．

Illustration of the separation method for determining initial parameter of dual-VG model for Daisen Andisol.

曲線と測定データの RMSE を計算し，値が小さい方のパラメータを dual-VG モデルのパラメータとして採用する手法を本論文では「2 点法」と名付ける．

2.4 モデルの比較
VG，dual-VG-CH， dual-VG の 3 つの水分特性モデ

うに，従来の手法である「CH 法」に加えて，新しく「分ルで非線形回帰をした．ここで一般的にモデルはパラ割法」を開発した．CH 法は，dual-VG-CH のパラメータメータの数が多いほどデータへの適合度が高くなる傾向

を初期値とする方法である．すなわち，式（3）で回帰しがあるが，適合度が同程度のモデルの中からはパラメー

たパラメータを初期値とし，dual-VG モデルで再度回帰して最適なパラメータを得る．
分割法は，水分特性全体を高水分領域と低水分領域に分割し，水分特性曲線全体で回帰をする方法であり，6 点以上の測定点を必要とする．水分特性測定データにおける θ の最大値 θmax と最小値 θmin の平均 θmed よりも小さい最大の θ の測定値を θm として，θ > θm のデータを
高水分領域，θ ≦ θm のデータを低水分領域と分割する．
ただし，高水分領域または低水分領域の測定点数が 3 点に満たない場合には，3 点となるように θm を変更する．
高水分領域を θs = θmax ，θr = θm として VG 式で回帰した曲線を VG1，低水分領域を θs = θm ，θr = 0 として

VG 式で回帰した曲線を VG2 とする．一例として，大仙黒ボク土の水分特性曲線についてこの手法により θmと VG1，VG2 を決定した曲線を Fig. 2 に示す．VG1 とVG2 のパラメータをそれぞれ（α1 ，n1 ），（α2 ，n2 ）とす
ると，VG1 と VG2 の曲線はそれぞれ
1 −1

θ1 = (θmax − θm )[1 + (α1 h)n1 ] n1

+ θm

(4)

タが少ないモデルを選ぶのが良いとされる（Gauch Jr.,

2002）．そこで，本論文ではパラメータを減らして良好な回帰結果を得ることを目的として， dual-VG（-CH）モ
デルにおいては，あらかじめ θr = 0 とした式を用いた．
θr = 0 とした理由は，① θr はフィッティングパラメータに過ぎず，線型結合モデルにおいて θr の値は透水性曲線の形に影響を与えないと Durner（1994）が言及している
こと，② Seki et al.（2022） において黒ボク土の水分特性曲線を θr = 0 とした dual-VG モデルで精度良く近似できることが示されたことから，θr をパラメータから除外してパラメータ数を減らすことで，モデルの複雑さが低下してより良いモデルになると判断したためである．

VG，dual-VG-CH， dual-VG の 3 つのモデルの中から，AIC（赤池情報量規準）によって最良のモデルを判定した．ここで AIC は

AIC = −2 log L + 2k

(7)

（L は最大尤度．k は自由パラメータの数．log は自然対数）で定義される．一般的にモデルはパラメータの数が多いほど適合度が高くなる傾向があるが，パラメータが

1 −1

θ2 = θm [1 + (α2 h)n2 ] n2

(5)となり，θs = θmax ，w = 1 − θm /θmax とすると dual-VG

多いとそれを決定するために多くのデータが必要となるなど，デメリットもある．統計モデルの指標である AIC
が最小のモデルを選ぶことで，モデルの複雑さ（パラモデル（式（2））は，メータの数）とデータの適合度とのバランスが取れた良[]1
1
n1 n1 −1n2 n2 −1
θ = θs w[1 + (α1 h) ]+ (1 − w)[1 + α2 h) ]いモデルの選定が可能になる．ここで標本の誤差が正規1
n1 n1 −1

= (θmax − θm )[1 + (α1 h) ]1n2 n2 −1

+ θm [1 + (α2 h) ]= θ1 + θ2 − θm(6)分布であると仮定すると

AIC = 2N log(σ ) + 2k
（σ は標準誤差．N は標本サイズ）となる．

(8)

論文：団粒構造が発達した土壌の水分特性曲線の回帰手法の改良
ただし，AIC は標本サイズ（データ数）無限を仮定しているため，標本サイズが小さい場合は AIC 最小化による

39

ときには，
（θs ，w，n1 ，n2 ）には Table 1 の初期値を設定することで大域的な探索をした．

モデル決定はパラメータ数を過大に見積もってしまうこ

3. 結果とになる．この対策として，標本サイズが有限の場合の次のような修正 AIC（AICc）が用いられる（Cavanaugh,

供試土壌の土性の測定結果を Fig. 3 に示す．Fig. 3 の1997）．

分類は日本の包括的土壌分類（小原ら, 2011）の大群に

AICc = 2n log(σ ) +2kN
N −k−12k(k + 1)= AIC +
N −k−1

よるものである．粘土含量が高い赤黄色土，比較的高め(9)n が十分に大きければ AICc は AIC に近似される．
VG，dual-VG-CH，dual-VG の 3 つのモデルで非線形回帰をして，修正 AIC が最小となるモデルがそのデー
タに対する最良のモデルであると判定をした．ここで，dual-VG モデルでは k = 6 であるため，修正 AIC が定義される N − k − 1 > 0 が成立するためには 8 個以上の測定点が必要となる．

2.5 大域探索の黒ボク土と褐色森林土，砂含量が高い未熟土，多様な土性の低地土と，多様な土性の土壌を供試したことがわかる．
水分特性曲線については，n 値に上限値を設け，分割法または 2 点法を用いる本論文で提案する手法を適用することで全体的に非常に適合度が高い近似曲線が得られた．それぞれの土壌で修正 AIC が最小のモデルを最適なモデルとして選択したところ，dual-VG モデルが最適である土壌が 72 個となり，本研究で供試した日本の農地
土壌の 73 % であった（Table 2）．UNSODA のデータで35 % であったことと比較すると，本研究で試験対象とし非線形回帰においては，パラメータの初期値から目的関数を最小化するために目的関数が低下する方向へとパラメータを変化させるため，パラメータ空間における目
的関数すなわち RMSE グラフの形状によっては，パラ

Table 1 大域探索初期パラメータ設定．mi = 1 − 1/ni ．
Initial parameter sets for the global search．mi = 1 − 1/ni .

メータが不自然な局所最適解に収束する可能性がある．

パラメータ

下限

上限

初期値の数

そこで，本論文で提案したアルゴリズム（2 点法）が大域

θ s/θmax
w1− log10 (α1 )m1− log10 (α2 )
m2

0.95
0.05
−1
0.01
−1
0.01

1.2
0.95
3
0.875
5
0.875

3
8
5
5
8
5

探索による最適解と比べて誤差が無視できないほどに大きい局所最適解に収束してないことを確かめるために，大域探索を実施し，その結果を 2 点法により得られた結果と比較した．大域探索は，Table 1 に示すパラメータの範囲内にそれぞれの個数のパラメータを初期値として設定した．すなわち，下限から上限の区間を p + 1 等分して，上限と下限以外の p 通りの分割点を初期値として設定した（初期値の総数は 24,000 個になる）．ここで，パラメータの初期値を計算する際には Table 1 に記載されているように元のパラメータを変換したパラメータ，たとえば ni から mi = 1 − 1/ni と変換された mi のようなパラメータを使っている．

黒ボク土
低地土
未熟土
褐色森林土
赤黄色土
Heavy
Clay

大域解の計算には時間がかかるため，unsatﬁt に実装されている多くの初期値からの解の探索を効率化する機能

45を利用した．すなわち，まずは緩やかな終了条件で速やかに多くの初期値からの回帰をして．RMSE が最小となる解から厳しい終了条件で回帰を続けることにより，計算時間を短縮している（Seki et al., 2023）．
また，パラメータ空間における RMSE の感度を確認することとした．特に本研究では，CH 法により α1 = α2を初期値とする場合と，分割法により α1 と α2 の初期値がそれぞれ高水分領域と低水分領域から選ばれる場合と

Sandy
Clay
25 Sandy
Loamy
Sand
15 Clay Loam
Sandy
Loam

5
85

Sand

45
Light
Clay
Clay Loam

Silty Clay Loam

Loam

Silty Loam

55

65

Sand / %の比較をするため，dual-VG モデルの 6 個のパラメータの中で，α1 と α2 に着目して，この 2 つのパラメータを変化させたときに，それ以外のパラメータ（θs ，w，n1 ，n2 ）を最適化したときの RMSE の等高線を描いた．その

Silty
Clay

Fig. 3 試料の土性三角座標．
Soil texture triangle.

75
85

土壌の物理性40

第 155 号

Table 2 モデルごとの修正 AIC 最小となった試料の数．
Numbers of soils with the least corrected AIC for respectivemodels.
土壌
日本農地

VG dual-VG-CH dual-VG 合計
黒ボク土
低地土
褐色森林土
未熟土
赤黄色土

1
9
5
0
1
16

0
1
2
7
0
10

34
23
12
2
1
72

35
33
19
9
2
98

274

93

198

565

(2023)

まで低下して θ がほぼ一定となる曲線となるが，dual-

VG-CH モデル（RMSE 0.00061，修正 AIC −141）では，
h = 100 cm 以降のなだらかな θ の低下を 2 つ目のサブ関数で表現できている．Dual-VG モデルでは dual-VG-CH
モデルよりも RMSE がわずかに低下するが（0.00058），
修正 AIC は大きくなる（−131）．したがって，Fig. 5 のデータに対しては dual-VG-CH モデルが最良のモデルであると判断できる．

合が多いことがわかる．特に，黒ボク土では 35 個の中Fig. 6 上の香川低地水田土（Clay loam）では，VG モデル（RMSE 0.019，修正 AIC −72）と dual-VG-CH モ
デル（RMSE 0.019，修正 AIC −65）に比べて，dual-VG
モデル（RMSE 0.0025，修正 AIC −99）では RMSE と
修正 AIC が最小となった．ここで，dual-VG モデルの近似において，従来の手法（CH 法）では適合度が良くなかった（RMSE 0.019，修正 AIC −54）ものの，本研究ので 34 個が dual-VG モデルが修正 AIC 最小となり，低地

手法（分割法）に改良したことによって良い適合が得ら土（低地水田土，灰色低地土）と褐色森林土でも dual-VG

れた．これはパラメータの初期値を変えた効果であり，モデルが修正 AIC 最小となることが多かった．一方，未この効果については次節で詳しく検証する．

合計

UNSODAた農地土壌においては，特に dual-VG モデルが有効な場

熟土（真砂土）では dual-VG-CH モデルが修正 AIC 最

ここで，n に上限を設ける効果を見るために，n に上小となることが多かった．Dual-VG モデルが最適である

限なしとした場合の CH 法の結果を示すと，長野県中野

72 個の日本の農地土壌では，RMSE は 0.0037 ± 0.0016，

市黒ボク土（Heavy clay）（Fig. 6 下）のように，ある h決定係数は 0.99887 ± 0.000934 であり，いずれもきわめの値で θ が直角に急降下する不自然な形の水分特性曲線て良い近似曲線となった．

が得られるものがあった．これはパラメータに制限を課

Fig. 4 は VG モデルが修正 AIC 最小となった香川低地
水田土（Sandy clay loam）である．RMSE は VG（0.0057）> dual-VG-CH モデル（0.0044）> dual-VG（0.0030）とdual-VG モデルが最小であるが，修正 AIC は dual-VG（−94.8）> dual-VG-CH（−97.54）> VG（−99.0）と VG
モデルが最小となった．すなわち，Fig. 4 のデータに対しては VG モデルが最良のモデルであると判断できる．
Fig. 5 は dual-VG-CH モデルが修正 AIC 最小となった
久留米真砂土（Sand）の水分特性曲線である．VG モデル（RMSE 0.0050，修正 AIC -102）では速やかに θr = 0.02

していなかったため n1 が 105 以上の極めて大きな値とFig. 4 香川低地水田土（Sandy clay loam）における VG モ
デルと dual-VG モデルの水分特性曲線近似比較．

Fig. 5 久留米真砂土（Sand）における dual-VG-CH モデルと VG モデルの水分特性曲線近似比較．
Fitted curves of dual-VG-CH and VG models for Kurume decomposed granite (Sand).

Fitted curves of VG and dual-VG models for Kagawa lowland
paddy ﬁeld soil (Sandy clay loam).

なったためである．本研究の手法では n1 と n2 に上限値を設定したことにより，適切なパラメータの範囲内で解が探索されるようになった．なお，長野県中野市黒ボク
土では CH 法（上限なし）の RMSE は 0.0108，本研究の手法では RMSE は 0.0021 である．また，dual-VG-CH
モデルが修正 AIC 最小となった 72 個の試料の中で，nに上限を設けない CH 法で n1 または n2 が 8 を超える値が得られた試料は 5 個であり，その中の 4 個は n1 > 105となった．

論文：団粒構造が発達した土壌の水分特性曲線の回帰手法の改良

4. 考察

41

近にデータが多い．すなわち，CH 法では α1 = α2 を回帰の初期値とすることから，初期値に近い局所最適解に4.1 CH 法と分割法の比較

収束し，α1 と α2 が離れた値になっているところに存在

Fig. 7 に，修正 AIC から dual-VG モデルが最適である

する解を探索できていない可能性がある．

と判断した 72 試料（Table 2）について，CH 法と分割法

これを検証するために感度解析を実施した．大仙黒ボによって得られた近似曲線の RMSE を比較する．なお，ク土（Heavy clay）では Fig. 9 のように CH 法と比べてこの比較ではいずれの手法でも n1 と n2 の上限値を設分割法が明らかに良い近似が得られている．そのパラ定している．1 : 1 直線上に多くのプロットがあるが，分メータと RMSE を Table 3 に示す．Fig. 10 にこの土壌割法によって大きく RMSE が低下した試料も多かった．

の感度解析の結果を示す．上の図に様々な組み合わせでこのことから，本研究の分割法の有効性が確認された．

α1 と α2 を固定し，他のパラメータを推定したときの最適化されたパラメータによる RMSE の等高線を示した．
図の三角のプロットが CH 法で得られた水分特性曲線の（α1 , α2 ），丸のプロットが分割法で得られた曲線の（α1 ,α2 ）である．下の図に，これら二つのプロットを結んだ

一方，CH 法によって RMSE が低下した試料もあり，CH法と分割法を両方計算する計算コストもそれほど高くはないことから, 2 つの方法によって RMSE の小さい方を選択する 2 点法もまた実用的な手法である．

Fig. 8 で CH 法と分割法の log10 (α1 /α2 ) を比較すると，CH 法では α1 = α2 すなわち log10 (α1 /α2 ) = 0 の付

Fig. 7 CH 法と分割法の RMSE 比較．
Comparison of RMSE for CH method and separation method.

Fig. 6 香川低地水田土（Clay loam）における従来の手法と本研究の手法における dual-VG モデルの水分特性曲線比較と，長野県中野市黒ボク土（Heavy clay）における CH 法（nに上限なし）と本研究の手法における dual-VG モデルの水分特性曲線比較．

Fitted curves of dual-VG model with the conventional and improved methods of ﬁtting for Kagawa lowland paddy ﬁeld soil
(Clay loam) and with the CH method (without upper limit ofn) and improved methods of ﬁtting for Nakano Andisol Heavyclay).

Fig. 8 CH 法と分割法の log10 |α1 /α2 | 比較．
Comparison of log10 |α1 /α2 | for CH method and separation
method.

土壌の物理性

42

第 155 号

(2023)

Table 3 大仙黒ボク土の CH 法と分割法による dual-VG モデルパラメータと RMSE．
Parameters and RMSE of CH method and separation method for Daisen Andisol.

CH 法
分割法

θs

w

α1

n1

α2

n2

RMSE

0.705
0.648

0.779
0.363

1.40
0.128

1.11
2.28

1.40
0.000846

1.11
1.25

0.02190.00309

直線上の RMSE の変化を示している．下の図から，CH
法の解からスタートして RMSE が最小になるように最適解を探索した場合，分割法の解の方向に向かう場合には RMSE がわずかに上昇している．従って，非線形回帰
で RMSE が最小になるようにする最適解の決定方法では，CH 法の初期値から分割法の解には到達することができないことが分かる．

4.2 大域法との比較
Fig. 11 で，2 点法（CH 法と分割法による方法）による解と大域法による解の RMSE を比較した．このように，日本の農地土壌においてはほぼすべての土壌で 2 点法によって大域法とほぼ等しい RMSE の解を得ること
Fig. 9 大仙黒ボク土（Heavy clay）の CH 法と分割法による
dual-VG モデルの近似曲線．
Fitted curves of dual-VG model with CH method and separation
method for Daisen Andisol (Heavy clay).

ができた．最も RMSE の差が大きかったものでもわず
か 0.0011 であり，Fig. 12 のように実用上は適合度がほぼ同じであると判断できた．UNSODA で dual-VG モデ
ルが適していた 198 個の試料について同様の比較をした
ところ，大域法と 2 点法の RMSE の差が最も大きかったものでもその差はわずか 0.0022 であり，その曲線を比較すると，Fig. 13 のように実用上は 2 点法で十分な適合度の曲線が得られた（Fig. 11 の 1 : 1 線上に存在する多くの UNSODA のプロットは日本の農地土壌のプロット

Fig. 10 大仙黒ボク土の dual-VG モデル近似における α1 と
α2 を固定して他のパラメータを非線形回帰で推定した時のRMSE の応答と CH 法と分割法により得られた α1 ，α2 ．
Response of RMSE when α1 and α2 are ﬁxed and other parameters are ﬁtted in the approximation of dual-VG model for
Daisen Andisol.

Fig. 11 日本の農地土壌と UNSODA 土壌における 2 点法と大域探索による解の RMSE 比較．
Comparison of RMSE for the 2-point method and the resultof the global search method for Japanese agricultural soils and
UNSODA soils.

論文：団粒構造が発達した土壌の水分特性曲線の回帰手法の改良

43

5. おわりに
日本の農地土壌，特に黒ボク土の水分特性曲線を近似するのに適している dual-VG モデルによって迅速に適合度の高い回帰曲線を得るアルゴリズムを開発し，SWRC

Fit に実装した．そのアルゴリズムとは，水分特性曲線
を高水分領域と低水分領域に分割し，それぞれを VG モ
デルで近似することで，dual-VG モデルの初期パラメー
タを得るというものである．本研究で開発したアルゴリ
ズムは，検証されたほぼすべての土壌において従来の手法と比較して極めて良い近似精度を示し，従来の手法とFig. 12 香川低地水田土（Clay loam）の 2 点法と大域探索による dual-VG モデルの近似曲線．
Fitted curves of dual-VG model with the 2-point method andthe global search method for Kagawa lowland paddy ﬁeld soil
(Clay loam).

あわせた 2 点法を使うことによって，大域解とほぼ等しい適合度の曲線が得られることが示された．本研究の成果を活用することで，離散的な水分量とマトリックポテ
ンシャルの測定データから自動的に dual-VG モデルのパラメータを決定することが可能となることから，黒ボク
土等の団粒構造が発達した土壌を含め，広い範囲で土壌
水分動態の数値シミュレーションの計算精度の向上が可能となると考えられる．また，現時点において化学性のデータに比べて事例が少ない土壌の物理性データベース
作成のための統一的なパラメータ決定法としての活用が期待される．

謝辞
本研究は生研支援センター「イノベーション創出強化
研究推進事業（JPJ007097）」による委託課題（代表：柳井洋介）の支援及び JSPS 科研費 JP22K05917（代表：岩Fig. 13 UNSODA 2130（Sandy loam，アメリカ）の 2 点法と大域探索による dual-VG モデルの近似曲線．
Fitted curves of dual-VG model with the 2-point method and theglobal search method for UNSODA 2130 (Sandy loam, USA).

上に隠れている）．
このように，本研究で提案した 2 点法は大域法とほぼ等しいデータとの適合度の解が得られた．大域法では，たとえば Table 1 の条件で計算するためには 24,000 通りの初期値からの計算が必要とされることと比較すると，2点法は極めて迅速に計算が可能である．著者の PC（Macmini，Apple M1 チップ，メモリ 16GB）で本研究の試験圃場 98 個のパラメータを決定してその平均時間を計測したところ，2 点法では 0.044 秒，大域法では 46 秒となり，大域法は 2 点法の 1000 倍の計算時間がかかった．大域法は一般公開の Web サーバー上で実行するにはサーバーの負荷が大きすぎ，また計算時間がかかりす

田幸良）の助成を受けた．

引用文献
Cavanaugh, J.E. (1997): Unifying the derivations of the
Akaike and corrected Akaike information criteria．Statistics and Probability Letters, 31(2): 201–208.
Durner, W. (1994): Hydraulic conductivity estimation for soils
with heterogeneous pore structure．Water Resources Research, 30(2): 211–223.
Gauch Jr., H.G. (2002): Scientiﬁc method in practice. Cambridge University Press, pp. 269–326, Cambridge, United
Kingdom.
池内隆夫, 佃晋太郎 (2020): 枠板式高うね栽培システム（かがわ型アスパラガス栽培システム）. 農業技術体系

野菜編 8–2, pp. 基 274 の 1 の 4 の 2–基 274 の 1 の 10. 農
文協, 戸田.

ぎるところ，SWRC Fit では 2 点法を採用することで，

Murano, H., Takata, Y. and Isoi, T. (2015): Origin of the soil

Web サーバーに大きな負荷をかけることなく実装することができた．そのため，SWRC Fit を利用することにより誰でも手軽に dual-VG モデルのパラメータを精度良く

texture classiﬁcation system used in Japan. Soil Science

得ることができるようになった．

and Plant Nutrition, 61: 688–697.
Nemes, A., Schaap, M., Leij, F.J. and Wösten, J.H.M. (2015):
UNSODA 2.0: Unsaturated Soil Hydraulic Database.

土壌の物理性

44

第 155 号

Database and program for indirect methods of estimating

modal unsaturated soil hydraulic properties. Vadose Zone

unsaturated hydraulic properties. US Salinity Laboratory ARS - USDA.
農 研 機 構 (2021):(2023)

Journal, 21(1): e20168.
Seki, K., Toride, N. and van Genuchten, M.Th. (2023): Eval-

日 本 土 壌 イ ン ベ ン ト リ ー. https://uation of a general model for multimodal unsaturated soil

soil-inventory.rad.naro.go.jp/（2023 年 7 月 26 日参照）

hydraulic properties. Journal of Hydrology and Hydrome-

小原 洋, 大倉利明, 高田裕介, 神山和則, 前島勇治, 浜崎忠雄

(2011): 包括的土壌分類 第 1 次試案. 農業環境技術研究

chanics, 71(1): 22–34.
寺沢四郎 (1972): 4. 保水性. 土壌物理性測定法委員会編 土壌所報告, 29: 1–73.

物理性測定法, pp. 137–159. 養賢堂, 東京.

坂井 勝, 取出伸夫 (2007): 砂質土と団粒土に対する水分移動

van Genuchten, M.Th. (1980): A closed-form equation for

特性関数の検討. 土壌の物理性, 107：63–77.

predicting the hydraulic conductivity of unsaturated soils.

Seki, K. (2007): SWRC ﬁt - a nonlinear ﬁtting program with

Soil Science Society of America Journal, 44: 892–898.

a water retention curve for soils having unimodal and bi-

van Genuchten, M.Th., Leij, F.J. and Yates, S.R. (1991):

modal pore structure. Hydrology and Earth System Sci-

The RETC code for quantifying the hydraulic functions ofences Discussions, 4: 407–437.

unsaturated soils, version 1.0. EPA Report 600/2-91/065,

関 勝寿 (2019): 水分特性曲線の回帰プログラム SWRC Fit

U.S. Salinity Laboratory, USDA, ARS, Riverside, Califor-(3) — 回帰アルゴリズム —. 東洋大学紀要自然科学篇,
63: 29–39.
Seki, K., Toride, N. and van Genuchten, M.Th. (2022):
Closed-form hydraulic conductivity equations for multi-

要nia.
This article is an open access article under theterms of the Creative Commons Attribution (CCBY) license (https://creativecommons.org/licenses/by/4.0/).

旨土壌の水分特性曲線の近似では van Genuchten の VG モデルが広く使われているが，団粒構造が発達した黒ボク土のような土壌では，VG モデルを足し合わせる Durner（1994) の dual-VG モデルがより適している．本研究では，SWRC Fit の dual-VG モデルによる非線形回帰のアルゴリズムを改良した．すな
わち，水分特性曲線を高水分領域と低水分領域に分割し，それぞれを VG モデルで近似して得たパラメータを dual-VG モデルの初期値として与えて近似するという手法である．日本全国のアスパラガス
圃場を中心とした試験圃場の実測データによりこの手法の精度を検証した．開発された手法により，検証されたすべての土壌において大域解とほぼ等しい適合度の曲線が得られることが示された．また，修正 AIC による VG モデル，dual-VG モデル，dual-VG-CH モデル（dual-VG モデルにおいて α1 =α2 としたモデル）の比較をしたところ，黒ボク土，低地土，褐色森林土において，dual-VG モデルが最も適している試料が多いことが示された．
キーワード：水分特性曲線，非線形回帰，dual-VG モデル