Parameter Variance of the Duncan Formula for Nonlinear Shear Strength of Coarse-Grained Soil

1. Introduction

Coarse-grained soil has the advantages of easy accessibility, low cost, and high permeability and is therefore widely utilized in engineering projects such as roads, earth-rock dams, land reclamation, and high-fill airports. The shear strength of coarse-grained soil is a critical technical parameter as it directly affects both project safety and investment benefits.

The guidelines for determining the standard soil shear strength vary among different standards. For instance, the Design Code for Rolled Earth-Rock Fill Dams (SL274-2020) and standard NB/T10872-2021. stipulate that the standard value should be the average of the minimum shear strength. These specifications also define the required number of test groups: for anti-seepage soil materials, at least 11 groups are required, and each group contains 4–6 confining pressure levels. For coarse-grained materials, NB/T10872-2021 mandates a minimum of six test groups, whereas SL274-2020 requires at least 11 groups for both fine-grained soil and coarse-grained soil. The Unified Standard for Reliability Design of Hydraulic Engineering Structures (GB5099-2013) recommends using the 0.1 quantile of the probability distribution as the standard value for the strength of geotechnical materials and artificial foundations. In addition, the Code for Investigation of Geotechnical Engineering (GB 50021-2001, 2009 Edition) specifies that the standard value of a geotechnical parameter should be based on the one-sided confidence limit of the confidence interval for the mean value of the parameter matrix, which is derived using statistical interval estimation theory. By simplifying this approach, the relationship between the standard value ${\varphi }_k$ and the mean value ${\varphi }_m$ can be expressed as follows:

$${\varphi }_k={\gamma }_s{\varphi }_m,{\gamma }_s=1\pm ({\displaystyle \frac{1.704}{n^{1/2}}}+{\displaystyle \frac{4.678}{n^2}})\delta$$

(1)

where $\delta$ is the coefficient of variation, which is the ratio of the standard deviation to the mean. n is the number of groups.

Apparently, it is equally important to determine the mean and variance (or standard deviation) of the soil shear strength for engineering practices.

Currently, the shear strength of coarse-grained soil is predominantly represented by the Duncan formula. In reliability analyses of the anti-sliding stability and foundation bearing capacity of coarse-grained soil, it is essential to predefine the variance of the parameters. Therefore, it is crucial to investigate the variance of the variables in the Duncan formula for determining the nonlinear shear strength of coarse-grained soil (Aghchai et al., 2020).

2. Nonlinear Shear Strength of Coarse-Grained Soil

The shear strength of soil is generally described by the Mohr-Coulomb theory:

$${\tau }_f={\sigma }_{n}tan\phi +c\mathrm{or}{\sigma }_1={\sigma }_3{tan}^2\left(\phi /2+{45}^{°}\right)+2ctan\left(\phi /2+{45}^{°}\right),$$

(2)

where ${\tau }_f$ is the shear strength of the soil, ${\sigma }_n$ is the normal stress on the shear plane, $c$ and $\phi$ are the cohesion (or interlocking force) and internal friction angle of the soil, and ${\sigma }_1\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }_3$ are the major and minor principal stresses at the limit equilibrium failure point in the triaxial test.

Typically, the indices c and $\phi$ are constants. However, for coarse-grained soil, the shear strength exhibits strong nonlinearity; hence, it is necessary to define these indices within a limited confining pressure range. In contrast to the linear relationship observed under low confining pressure conditions, under high confining pressure conditions, the strength envelope deviates from a straight line. To characterize the shear strength of coarse-grained soil across a wide confining pressure range, the following approaches have been employed in engineering practices. The first method involves a segmented description, in which different values of c and $\phi$ are used for low (typically 100–800 kPa) and high confining pressure conditions. Another approach directly models the curved envelope of the shear strength, such as the Demeler formula ${\tau }_f=A{\sigma }^m$ and the Duncan formula

$$\mathsf{\phi }={\mathsf{\phi }}_0-∆\mathsf{\phi }\mathrm{l}\mathrm{g}\left({\mathsf{\sigma }}_3/{\mathrm{P}}_{\mathrm{a}}\right),$$

(3)

where $A$ and $m$ are the shear strength coefficient and exponent. ${\mathsf{\phi }}_0$ is the internal friction angle when the confining pressure is equal to the standard atmosphere, and $∆\mathsf{\phi }$ is the decrease in the friction angle when the confining pressure ${\mathsf{\sigma }}_3$ increases by 10 times. According to the Duncan formula, the friction angle is expressed as follows: $\mathsf{\phi }={\mathrm{s}\mathrm{i}\mathrm{n}}^{-1}\left({\displaystyle \frac{{\mathsf{\sigma }}_1-{\mathsf{\sigma }}_3}{{\mathsf{\sigma }}_1+{\mathsf{\sigma }}_3}}\right)$.

There are two methods for determining nonlinear strength indices. In the first method, the strength indices ${\mathsf{\phi }}_0$ and $∆\mathsf{\phi }$ are determined based on the test groups, and then, their minimum average values are determined. In the second method, the minimum average values of $\left({\sigma }_1-{\sigma }_3\right)$ and $\left({\sigma }_1+{\sigma }_3\right)$ under each confining pressure are determined first, and then, the friction angle is calculated based on these minimum average values. Finally, the strength indices are determined based on the Duncan formula. The parameters derived using these two methods differ and both methods suffer from an insufficient number of grouped samples, leading to large errors and a low accuracy for intermediate transition variables. To address this issue, all test data should be utilized, and the least squares method should be used to directly estimate the nonlinear strength parameters (Chen et al., 2019; Cao et al., 2007).

Linear regression for estimating Duncan’s nonlinear shear strength parameters of coarse-grained soil.

$$\mathsf{\phi }={\mathsf{\phi }}_0-∆\mathsf{\phi }\mathrm{l}\mathrm{g}\left({\displaystyle \frac{{\sigma }_3}{P_a}}\right)$$

(4)

where $\mathrm{l}\mathrm{g}\left({\mathsf{\sigma }}_3/{\mathrm{P}}_{\mathrm{a}}\right)$ is the independent variable, which is a non-random variable. $\mathsf{\Phi }$ is the dependent variable, which is a random variable. ${\mathsf{\phi }}_0$ and −$∆\mathsf{\phi }$ are the parameters to be determined.

The prerequisite of this method is that the regression residuals meet the requirements of classic least squares estimation.

3. Classic Least Squares Estimation

The linear regression model is as follows:

$${\left(\phi \right)}_i={\phi }_0-∆\phi {\left[lg\left({\sigma }_3/P_a\right)\right]}_i+e_i,(i=\mathrm{1,2},\dots ,n),$$

(5)

where $e_i$ is the residual term or the error term, which represents the influences of factors other than $\left[lg\left({\sigma }_3/P_a\right)\right]$ on the friction angle $\phi$ and the measurement error.

If the residual term $e_i(i=\mathrm{1,2},\dots ,n)$satisfies the following three conditions, then the least squares regression is valid: (a) the mathematical expectation is $E\left(e_i\right)=0$, (b) the residuals under different confining pressures have equal variances, that is, $\mathrm{V}\mathrm{a}\mathrm{r}\left(e_i\right)={\mathsf{\sigma }}^2$, ($i=\mathrm{1,2},\dots ,n$), and (c) the residuals under different confining pressures are uncorrelated, that is, $\mathrm{C}\mathrm{o}\mathrm{v}(e_i,e_j)=0$, $i\ne j$. The above is the Gauss-Markov hypothesis.

The least squares method finds estimates of ${\phi }_0$ and −$∆\phi$ that minimize the square of the length of the residual vector $e={\left(\phi \right)}_i-\left\{{\mathit{\phi }}_0-∆\mathit{\phi }{\left[\mathit{l}\mathit{g}\left({\mathit{\sigma }}_3/{\mathit{P}}_{\mathit{a}}\right)\right]}_{\mathit{i}}\right\}$, i.e.,${‖{\left(\phi \right)}_i-\left\{{\phi }_0-∆\phi {\left[lg\left({\sigma }_3/P_a\right)\right]}_i\right\}‖}^2$.

By expanding this formula and taking partial derivatives with respect to ${\phi }_0$ and −$∆\phi$ and setting them equal to zero, ${\phi }_0$ and −$∆\phi$ can be obtained.

$$\begin{gathered} {\left[lg\left({\sigma }_3/P_a\right)\right]}_i=x_i,{\left(\phi \right)}_i=y_i,{\displaystyle \frac{1}{\mathrm{n}}}\sum _{i=1}^n{x}_i=\overline{x},{\displaystyle \frac{1}{\mathrm{n}}}\sum _{i=1}^n{y}_i=\overline{y}, \\ \begin{array}{c}-∆\phi ={\displaystyle \frac{\sum _{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sum _{i=1}^n{(x_i-\overline{x})}^2}}\\ {\phi }_0=\overline{y}+∆\phi \overline{x}.\end{array}, \end{gathered}$$

(6)

By substituting n sets of observation data $\left({\left[lg\left({\sigma }_3/P_a\right)\right]}_i,{\left(\phi \right)}_i\right)$ into the above equation, we can obtain estimated values of ${\phi }_0$ and −$∆\phi$.

4. Homogeneity of Variance and Correlation Test of the Residuals of the Duncan Formula Estimated VIA Least Squares Regression

Currently, the prerequisite of the least squares method is generally not tested when determining the shear strength parameters of soil. Chen et al. (2005) noticed the variance of the regression residuals and proposed that the weighted least squares method be applied to improve the regression. The weighting coefficient is $w_i=1/{\left(\alpha +\beta x_i\right)}^2$, and $\alpha \mathrm{a}\mathrm{n}\mathrm{d}\beta$ are obtained via linear regression of the residuals and independent variables. This approach partially eliminates the variance of the residuals, but it does not fundamentally solve the problem. In the regression method recommended by the specification, i.e., $q={\displaystyle \frac{{\sigma }_1-{\sigma }_3}{2}},\mathrm{a}\mathrm{n}\mathrm{d}p={\displaystyle \frac{{\sigma }_1+{\sigma }_3}{2}}$, the independent variable $p={\displaystyle \frac{{\sigma }_1+{\sigma }_3}{2}}$ is also a random variable. Based on this characteristic, Yu et al. (2012) proposed that the orthogonal least squares method be applied to correct the shear strength estimated based on regression using $p$-$q$. The above methods only involve regressions of the parameters of the Mohr-Coulomb formula, which only partially eliminates the variance of the regression residuals. To the best of our knowledge, the variance and correlation of the residuals in nonlinear regression of the shear strength of coarse-grained soil has not been studied (Chen et al., 2007; Gong et al., 2020; Gong et al., 2018).

4.1. Regression Residuals and the Covariance Matrix

Based on the test results for several coarse-grained materials from the sites of engineering projects, the nonlinear shear strength was determined to examine the variance and correlation of the residuals in the nonlinear strength regression. Specifically, in the triaxial test on coarse-grained granite from Area I of earth-rock dam project A, twelve sets of parallel tests were conducted under six confining pressure conditions (i.e., 100, 500, 1000, 1500, 2000, and 2500 kPa). The regression equation is shown in Figure 1, and the regression residuals are shown in Figure 2. As can be seen from Figure 2, the distribution of the residuals varies across the different confining pressures, indicating that the variance of the residuals is not uniform. To determine whether the residuals under the six confining pressures are correlated, it is necessary to test the covariance matrix of the residuals. The regression residuals of the friction angle from the triaxial tests were calculated, and the residual covariance matrix was obtained:

$$\mathrm{C}\mathrm{o}\mathrm{v}\left({\mathrm{e}}_{\mathrm{i}},{\mathrm{e}}_{\mathrm{j}}\right)=\left[\begin{array}{cc}\begin{array}{ccc}3.685& 0.969& 0.331\\ 0.969& 1.147& 0.308\\ 0.331& 0.308& 0.402\end{array}& \begin{array}{ccc}1.591& 0.82& 0.688\\ 0.747& 0.438& 0.353\\ 0.641& 0.432& 0.188\end{array}\\ \begin{array}{ccc}1.591& 0.747& 0.641\\ 0.82& 0.438& 0.432\\ 0.688& 0.353& 0.188\end{array}& \begin{array}{ccc}4.02& 2.432& 0.754\\ 2.432& 1.607& 0.468\\ 0.754& 0.468& 0.286\end{array}\end{array}\right]$$

(7)

It can be seen that the main diagonal elements of the residual covariance matrix are not equal, and the non-diagonal elements are not zero, suggesting variance homogeneity and correlation of the residuals. In order to accurately determine the homogeneity and correlation of the regression residuals, hypothesis testing is required.

First, the Kolmogorov–Smirnov (K-S) method was applied to test whether the residuals as a whole satisfy the assumption of a normal distribution. Next, assuming that the ratio of the sample variances of the regression residuals at each confining pressure level follows the F-distribution, the sample variance at each confining pressure was calculated to test whether the residual variances were consistent across different confining pressures under a certain probability level (quantile). This step corresponds to the variance homogeneity test. Finally, the correlation of the regression residuals under different confining pressure levels was examined.

4.2. K-S Test for Normal Distribution of Regression Residuals

(1) Null hypothesis $H_0$: $F\left(x\right)=F_0\left(x\right);$and alternative hypothesis $H_1$: $F\left(x\right)\ne F_0\left(x\right)$, where $F_0\left(x\right)$ is the given normal distribution function. (2) The test statistic $D_n=sup\left|F_n\left(x\right)-F\left(x\right)\right|$ is selected as the principal variable. When the null hypothesis is true, the test statistic is $D_n=sup\left|F_n\left(x\right)-F_0\left(x\right)\right|$. At the significance level $\mathsf{\alpha }$, if $\sqrt{n}{D}_n>k$, then the null hypothesis is rejected, where $k$ is an appropriate constant. Kolmogorov derived the limiting distribution of $\sqrt{n}{D}_n$ and tabulated it. (3) The rejection region is determined. For a given significance level $\mathsf{\alpha }$ that satisfies $\mathrm{P}(\sqrt{n}{D}_n\ge {\mathrm{t}}_{\mathsf{\alpha }})=\mathsf{\alpha }$, the quantile ${\mathrm{t}}_{\mathsf{\alpha }}$ is the critical value, and the rejection region is $\left[t_{\alpha },+\infty )\right.$.

Table 1 presents the regression residuals for the coarse-grained granite from Area I of project A under various confining pressure conditions. The mean and standard deviation of all of the residuals are 2.6659 $\times$10⁻⁴ and 1.6573. For a significance level of 0.01, $P\left(\sqrt{n}{D}_n\ge t_{0.01}\right)=$0.01, and t $=$1.63, which is taken as an approximate value of $t_{0.01}$. Then, the rejection region is $[1.63,+\mathrm{\infty }$). $D_n=max(d_{01},d_{02})=$ 0.1782. Because $\sqrt{n}{D}_n=\sqrt{72}\times$ 0.1782 = 1.5119, and it falls outside the rejection region [1.63, +$\mathrm{\infty }$), we accept the hypothesis that the regression residuals of the shear strength are normally distributed.

4.3. Homogeneity Test of Regression Residuals Under Different Confining Pressure Conditions

Theoretical studies have shown that the F-test is sensitive to deviations from homogeneity of variance, making it suitable for testing the variance of residuals under different confining pressures. The null hypothesis H₀ and alternative hypothesis H₁ are $H_0:{\sigma }_1^2={\sigma }_2^2=\dots ={\sigma }_k^2$ and $H_1:{\sigma }_1^2$, ${\sigma }_2^2$, $\dots$, ${\sigma }_k^2$, and they are not all equal.

Determine the rejection region.

Let $\mathsf{\alpha }=0.05$,

$$P(F\le F_{1-{\displaystyle \frac{\alpha }{2}}}(n_1-1,n_2-1\left)\right)={\displaystyle \frac{\alpha }{2}},P\left(F\ge F_{{\displaystyle \frac{\alpha }{2}}}\left(n_1-1,n_2-1\right)\right)={\displaystyle \frac{\alpha }{2}},$$

(8)

$$P({\displaystyle \frac{S_1^2}{S_2^2}}\le F_{0.975}(\mathrm{11,11}\left)\right)=0.025,$$

$$F_{0.975}\left(\mathrm{11,11}\right)={\displaystyle \frac{1}{F_{0.025}\left(\mathrm{11,11}\right)}}={\displaystyle \frac{1}{3.48}}=0.287356,$$

and $F_{0.025}\left(\mathrm{11,11}\right)=3.48$.

Then, the rejection region is (0, 0.2874] and [3.48, +$\mathrm{\infty })$.

Table 2 shows the sample variance ratios of the regression residuals under various confining pressures and whether they are in the rejection region. $\surd$ indicates that the variance ratio is in the acceptance region and the null hypothesis is accepted, while $\times$indicates that the variance ratio is in the rejection region and the null hypothesis is not accepted.

It can be seen that except for the diagonal elements where the variance ratio is equal to 1, seven of the sample variance ratios of the residuals under different confining pressures are in the acceptance region and eight are in the rejection region. Therefore, the variance homogeneity hypothesis is not accepted.

4.4. Residual Independence Test Under Different Confining Pressure Conditions

Assuming that the regression residuals under different confining pressures are two different random variables U and V, the independence test is equivalent to testing their correlation coefficient $\mathsf{\rho }$. The null hypothesis H₀ and alternative hypothesis H₁ are $H_0:\rho =0,\mathrm{a}\mathrm{n}\mathrm{d}{H}_1:\rho \ne$ 0. The sample correlation coefficient R is used as the point estimator of $\rho$ . R is shown in Table 3.

$$R={\displaystyle \frac{\sum _{i=1}^n(U_i-\overline{U})(V_i-\overline{V})}{\sqrt{\sum _{i=1}^n{(U_i-\overline{U})}^2}\sqrt{\sum _{i=1}^n{(V_i-\overline{V})}^2}}}.$$

(9)

When $H_0$ is true,

$$t={\displaystyle \frac{R}{\sqrt{1-R^2}}}\sqrt{n-2}.$$

(10)

It follows a t-distribution with degrees of freedom of $n-2$. t can be found in Table 4. Given a significance level $\mathsf{\alpha }$,

$$P\left(\left|t\right|\ge t_{{\displaystyle \frac{\alpha }{2}}}\left(n-2\right)\right)=\alpha .$$

(11)

If the observed value of t satisfies $\left|t\right|\ge t_{{\displaystyle \frac{\alpha }{2}}}\left(n-2\right)$, then the null hypothesis $H_0$ is rejected.

$$\mathsf{\alpha }=0.1.$$

$\mathrm{P}(\left|\mathrm{t}\right|\ge {\mathrm{t}}_{0.05}(10\left)\right)=0.1$, and ${\mathrm{t}}_{0.05}\left(10\right)$=1.812.

The rejection regions are (0, −1.812] and [1.812, +$\mathrm{\infty })$.

It can be seen that seven of the above t values are within the rejection region. Therefore, the independence hypothesis of the residuals under each confining pressure is rejected.

The regression residual analysis of the nonlinear shear strength based on the Duncan formula for the coarse-grained soils from three engineering projects (projects A, B, and C) indicates that the diagonal elements of the covariance matrix of the least squares regression residuals for the nonlinear shear strength friction angle are not equal, and the non-diagonal elements are not zero. Although the probability distribution of the residuals passes the normality test, the residuals under different confining pressures exhibit homogeneity of variance and irrelevance. This confirms the occurrence of residual variance heteroscedasticity and autocorrelation across the different confining pressures.

5. Elimination of Residual Variance Heteroscedasticity and Autocorrelation

If the residuals exhibit heteroscedasticity and autocorrelation, the regression equation no longer meets the prerequisites of least squares regression. In this case, if the least squares method is used, although the estimated values of the parameters are unbiased, the variances of the parameters could deviate from the actual value. A solution is to replace the regression variables. For the linear regression equation

$$\mathit{y}=\mathit{X}\mathit{\beta }+\mathit{e},\mathit{E}\left(\mathit{e}\right)=0,\mathit{C}\mathit{o}\mathit{v}\left(\mathit{e}\right)={\mathit{\sigma }}^2\mathit{V},$$

(12)

where $\mathrm{X}=\left[\begin{array}{ccc}1& 1& \begin{array}{cc}\dots & 1\end{array}\\ {\mathrm{x}}_1& {\mathrm{x}}_2& \begin{array}{cc}\dots & {\mathrm{x}}_{\mathrm{n}}\end{array}\end{array}\right]\mathrm{\text{'}}$, $\mathrm{y}=\left[\begin{array}{ccc}\begin{array}{cc}{\mathsf{\phi }}_1& {\mathsf{\phi }}_2\end{array}& \dots & {\mathsf{\phi }}_{\mathrm{n}}\end{array}\right]$‘. $\mathit{\beta }=\left[\begin{array}{cc}{\mathit{\phi }}_0& -∆\mathit{\phi }\end{array}\right],{\mathrm{x}}_1=\mathrm{l}\mathrm{g}\left({\mathsf{\sigma }}_3^1/{\mathrm{P}}_{\mathrm{a}}\right),\dots ,{\mathrm{x}}_{\mathrm{i}}=\mathrm{l}\mathrm{g}\left({\displaystyle \frac{{\sigma }_3^i}{P_a}}\right),\dots ,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{x}}_{\mathrm{n}}=\mathrm{l}\mathrm{g}\left({\mathsf{\sigma }}_3^{\mathrm{n}}/{\mathrm{P}}_{\mathrm{a}}\right)$, ${\mathsf{\sigma }}_3^{\mathrm{i}}$ and ${\mathsf{\phi }}_{\mathrm{i}}$ are the confining pressure and friction angle of the $\mathrm{i}$th test. The symbol ‘ denotes matrix transposition. $V$ is the covariance matrix of the residual vector under different confining pressures. ${\mathsf{\sigma }}^2$ is the variance of the residual vector.

It can be seen that the covariance matrix of the residual is not an identity matrix. $V$ is nonsingular and positive definite, and there exists a nonsingular symmetric matrix $K$ such that $K‘K=K^{-1}K=V$, that is, the matrix $K$ is the square root matrix of $V$.

By defining new variables (Wang et al., 2021; Douglas et al., 2022), namely,

$$\mathit{z}=K^{-1}y,B=K^{-1}X,\mathrm{and}g=K^{-1}e$$

(13)

the regression equation $\mathit{y}=\mathit{X}\mathit{\beta }+\mathit{\epsilon }$ becomes ${\mathit{K}}^{-1}\mathit{y}={\mathit{K}}^{-1}\mathit{X}\mathit{\beta }+{\mathit{K}}^{-1}\mathit{e}$, i.e.,

$$\mathit{z}=\mathit{B}\mathit{\beta }+\mathit{g}.$$

(14)

After this transformation, the model residuals have a zero mean, i.e., $\mathit{E}\left(\mathit{g}\right)={\mathit{K}}^{-1}\mathit{E}\left(\mathit{e}\right)$=0, and the covariance matrix of $\mathit{g}$ is

$$\begin{gathered} \mathit{V}\mathit{a}\mathit{r}\left(\mathit{g}\right)=\left\{\left[\mathit{g}-\mathit{E}\left(\mathit{g}\right)\right]{\left[\mathit{g}-\mathit{E}\left(\mathit{g}\right)\right]}^{\mathrm{\text{'}}}\right\}=\mathit{E}\left(\mathit{g}{\mathit{g}}^{\mathrm{\text{'}}}\right) \\ =\mathit{E}\left({\mathit{K}}^{-1}\mathit{e}{\mathit{e}}^{\mathrm{\text{'}}}{\mathit{K}}^{-1}\right)={\mathit{K}}^{-1}\mathit{E}\left(\mathit{e}{\mathit{e}}^{\mathrm{\text{'}}}\right){\mathit{K}}^{-1} \\ ={\sigma }^2{\mathit{K}}^{-1}V{\mathit{K}}^{-1}={\sigma }^2{\mathit{K}}^{-1}\mathit{K}\mathit{K}{\mathit{K}}^{-1}={\sigma }^{2}I. \end{gathered}$$

Therefore, $\mathit{g}$ has a zero mean, constant variance, and uncorrelated variance. Equation (14) satisfies the general assumptions of the classic least squares estimation, and the least squares function is

$$\mathit{Q}\left(\mathit{\beta }\right)={\mathit{g}}^{\mathrm{\text{'}}}\mathit{g}={\mathit{e}}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\mathit{e}={\left(\mathit{y}-\mathit{X}\mathit{\beta }\right)}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\left(\mathit{y}-\mathit{X}\mathit{\beta }\right).$$

The canonical equation for least squares is

$${\mathit{X}}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\mathit{X}{\widehat{\mathit{\beta }}}^{\mathrm{*}}={\mathit{X}}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\mathit{y}.$$

The solution to the equation is

$${\widehat{\mathit{\beta }}}^{\mathrm{*}}={\left({\mathit{X}}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\mathit{X}\right)}^{-1}{\mathit{X}}^{\mathrm{\text{'}}}{\mathit{V}}^{-1}\mathit{y}.$$

(15)

6. Examples

The triaxial test data for five types of coarse-grained materials used in three engineering projects (A, B, and C) were used to perform classic least squares regression and regression after eliminating the heteroscedasticity and autocorrelation of the residual variance. After obtaining the regression equation, the covariance matrix of the parameters was calculated (Table 5).

It can be seen that for each coarse-grained material, the regression parameters obtained using the two methods are similar, yet the variances of the parameters are reduced when the proposed method is applied.

For example, for the coarse-grained material from Area II of Project A, the regression equation from the classic least squares method and the covariance matrix between ${\phi }_0\mathrm{a}\mathrm{n}\mathrm{d}∆\mathsf{\phi }$ are

$$\mathsf{\phi }=50.309-11.680\mathrm{l}\mathrm{g}\left({\displaystyle \frac{{\sigma }_3}{P_a}}\right)$$

(16)

$$\mathrm{and}\mathrm{C}\mathrm{o}\mathrm{v}\left({\phi }_0,∆\mathsf{\phi }\right)=\left[\begin{array}{cc}1.854& -1.763\\ -1.763& 2.265\end{array}\right]$$

(17)

After eliminating the residual variance heteroscedasticity and autocorrelation, the regression equation and covariance matrix are

$$\mathsf{\phi }=50.709-11.721\mathrm{l}\mathrm{g}\left({\mathsf{\sigma }}_3/{\mathrm{P}}_{\mathrm{a}}\right)$$

(18)

$$\mathrm{and}\mathrm{C}\mathrm{o}\mathrm{v}\left({\phi }_0,∆\mathsf{\phi }\right)=\left[\begin{array}{cc}0.569& -0.312\\ -0.312& 0.309\end{array}\right]$$

(19)

It can be seen that after eliminating the heteroscedasticity and autocorrelation, the parameter variance is reduced, providing more accurate estimates of the shear strength parameters of the coarse-grained soil.

7. Conclusions

The shear strength of coarse-grained soil has strong non-linearity and is commonly determined using the Duncan formula. However, the classic grouping method for estimating the parameters has the limitations of small sample sizes and large errors. When the classic least squares method is applied to estimate the Duncan shear strength parameters based on triaxial test data, it has been found that the residual variance exhibited heteroscedasticity and autocorrelation, which violates the prerequisites of the least squares method. In this study, we introduced a generalized least squares method to eliminate these issues. By analyzing triaxial test data from five types of coarse-grained soil materials used in three engineering projects, it was found that while the mean values of the shear strength parameters remained relatively unchanged compared to those obtained using the classic least squares method, the variance was significantly reduced, demonstrating the effectiveness of the proposed method. These improvements enable a more accurate characterization of the nonlinear shear strength of coarse-grained soil materials, enhancing the precision of geotechnical calculations.