Thank you for your contributions.I have a question which I hope you can help me to answer
   The dual_channel_image.shape I saved with sitk is (275, 386, 386, 386, 2), why the shape I read with nibabel is (386, 386, 275, 1, 2), the result is one dimension more? 

        dual_channel_image = np.stack((image_R, image_L), axis=-1)
        heatmap_sitk = sitk.GetImageFromArray(dual_channel_image)
        output_file = f"/mnt/nas/temp/sch/landmark/{data[i]}/{data[i]}_666channel.nii.gz"
        heatmap_sitk.SetSpacing(spacing)
        heatmap_sitk.SetOrigin(Origin)
        sitk.WriteImage(heatmap_sitk, output_file)
        img = nib.load(output_file)
        data_array = img.get_fdata()
        data_shape = img.shape
        print("Data shape:", data_shape)
        print("OKOKOK")